Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

G4QGSDiffractiveExcitation.cc @ 1198

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

Line
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: G4QGSDiffractiveExcitation.cc,v 1.1 2007/05/25 07:30:47 gunter Exp $
28	// ------------------------------------------------------------
29	// GEANT 4 class implemetation file
30	//
31	// ---------------- G4QGSDiffractiveExcitation --------------
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	// ---------------------------------------------------------------------
40
41	// Modified:
42	// 25-05-07 : G.Folger
43	// move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation
44	//
45
46	#include "globals.hh"
47	#include "Randomize.hh"
48
49	#include "G4QGSDiffractiveExcitation.hh"
50	#include "G4LorentzRotation.hh"
51	#include "G4ThreeVector.hh"
52	#include "G4ParticleDefinition.hh"
53	#include "G4VSplitableHadron.hh"
54	#include "G4ExcitedString.hh"
55	//#include "G4ios.hh"
56
57	G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation() // Uzhi
58	{
59	}
60
61	G4bool G4QGSDiffractiveExcitation::
62	ExciteParticipants(G4VSplitableHadron projectile, G4VSplitableHadron target) const
63	{
64
65	G4LorentzVector Pprojectile=projectile->Get4Momentum();
66
67	// -------------------- Projectile parameters -----------------------------------
68	G4bool PutOnMassShell=0;
69
70	// G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
71	G4double M0projectile = Pprojectile.mag(); // Without de-excitation
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 AveragePt2;
85
86	if( absPDGcode > 1000 ) //------Projectile is baryon --------
87	{
88	ProjectileDiffCut = 1.1; // GeV
89	AveragePt2 = 0.3; // GeV^2
90	}
91	else if( absPDGcode == 211 \|\| PDGcode == 111) //------Projectile is Pion -----------
92	{
93	ProjectileDiffCut = 1.0; // GeV
94	AveragePt2 = 0.3; // GeV^2
95	}
96	else if( absPDGcode == 321 \|\| PDGcode == -311) //------Projectile is Kaon -----------
97	{
98	ProjectileDiffCut = 1.1; // GeV
99	AveragePt2 = 0.3; // GeV^2
100	}
101	else //------Projectile is undefined, Nucleon assumed
102	{
103	ProjectileDiffCut = 1.1; // GeV
104	AveragePt2 = 0.3; // GeV^2
105	};
106
107	ProjectileDiffCut = ProjectileDiffCut * GeV;
108	AveragePt2 = AveragePt2 * GeV*GeV;
109
110	// -------------------- Target parameters ----------------------------------------------
111	G4LorentzVector Ptarget=target->Get4Momentum();
112
113	G4double M0target = Ptarget.mag();
114
115	if(M0target < target->GetDefinition()->GetPDGMass())
116	{
117	PutOnMassShell=1;
118	M0target=target->GetDefinition()->GetPDGMass();
119	}
120
121	G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
122
123	G4double NuclearNucleonDiffCut = 1.1*GeV;
124
125	G4double ProjectileDiffCut2 = ProjectileDiffCut * ProjectileDiffCut;
126	G4double NuclearNucleonDiffCut2 = NuclearNucleonDiffCut * NuclearNucleonDiffCut;
127
128	// Transform momenta to cms and then rotate parallel to z axis;
129
130	G4LorentzVector Psum;
131	Psum=Pprojectile+Ptarget;
132
133	G4LorentzRotation toCms(-1*Psum.boostVector());
134
135	G4LorentzVector Ptmp=toCms*Pprojectile;
136
137	if ( Ptmp.pz() <= 0. )
138	{
139	// "String" moving backwards in CMS, abort collision !!
140	//G4cout << " abort Collision!! " << G4endl;
141	return false;
142	}
143
144	toCms.rotateZ(-1*Ptmp.phi());
145	toCms.rotateY(-1*Ptmp.theta());
146
147	G4LorentzRotation toLab(toCms.inverse());
148
149	Pprojectile.transform(toCms);
150	Ptarget.transform(toCms);
151
152	G4double Pt2;
153	G4double ProjMassT2, ProjMassT;
154	G4double TargMassT2, TargMassT;
155	G4double PZcms2, PZcms;
156	G4double PMinusNew, TPlusNew;
157
158	G4double S=Psum.mag2();
159	G4double SqrtS=std::sqrt(S);
160
161	if(SqrtS < 2200*MeV) {return false;} // The model cannot work for pp-interactions
162	// at Plab < 1.3 GeV/c. Uzhi
163
164	PZcms2=(SS+Mprojectile2Mprojectile2+Mtarget2*Mtarget2-
165	2SMprojectile2-2SMtarget2-2Mprojectile2Mtarget2)/4./S;
166	if(PZcms2 < 0)
167	{return false;} // It can be in an interaction with off-shell nuclear nucleon
168
169	PZcms = std::sqrt(PZcms2);
170
171	if(PutOnMassShell)
172	{
173	if(Pprojectile.z() > 0.)
174	{
175	Pprojectile.setPz( PZcms);
176	Ptarget.setPz( -PZcms);
177	}
178	else
179	{
180	Pprojectile.setPz(-PZcms);
181	Ptarget.setPz( PZcms);
182	};
183
184	Pprojectile.setE(std::sqrt(Mprojectile2+
185	Pprojectile.x()*Pprojectile.x()+
186	Pprojectile.y()*Pprojectile.y()+
187	PZcms2));
188	Ptarget.setE(std::sqrt( Mtarget2 +
189	Ptarget.x()*Ptarget.x()+
190	Ptarget.y()*Ptarget.y()+
191	PZcms2));
192	}
193
194	G4double maxPtSquare = PZcms2;
195
196	//G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
197	//G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
198	// G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
199
200	// G4cout << " Projectile Xplus / Xminus : " <<
201	// Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
202	// G4cout << " Target Xplus / Xminus : " <<
203	// Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
204
205	G4LorentzVector Qmomentum;
206	G4double Qminus, Qplus;
207
208	// /* Vova
209	G4int whilecount=0;
210	do {
211	// Generate pt
212
213	if (whilecount++ >= 500 && (whilecount%100)==0)
214	// G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
215	// << ", loop count/ maxPtSquare : "
216	// << whilecount << " / " << maxPtSquare << G4endl;
217	if (whilecount > 1000 )
218	{
219	Qmomentum=G4LorentzVector(0.,0.,0.,0.);
220	return false; // Ignore this interaction
221	}
222
223	Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
224
225	//G4cout << "generated Pt " << Qmomentum << G4endl;
226	//G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
227	//G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
228
229	// Momentum transfer
230	/* // Uzhi
231	G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
232	G4double Xmax=1.;
233	G4double Xplus =ChooseX(Xmin,Xmax);
234	G4double Xminus=ChooseX(Xmin,Xmax);
235
236	// G4cout << " X-plus " << Xplus << G4endl;
237	// G4cout << " X-minus " << Xminus << G4endl;
238
239	G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
240	G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
241	G4double Qminus= pt2 / Xplus /Pprojectile.plus();
242	/ // Uzhi
243
244	Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
245	ProjMassT2=Mprojectile2+Pt2;
246	ProjMassT =std::sqrt(ProjMassT2);
247
248	TargMassT2=Mtarget2+Pt2;
249	TargMassT =std::sqrt(TargMassT2);
250
251	PZcms2=(SS+ProjMassT2ProjMassT2+
252	TargMassT2*TargMassT2-
253	2.SProjMassT2-2.STargMassT2-
254	2.ProjMassT2TargMassT2)/4./S;
255	if(PZcms2 < 0 ) {PZcms2=0;};
256	PZcms =std::sqrt(PZcms2);
257
258	G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
259	G4double PMinusMax=SqrtS-TargMassT;
260
261	PMinusNew=ChooseP(PMinusMin,PMinusMax);
262	Qminus=PMinusNew-Pprojectile.minus();
263
264	G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
265	G4double TPlusMax=SqrtS-ProjMassT;
266
267	TPlusNew=ChooseP(TPlusMin, TPlusMax);
268	Qplus=-(TPlusNew-Ptarget.plus());
269
270	Qmomentum.setPz( (Qplus-Qminus)/2 );
271	Qmomentum.setE( (Qplus+Qminus)/2 );
272
273	//G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
274	// G4cout << "pt2" << pt2 << G4endl;
275	// G4cout << "Qmomentum " << Qmomentum << G4endl;
276	// G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
277	// " / " << (Ptarget-Qmomentum).mag() << G4endl;
278	/* // Uzhi
279	} while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 \|\|
280	(Ptarget-Qmomentum).mag2() <= Mtarget2 );
281	/ // Uzhi
282
283
284	} while (( (Pprojectile+Qmomentum).mag2() < Mprojectile2 \|\| // Uzhi No without excitation
285	(Ptarget -Qmomentum).mag2() < Mtarget2 ) \|\| // Uzhi
286	( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // Uzhi No double Diffraction
287	(Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );// Uzhi
288
289	if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi Projectile diffraction
290	{
291	G4double TMinusNew=SqrtS-PMinusNew;
292	Qminus=Ptarget.minus()-TMinusNew;
293	TPlusNew=TargMassT2/TMinusNew;
294	Qplus=Ptarget.plus()-TPlusNew;
295
296	Qmomentum.setPz( (Qplus-Qminus)/2 );
297	Qmomentum.setE( (Qplus+Qminus)/2 );
298	}
299	else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi Target diffraction
300	{
301	G4double PPlusNew=SqrtS-TPlusNew;
302	Qplus=PPlusNew-Pprojectile.plus();
303	PMinusNew=ProjMassT2/PPlusNew;
304	Qminus=PMinusNew-Pprojectile.minus();
305
306	Qmomentum.setPz( (Qplus-Qminus)/2 );
307	Qmomentum.setE( (Qplus+Qminus)/2 );
308	};
309
310	Pprojectile += Qmomentum;
311	Ptarget -= Qmomentum;
312
313	// Vova
314
315	/*
316	Pprojectile.setPz(0.);
317	Pprojectile.setE(SqrtS-M0target);
318
319	Ptarget.setPz(0.);
320	Ptarget.setE(M0target);
321	*/
322
323	//G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
324	//G4cout << "Ptarget with Q : " << Ptarget << G4endl;
325
326	// G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
327	// G4cout << "Target back: " << toLab * Ptarget << G4endl;
328
329	// Transform back and update SplitableHadron Participant.
330	Pprojectile.transform(toLab);
331	Ptarget.transform(toLab);
332
333	//G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
334	//G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
335
336	//G4cout << "Target mass " << Ptarget.mag() << G4endl;
337
338	target->Set4Momentum(Ptarget);
339
340	//G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
341
342	projectile->Set4Momentum(Pprojectile);
343
344	return true;
345	}
346
347
348	G4ExcitedString * G4QGSDiffractiveExcitation::
349	String(G4VSplitableHadron * hadron, G4bool isProjectile) const
350	{
351	hadron->SplitUp();
352	G4Parton *start= hadron->GetNextParton();
353	if ( start==NULL)
354	{ G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
355	return NULL;
356	}
357	G4Parton *end = hadron->GetNextParton();
358	if ( end==NULL)
359	{ G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
360	return NULL;
361	}
362
363	G4ExcitedString * string;
364	if ( isProjectile )
365	{
366	string= new G4ExcitedString(end,start, +1);
367	} else {
368	string= new G4ExcitedString(start,end, -1);
369	}
370
371	string->SetPosition(hadron->GetPosition());
372
373	// momenta of string ends
374	G4double ptSquared= hadron->Get4Momentum().perp2();
375	G4double transverseMassSquared= hadron->Get4Momentum().plus()
376	* hadron->Get4Momentum().minus();
377
378
379	G4double maxAvailMomentumSquared=
380	sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
381
382	G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ??????????????????
383	G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
384
385	G4LorentzVector Pstart(G4LorentzVector(pt,0.));
386	G4LorentzVector Pend;
387	Pend.setPx(hadron->Get4Momentum().px() - pt.x());
388	Pend.setPy(hadron->Get4Momentum().py() - pt.y());
389
390	G4double tm1=hadron->Get4Momentum().minus() +
391	( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
392
393	G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
394	4. * Pend.perp2() * hadron->Get4Momentum().minus()
395	/ hadron->Get4Momentum().plus() ));
396
397	G4int Sign= isProjectile ? -1 : 1;
398
399	G4double endMinus = 0.5 * (tm1 + Sign*tm2);
400	G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
401
402	G4double startPlus= Pstart.perp2() / startMinus;
403	G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
404
405	Pstart.setPz(0.5*(startPlus - startMinus));
406	Pstart.setE(0.5*(startPlus + startMinus));
407
408	Pend.setPz(0.5*(endPlus - endMinus));
409	Pend.setE(0.5*(endPlus + endMinus));
410
411	start->Set4Momentum(Pstart);
412	end->Set4Momentum(Pend);
413
414	#ifdef G4_FTFDEBUG
415	G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
416	G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
417	G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
418	G4cout << " sum of ends " << Pstart+Pend << G4endl;
419	G4cout << " Original " << hadron->Get4Momentum() << G4endl;
420	#endif
421
422	return string;
423	}
424
425
426	// --------- private methods ----------------------
427
428	G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const // Uzhi
429	{
430	// choose an x between Xmin and Xmax with P(x) ~ 1/x
431
432	// to be improved...
433
434	G4double range=Pmax-Pmin; // Uzhi
435
436	if ( Pmin <= 0. \|\| range <=0. )
437	{
438	G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
439	throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
440	}
441
442	G4double P;
443	/* // Uzhi
444	do {
445	x=Xmin + G4UniformRand() * range;
446	} while ( Xmin/x < G4UniformRand() );
447	*/ // Uzhi
448
449	P=Pmin * std::pow(Pmax/Pmin,G4UniformRand()); // Uzhi
450
451	//debug-hpw cout << "DiffractiveX "<<x<<G4endl;
452	return P;
453	}
454
455	G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const // Uzhi
456	{ // @@ this method is used in FTFModel as well. Should go somewhere common!
457
458	G4double Pt2;
459	/* // Uzhi
460	do {
461	pt2=widthSquare * std::log( G4UniformRand() );
462	} while ( pt2 > maxPtSquare);
463	*/ // Uzhi
464
465	Pt2 = -AveragePt2 * std::log(1. + G4UniformRand() * (std::exp(-maxPtSquare/AveragePt2)-1.));// Uzhi
466
467	G4double Pt=std::sqrt(Pt2);
468
469	G4double phi=G4UniformRand() * twopi;
470
471	return G4ThreeVector (Ptstd::cos(phi), Ptstd::sin(phi), 0.);
472	}
473
474	G4QGSDiffractiveExcitation::G4QGSDiffractiveExcitation(const G4QGSDiffractiveExcitation &)
475	{
476	throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation copy contructor not meant to be called");
477	}
478
479
480	G4QGSDiffractiveExcitation::~G4QGSDiffractiveExcitation()
481	{
482	}
483
484
485	const G4QGSDiffractiveExcitation & G4QGSDiffractiveExcitation::operator=(const G4QGSDiffractiveExcitation &)
486	{
487	throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation = operator meant to be called");
488	return *this;
489	}
490
491
492	int G4QGSDiffractiveExcitation::operator==(const G4QGSDiffractiveExcitation &) const
493	{
494	throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation == operator meant to be called");
495	return false;
496	}
497
498	int G4QGSDiffractiveExcitation::operator!=(const G4QGSDiffractiveExcitation &) const
499	{
500	throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation != operator meant to be called");
501	return true;
502	}

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

Context Navigation

source: trunk/source/processes/hadronic/models/parton_string/qgsm/src/G4QGSDiffractiveExcitation.cc @ 1198

Download in other formats: