Context Navigation

source: trunk/source/processes/hadronic/models/parton_string/diffraction/src/G4DiffractiveExcitation.cc @ 847

Last change on this file since 847 was 819, checked in by garnier, 16 years ago
import all except CVS
File size: 21.9 KB

Rev	Line
[819]	1	//
	2	// ********************************************************************
	3	// * License and Disclaimer *
	4	// * *
	5	// * The Geant4 software is copyright of the Copyright Holders of *
	6	// * the Geant4 Collaboration. It is provided under the terms and *
	7	// * conditions of the Geant4 Software License, included in the file *
	8	// * LICENSE and available at http://cern.ch/geant4/license . These *
	9	// * include a list of copyright holders. *
	10	// * *
	11	// * Neither the authors of this software system, nor their employing *
	12	// * institutes,nor the agencies providing financial support for this *
	13	// * work make any representation or warranty, express or implied, *
	14	// * regarding this software system or assume any liability for its *
	15	// * use. Please see the license in the file LICENSE and URL above *
	16	// * for the full disclaimer and the limitation of liability. *
	17	// * *
	18	// * This code implementation is the result of the scientific and *
	19	// * technical work of the GEANT4 collaboration. *
	20	// * By using, copying, modifying or distributing the software (or *
	21	// * any work based on the software) you agree to acknowledge its *
	22	// * use in resulting scientific publications, and indicate your *
	23	// * acceptance of all terms of the Geant4 Software license. *
	24	// ********************************************************************
	25	//
	26	//
	27	// $Id: G4DiffractiveExcitation.cc,v 1.1 2007/05/25 06:56:53 vuzhinsk Exp $
	28	// ------------------------------------------------------------
	29	// GEANT 4 class implemetation file
	30	//
	31	// ---------------- G4DiffractiveExcitation --------------
	32	// by Gunter Folger, October 1998.
	33	// diffractive Excitation used by strings models
	34	// Take a projectile and a target
	35	// excite the projectile and target
	36	// Essential changed by V. Uzhinsky in November - December 2006
	37	// in order to put it in a correspondence with original FRITIOF
	38	// model. Variant of FRITIOF with nucleon de-excitation is implemented.
	39	// Other changes by V.Uzhinsky in May 2007 were introduced to fit
	40	// meson-nucleon interactions
	41	// ---------------------------------------------------------------------
	42
	43
	44	#include "globals.hh"
	45	#include "Randomize.hh"
	46
	47	#include "G4DiffractiveExcitation.hh"
	48	#include "G4LorentzRotation.hh"
	49	#include "G4ThreeVector.hh"
	50	#include "G4ParticleDefinition.hh"
	51	#include "G4VSplitableHadron.hh"
	52	#include "G4ExcitedString.hh"
	53	//#include "G4ios.hh"
	54
	55	G4DiffractiveExcitation::G4DiffractiveExcitation() // Uzhi
	56	{
	57	}
	58
	59	G4bool G4DiffractiveExcitation::
	60	ExciteParticipants(G4VSplitableHadron projectile, G4VSplitableHadron target) const
	61	{
	62
	63	G4LorentzVector Pprojectile=projectile->Get4Momentum();
	64
	65	// -------------------- Projectile parameters -----------------------------------
	66	G4bool PutOnMassShell=0;
	67
	68	// With de-excitation
	69	// G4double M0projectile=projectile->GetDefinition()->GetPDGMass();
	70	// Without de-excitation
	71	G4double M0projectile = Pprojectile.mag();
	72
	73	if(M0projectile < projectile->GetDefinition()->GetPDGMass())
	74	{
	75	PutOnMassShell=1;
	76	M0projectile=projectile->GetDefinition()->GetPDGMass();
	77	}
	78
	79	G4double Mprojectile2 = M0projectile * M0projectile;
	80
	81	G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
	82	G4int absPDGcode=std::abs(PDGcode);
	83	G4double ProjectileDiffCut;
	84	G4double ProjectileBackgroundWeight; // Uzhi May 2007
	85
	86	G4double TargetBackgroundWeight; // Uzhi May 2007
	87
	88	G4double AveragePt2;
	89
	90	if( absPDGcode > 1000 ) //------Projectile is baryon --------
	91	{
	92	ProjectileDiffCut = 1.1; // GeV
	93	ProjectileBackgroundWeight=0.; // Uzhi May 2007
	94	AveragePt2 = 0.3; // GeV^2
	95	}
	96	else if( absPDGcode == 211 \|\| PDGcode == 111) //------Projectile is Pion -----------
	97	{
	98	ProjectileDiffCut = 0.6; // GeV Uzhi May 2007 1.0->0.6
	99	ProjectileBackgroundWeight=0.9; // Uzhi May 2007
	100	AveragePt2 = 0.3; // GeV^2
	101	}
	102	else if( absPDGcode == 321 \|\| PDGcode == -311) //------Projectile is Kaon -----------
	103	{
	104	ProjectileDiffCut = 0.7; // GeV 1.1
	105	ProjectileBackgroundWeight=0.5; // Uzhi May 2007 ???
	106	AveragePt2 = 0.3; // GeV^2
	107	}
	108	else //------Projectile is undefined,
	109	//------Nucleon assumed
	110	{
	111	ProjectileDiffCut = 1.1; // GeV
	112	ProjectileBackgroundWeight=0.; // Uzhi May 2007
	113	AveragePt2 = 0.3; // GeV^2
	114	};
	115
	116	ProjectileDiffCut = ProjectileDiffCut * GeV;
	117	AveragePt2 = AveragePt2 * GeV*GeV;
	118
	119	// -------------------- Target parameters ----------------------------------------------
	120	G4LorentzVector Ptarget=target->Get4Momentum();
	121
	122	G4double M0target = Ptarget.mag();
	123
	124	if(M0target < target->GetDefinition()->GetPDGMass())
	125	{
	126	PutOnMassShell=1;
	127	M0target=target->GetDefinition()->GetPDGMass();
	128	}
	129
	130	G4double Mtarget2 = M0target * M0target; //Ptarget.mag2();
	131	// for AA-inter.
	132
	133	G4double NuclearNucleonDiffCut = 1.1*GeV;
	134	TargetBackgroundWeight=0.; // Uzhi May 2007
	135
	136	G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
	137	G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
	138
	139	// Transform momenta to cms and then rotate parallel to z axis;
	140
	141	G4LorentzVector Psum;
	142	Psum=Pprojectile+Ptarget;
	143
	144	G4LorentzRotation toCms(-1*Psum.boostVector());
	145
	146	G4LorentzVector Ptmp=toCms*Pprojectile;
	147
	148	if ( Ptmp.pz() <= 0. )
	149	{
	150	// "String" moving backwards in CMS, abort collision !!
	151	//G4cout << " abort Collision!! " << G4endl;
	152	return false;
	153	}
	154
	155	toCms.rotateZ(-1*Ptmp.phi());
	156	toCms.rotateY(-1*Ptmp.theta());
	157
	158	G4LorentzRotation toLab(toCms.inverse());
	159
	160	Pprojectile.transform(toCms);
	161	Ptarget.transform(toCms);
	162
	163	G4double Pt2;
	164	G4double ProjMassT2, ProjMassT;
	165	G4double TargMassT2, TargMassT;
	166	G4double PZcms2, PZcms;
	167	G4double PMinusNew, TPlusNew;
	168
	169	G4double S=Psum.mag2();
	170	G4double SqrtS=std::sqrt(S);
	171
	172	if(absPDGcode > 1000 && SqrtS < 2200*MeV)
	173	{return false;} // The model cannot work for
	174	// p+p-interactions
	175	// at Plab < 1.3 GeV/c.
	176
	177	if(( absPDGcode == 211 \|\| PDGcode == 111) && SqrtS < 1600*MeV)
	178	{return false;} // The model cannot work for
	179	// Pi+p-interactions
	180	// at Plab < 1. GeV/c.
	181
	182	if(( absPDGcode == 321 \|\| PDGcode == -311) && SqrtS < 1600*MeV)
	183	{return false;} // The model cannot work for
	184	// K+p-interactions
	185	// at Plab < ??? GeV/c. ???
	186
	187	PZcms2=(SS+Mprojectile2Mprojectile2+Mtarget2*Mtarget2-
	188	2SMprojectile2-2SMtarget2-2Mprojectile2Mtarget2)/4./S;
	189	if(PZcms2 < 0)
	190	{return false;} // It can be in an interaction with off-shell nuclear nucleon
	191
	192	PZcms = std::sqrt(PZcms2);
	193
	194	if(PutOnMassShell)
	195	{
	196	if(Pprojectile.z() > 0.)
	197	{
	198	Pprojectile.setPz( PZcms);
	199	Ptarget.setPz( -PZcms);
	200	}
	201	else
	202	{
	203	Pprojectile.setPz(-PZcms);
	204	Ptarget.setPz( PZcms);
	205	};
	206
	207	Pprojectile.setE(std::sqrt(Mprojectile2+
	208	Pprojectile.x()*Pprojectile.x()+
	209	Pprojectile.y()*Pprojectile.y()+
	210	PZcms2));
	211	Ptarget.setE(std::sqrt( Mtarget2 +
	212	Ptarget.x()*Ptarget.x()+
	213	Ptarget.y()*Ptarget.y()+
	214	PZcms2));
	215	}
	216
	217	G4double maxPtSquare = PZcms2;
	218
	219	//G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
	220	//G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
	221	// G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
	222
	223	// G4cout << " Projectile Xplus / Xminus : " <<
	224	// Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
	225	// G4cout << " Target Xplus / Xminus : " <<
	226	// Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
	227
	228	G4LorentzVector Qmomentum;
	229	G4double Qminus, Qplus;
	230
	231	G4int whilecount=0;
	232	do {
	233	// Generate pt
	234
	235	if (whilecount++ >= 500 && (whilecount%100)==0)
	236	// G4cout << "G4DiffractiveExcitation::ExciteParticipants possibly looping"
	237	// << ", loop count/ maxPtSquare : "
	238	// << whilecount << " / " << maxPtSquare << G4endl;
	239	if (whilecount > 1000 )
	240	{
	241	Qmomentum=G4LorentzVector(0.,0.,0.,0.);
	242	return false; // Ignore this interaction
	243	}
	244
	245	Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
	246
	247	//G4cout << "generated Pt " << Qmomentum << G4endl;
	248	//G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
	249	//G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
	250
	251	// Momentum transfer
	252	/* // Uzhi
	253	G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
	254	G4double Xmax=1.;
	255	G4double Xplus =ChooseX(Xmin,Xmax);
	256	G4double Xminus=ChooseX(Xmin,Xmax);
	257
	258	// G4cout << " X-plus " << Xplus << G4endl;
	259	// G4cout << " X-minus " << Xminus << G4endl;
	260
	261	G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
	262	G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
	263	G4double Qminus= pt2 / Xplus /Pprojectile.plus();
	264	/ // Uzhi
	265
	266	Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
	267	ProjMassT2=Mprojectile2+Pt2;
	268	ProjMassT =std::sqrt(ProjMassT2);
	269
	270	TargMassT2=Mtarget2+Pt2;
	271	TargMassT =std::sqrt(TargMassT2);
	272
	273	PZcms2=(SS+ProjMassT2ProjMassT2+
	274	TargMassT2*TargMassT2-
	275	2.SProjMassT2-2.STargMassT2-
	276	2.ProjMassT2TargMassT2)/4./S;
	277	if(PZcms2 < 0 ) {PZcms2=0;};
	278	PZcms =std::sqrt(PZcms2);
	279
	280	G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
	281	G4double PMinusMax=SqrtS-TargMassT;
	282
	283	if(G4UniformRand() < ProjectileBackgroundWeight) // Uzhi May 2007
	284	{
	285	PMinusNew=PMinusMax-(PMinusMax-PMinusMin)*G4UniformRand(); // Uzhi May 2007
	286	} // Uzhi May 2007
	287	else // Uzhi May 2007
	288	{ // Uzhi May 2007
	289	PMinusNew=ChooseP(PMinusMin, PMinusMax); // Uzhi May 2007
	290	}; // Uzhi May 2007
	291	// PMinusNew=ChooseP(PMinusMin, PMinusMax); // Uzhi May 2007
	292	Qminus=PMinusNew-Pprojectile.minus();
	293
	294	G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
	295	G4double TPlusMax=SqrtS-PMinusNew; // Uzhi May 2007
	296	// G4double TPlusMax=SqrtS-ProjMassT; // Uzhi May 2007
	297
	298	if(G4UniformRand() < TargetBackgroundWeight) // Uzhi May 2007
	299	{
	300	TPlusNew=TPlusMax-(TPlusMax-TPlusMin)*G4UniformRand(); // Uzhi May 2007
	301	} // Uzhi May 2007
	302	else // Uzhi May 2007
	303	{ // Uzhi May 2007
	304	TPlusNew=ChooseP(TPlusMin, TPlusMax); // Uzhi May 2007
	305	}; // Uzhi May 2007
	306
	307	// TPlusNew=ChooseP(TPlusMin, TPlusMax); // Uzhi May 2007
	308	Qplus=-(TPlusNew-Ptarget.plus());
	309
	310	Qmomentum.setPz( (Qplus-Qminus)/2 );
	311	Qmomentum.setE( (Qplus+Qminus)/2 );
	312
	313	//G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
	314	// G4cout << "pt2" << pt2 << G4endl;
	315	// G4cout << "Qmomentum " << Qmomentum << G4endl;
	316	// G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
	317	// " / " << (Ptarget-Qmomentum).mag() << G4endl;
	318
	319	} while (
	320	( (Pprojectile+Qmomentum).mag2() < Mprojectile2 \|\| // Uzhi No without excitation
	321	(Ptarget -Qmomentum).mag2() < Mtarget2 ) \|\| // Uzhi
	322	( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction
	323	(Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi
	324
	325	if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi
	326	//Projectile diffraction
	327	{
	328	G4double TMinusNew=SqrtS-PMinusNew;
	329	Qminus=Ptarget.minus()-TMinusNew;
	330	TPlusNew=TargMassT2/TMinusNew;
	331	Qplus=Ptarget.plus()-TPlusNew;
	332
	333	Qmomentum.setPz( (Qplus-Qminus)/2 );
	334	Qmomentum.setE( (Qplus+Qminus)/2 );
	335	}
	336	else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi
	337	//Target diffraction
	338	{
	339	G4double PPlusNew=SqrtS-TPlusNew;
	340	Qplus=PPlusNew-Pprojectile.plus();
	341	PMinusNew=ProjMassT2/PPlusNew;
	342	Qminus=PMinusNew-Pprojectile.minus();
	343
	344	Qmomentum.setPz( (Qplus-Qminus)/2 );
	345	Qmomentum.setE( (Qplus+Qminus)/2 );
	346	};
	347
	348	Pprojectile += Qmomentum;
	349	Ptarget -= Qmomentum;
	350
	351	/*
	352	Pprojectile.setPz(0.);
	353	Pprojectile.setE(SqrtS-M0target);
	354
	355	Ptarget.setPz(0.);
	356	Ptarget.setE(M0target);
	357	*/
	358
	359	//G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
	360	//G4cout << "Ptarget with Q : " << Ptarget << G4endl;
	361
	362	// G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
	363	// G4cout << "Target back: " << toLab * Ptarget << G4endl;
	364
	365	// Transform back and update SplitableHadron Participant.
	366	Pprojectile.transform(toLab);
	367	Ptarget.transform(toLab);
	368
	369	//G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
	370	//G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
	371
	372	//G4cout << "Target mass " << Ptarget.mag() << G4endl;
	373
	374	target->Set4Momentum(Ptarget);
	375
	376	//G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
	377
	378	projectile->Set4Momentum(Pprojectile);
	379
	380	return true;
	381	}
	382
	383
	384	G4ExcitedString * G4DiffractiveExcitation::
	385	String(G4VSplitableHadron * hadron, G4bool isProjectile) const
	386	{
	387	hadron->SplitUp();
	388	G4Parton *start= hadron->GetNextParton();
	389	if ( start==NULL)
	390	{ G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
	391	return NULL;
	392	}
	393	G4Parton *end = hadron->GetNextParton();
	394	if ( end==NULL)
	395	{ G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
	396	return NULL;
	397	}
	398
	399	G4ExcitedString * string;
	400	if ( isProjectile )
	401	{
	402	string= new G4ExcitedString(end,start, +1);
	403	} else {
	404	string= new G4ExcitedString(start,end, -1);
	405	}
	406
	407	string->SetPosition(hadron->GetPosition());
	408
	409	// momenta of string ends
	410	G4double ptSquared= hadron->Get4Momentum().perp2();
	411	G4double transverseMassSquared= hadron->Get4Momentum().plus()
	412	* hadron->Get4Momentum().minus();
	413
	414
	415	G4double maxAvailMomentumSquared=
	416	sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
	417
	418	G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ??????????????????
	419	G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
	420
	421	G4LorentzVector Pstart(G4LorentzVector(pt,0.));
	422	G4LorentzVector Pend;
	423	Pend.setPx(hadron->Get4Momentum().px() - pt.x());
	424	Pend.setPy(hadron->Get4Momentum().py() - pt.y());
	425
	426	G4double tm1=hadron->Get4Momentum().minus() +
	427	( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
	428
	429	G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
	430	4. * Pend.perp2() * hadron->Get4Momentum().minus()
	431	/ hadron->Get4Momentum().plus() ));
	432
	433	G4int Sign= isProjectile ? -1 : 1;
	434
	435	G4double endMinus = 0.5 * (tm1 + Sign*tm2);
	436	G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
	437
	438	G4double startPlus= Pstart.perp2() / startMinus;
	439	G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
	440
	441	Pstart.setPz(0.5*(startPlus - startMinus));
	442	Pstart.setE(0.5*(startPlus + startMinus));
	443
	444	Pend.setPz(0.5*(endPlus - endMinus));
	445	Pend.setE(0.5*(endPlus + endMinus));
	446
	447	start->Set4Momentum(Pstart);
	448	end->Set4Momentum(Pend);
	449
	450	#ifdef G4_FTFDEBUG
	451	G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
	452	G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
	453	G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
	454	G4cout << " sum of ends " << Pstart+Pend << G4endl;
	455	G4cout << " Original " << hadron->Get4Momentum() << G4endl;
	456	#endif
	457
	458	return string;
	459	}
	460
	461
	462	// --------- private methods ----------------------
	463
	464	G4double G4DiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
	465	{
	466	// choose an x between Xmin and Xmax with P(x) ~ 1/x
	467
	468	// to be improved...
	469
	470	G4double range=Pmax-Pmin; // Uzhi
	471
	472	if ( Pmin <= 0. \|\| range <=0. )
	473	{
	474	G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
	475	throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation::ChooseP : Invalid arguments ");
	476	}
	477
	478	G4double P;
	479	/* // Uzhi
	480	do {
	481	x=Xmin + G4UniformRand() * range;
	482	} while ( Xmin/x < G4UniformRand() );
	483	*/ // Uzhi
	484
	485	P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); // Uzhi
	486
	487	//debug-hpw cout << "DiffractiveX "<<x<<G4endl;
	488	return P;
	489	}
	490
	491	G4ThreeVector G4DiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
	492	{ // @@ this method is used in FTFModel as well. Should go somewhere common!
	493
	494	G4double Pt2;
	495	/* // Uzhi
	496	do {
	497	pt2=widthSquare * std::log( G4UniformRand() );
	498	} while ( pt2 > maxPtSquare);
	499	*/ // Uzhi
	500
	501	Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
	502
	503	G4double Pt=std::sqrt(Pt2);
	504
	505	G4double phi=G4UniformRand() * twopi;
	506
	507	return G4ThreeVector (Ptstd::cos(phi), Ptstd::sin(phi), 0.);
	508	}
	509
	510	G4DiffractiveExcitation::G4DiffractiveExcitation(const G4DiffractiveExcitation &)
	511	{
	512	throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation copy contructor not meant to be called");
	513	}
	514
	515
	516	G4DiffractiveExcitation::~G4DiffractiveExcitation()
	517	{
	518	}
	519
	520
	521	const G4DiffractiveExcitation & G4DiffractiveExcitation::operator=(const G4DiffractiveExcitation &)
	522	{
	523	throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation = operator meant to be called");
	524	return *this;
	525	}
	526
	527
	528	int G4DiffractiveExcitation::operator==(const G4DiffractiveExcitation &) const
	529	{
	530	throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation == operator meant to be called");
	531	return false;
	532	}
	533
	534	int G4DiffractiveExcitation::operator!=(const G4DiffractiveExcitation &) const
	535	{
	536	throw G4HadronicException(__FILE__, __LINE__, "G4DiffractiveExcitation != operator meant to be called");
	537	return true;
	538	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: