Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QGSDiffractiveExcitation.cc
Go to the documentation of this file.
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 100828 2016-11-02 15:25:59Z gcosmo $
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 "G4PhysicalConstants.hh"
48 #include "G4SystemOfUnits.hh"
49 #include "Randomize.hh"
50 
52 #include "G4LorentzRotation.hh"
53 #include "G4ThreeVector.hh"
54 #include "G4ParticleDefinition.hh"
55 #include "G4VSplitableHadron.hh"
56 #include "G4ExcitedString.hh"
57 //#include "G4ios.hh"
58 
59 #include "G4Exp.hh"
60 #include "G4Log.hh"
61 #include "G4Pow.hh"
62 
63 
65 {
66 }
67 
69 {
70 }
71 
72 
75 {
76  G4LorentzVector Pprojectile=projectile->Get4Momentum();
77 
78  // -------------------- Projectile parameters -----------------------------------
79  G4bool PutOnMassShell=0;
80 
81  //G4double M0projectile=projectile->GetDefinition()->GetPDGMass(); // With de-excitation
82  G4double M0projectile = Pprojectile.mag(); // Without de-excitation
83 
84  if(M0projectile < projectile->GetDefinition()->GetPDGMass())
85  {
86  PutOnMassShell=1;
87  M0projectile=projectile->GetDefinition()->GetPDGMass();
88  }
89 
90  G4double Mprojectile2 = M0projectile * M0projectile;
91 
92  G4int PDGcode=projectile->GetDefinition()->GetPDGEncoding();
93  G4int absPDGcode=std::abs(PDGcode);
94  G4double ProjectileDiffCut;
95  G4double AveragePt2;
96 
97  if( absPDGcode > 1000 ) //------Projectile is baryon --------
98  {
99  ProjectileDiffCut = 1.1; // GeV
100  AveragePt2 = 0.3; // GeV^2
101  }
102  else if( absPDGcode == 211 || PDGcode == 111) //------Projectile is Pion -----------
103  {
104  ProjectileDiffCut = 1.0; // GeV
105  AveragePt2 = 0.3; // GeV^2
106  }
107  else if( absPDGcode == 321 || PDGcode == -311) //------Projectile is Kaon -----------
108  {
109  ProjectileDiffCut = 1.1; // GeV
110  AveragePt2 = 0.3; // GeV^2
111  }
112  else //------Projectile is undefined, Nucleon assumed
113  {
114  ProjectileDiffCut = 1.1; // GeV
115  AveragePt2 = 0.3; // GeV^2
116  };
117 
118  ProjectileDiffCut = ProjectileDiffCut * GeV;
119  AveragePt2 = AveragePt2 * GeV*GeV;
120 
121  // -------------------- Target parameters ----------------------------------------------
122  G4LorentzVector Ptarget=target->Get4Momentum();
123 
124  G4double M0target = Ptarget.mag();
125 
126  if(M0target < target->GetDefinition()->GetPDGMass())
127  {
128  PutOnMassShell=1;
129  M0target=target->GetDefinition()->GetPDGMass();
130  }
131 
132  G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
133 
134  G4double NuclearNucleonDiffCut = 1.1*GeV;
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(SqrtS < 2200*MeV) {return false;} // The model cannot work for pp-interactions
173  // at Plab < 1.3 GeV/c. Uzhi
174 
175  PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
176  2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
177  if(PZcms2 < 0)
178  {return false;} // It can be in an interaction with off-shell nuclear nucleon
179 
180  PZcms = std::sqrt(PZcms2);
181 
182  if(PutOnMassShell)
183  {
184  if(Pprojectile.z() > 0.)
185  {
186  Pprojectile.setPz( PZcms);
187  Ptarget.setPz( -PZcms);
188  } else
189  {
190  Pprojectile.setPz(-PZcms);
191  Ptarget.setPz( PZcms);
192  };
193 
194  Pprojectile.setE(std::sqrt(Mprojectile2+
195  Pprojectile.x()*Pprojectile.x()+
196  Pprojectile.y()*Pprojectile.y()+
197  PZcms2));
198  Ptarget.setE(std::sqrt( Mtarget2 +
199  Ptarget.x()*Ptarget.x()+
200  Ptarget.y()*Ptarget.y()+
201  PZcms2));
202  }
203 
204  G4double maxPtSquare = PZcms2;
205 
206  //G4cout << "Pprojectile aft boost : " << Pprojectile << G4endl;
207  //G4cout << "Ptarget aft boost : " << Ptarget << G4endl;
208  // G4cout << "cms aft boost : " << (Pprojectile+ Ptarget) << G4endl;
209  // G4cout << " Projectile Xplus / Xminus : " <<
210  // Pprojectile.plus() << " / " << Pprojectile.minus() << G4endl;
211  // G4cout << " Target Xplus / Xminus : " <<
212  // Ptarget.plus() << " / " << Ptarget.minus() << G4endl;
213 
214  G4LorentzVector Qmomentum;
215  G4double Qminus, Qplus;
216 
217  G4int whilecount=0;
218  do {
219  // Generate pt
220 
221  if (whilecount++ >= 500 && (whilecount%100)==0)
222  //G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
223  // << ", loop count/ maxPtSquare : "
224  // << whilecount << " / " << maxPtSquare << G4endl;
225 
226  if (whilecount > 1000 )
227  {
228  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
229  return false; // Ignore this interaction
230  }
231 
232  Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
233 
234  //G4cout << "generated Pt " << Qmomentum << G4endl;
235  //G4cout << "Pprojectile with pt : " << Pprojectile+Qmomentum << G4endl;
236  //G4cout << "Ptarget with pt : " << Ptarget-Qmomentum << G4endl;
237 
238  // Momentum transfer
239  /*
240  G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
241  G4double Xmax=1.;
242  G4double Xplus =ChooseX(Xmin,Xmax);
243  G4double Xminus=ChooseX(Xmin,Xmax);
244 
245  //G4cout << " X-plus " << Xplus << G4endl;
246  //G4cout << " X-minus " << Xminus << G4endl;
247 
248  G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
249  G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
250  G4double Qminus= pt2 / Xplus /Pprojectile.plus();
251  */
252 
253  Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
254  ProjMassT2=Mprojectile2+Pt2;
255  ProjMassT =std::sqrt(ProjMassT2);
256 
257  TargMassT2=Mtarget2+Pt2;
258  TargMassT =std::sqrt(TargMassT2);
259 
260  PZcms2=(S*S+ProjMassT2*ProjMassT2+
261  TargMassT2*TargMassT2-
262  2.*S*ProjMassT2-2.*S*TargMassT2-
263  2.*ProjMassT2*TargMassT2)/4./S;
264  if(PZcms2 < 0 ) {PZcms2=0;};
265  PZcms =std::sqrt(PZcms2);
266 
267  G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
268  G4double PMinusMax=SqrtS-TargMassT;
269 
270  PMinusNew=ChooseP(PMinusMin,PMinusMax);
271  Qminus=PMinusNew-Pprojectile.minus();
272 
273  G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
274  G4double TPlusMax=SqrtS-ProjMassT;
275 
276  TPlusNew=ChooseP(TPlusMin, TPlusMax);
277  Qplus=-(TPlusNew-Ptarget.plus());
278 
279  Qmomentum.setPz( (Qplus-Qminus)/2 );
280  Qmomentum.setE( (Qplus+Qminus)/2 );
281 
282  //G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
283  //G4cout << "pt2" << pt2 << G4endl;
284  //G4cout << "Qmomentum " << Qmomentum << G4endl;
285  //G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
286  // " / " << (Ptarget-Qmomentum).mag() << G4endl;
287  /*
288  } while ( (Pprojectile+Qmomentum).mag2() <= Mprojectile2 ||
289  (Ptarget-Qmomentum).mag2() <= Mtarget2 );
290  */
291  } while ( /* Loop checking, 26.10.2015, A.Ribon */
292  ( (Pprojectile+Qmomentum).mag2() < Mprojectile2 || // No without excitation
293  (Ptarget -Qmomentum).mag2() < Mtarget2 ) ||
294  ( (Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2 && // No double Diffraction
295  (Ptarget -Qmomentum).mag2() < NuclearNucleonDiffCut2) );
296 
297  if((Ptarget-Qmomentum).mag2() < NuclearNucleonDiffCut2) // Uzhi Projectile diffraction
298  {
299  G4double TMinusNew=SqrtS-PMinusNew;
300  Qminus=Ptarget.minus()-TMinusNew;
301  TPlusNew=TargMassT2/TMinusNew;
302  Qplus=Ptarget.plus()-TPlusNew;
303 
304  Qmomentum.setPz( (Qplus-Qminus)/2 );
305  Qmomentum.setE( (Qplus+Qminus)/2 );
306  }
307  else if((Pprojectile+Qmomentum).mag2() < ProjectileDiffCut2) // Uzhi Target diffraction
308  {
309  G4double PPlusNew=SqrtS-TPlusNew;
310  Qplus=PPlusNew-Pprojectile.plus();
311  PMinusNew=ProjMassT2/PPlusNew;
312  Qminus=PMinusNew-Pprojectile.minus();
313 
314  Qmomentum.setPz( (Qplus-Qminus)/2 );
315  Qmomentum.setE( (Qplus+Qminus)/2 );
316  };
317 
318  Pprojectile += Qmomentum;
319  Ptarget -= Qmomentum;
320 
321  // Vova
322 
323  /*
324  Pprojectile.setPz(0.);
325  Pprojectile.setE(SqrtS-M0target);
326  Ptarget.setPz(0.);
327  Ptarget.setE(M0target);
328  */
329 
330  //G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
331  //G4cout << "Ptarget with Q : " << Ptarget << G4endl;
332  //G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
333  //G4cout << "Target back: " << toLab * Ptarget << G4endl;
334 
335  // Transform back and update SplitableHadron Participant.
336  Pprojectile.transform(toLab);
337  Ptarget.transform(toLab);
338 
339  //G4cout << "Pprojectile with Q M: " << Pprojectile<<" "<< Pprojectile.mag() << G4endl;
340  //G4cout << "Ptarget with Q M: " << Ptarget <<" "<< Ptarget.mag() << G4endl;
341  //G4cout << "Target mass " << Ptarget.mag() << G4endl;
342 
343  target->Set4Momentum(Ptarget);
344 
345  //G4cout << "Projectile mass " << Pprojectile.mag() << G4endl;
346 
347  projectile->Set4Momentum(Pprojectile);
348 
349  return true;
350 }
351 
352 
354 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
355 {
356  hadron->SplitUp();
357  G4Parton *start= hadron->GetNextParton();
358  if ( start==NULL) {
359  G4cout << " G4FTFModel::String() Error:No start parton found"<< G4endl;
360  return NULL;
361  }
362  G4Parton *end = hadron->GetNextParton();
363  if ( end==NULL) {
364  G4cout << " G4FTFModel::String() Error:No end parton found"<< G4endl;
365  return NULL;
366  }
367 
368  G4ExcitedString * string;
369  if ( isProjectile ) {
370  string= new G4ExcitedString(end,start, +1);
371  } else {
372  string= new G4ExcitedString(start,end, -1);
373  }
374 
375  string->SetPosition(hadron->GetPosition());
376 
377  // momenta of string ends
378  G4double ptSquared= hadron->Get4Momentum().perp2();
379  G4double transverseMassSquared= hadron->Get4Momentum().plus()
380  * hadron->Get4Momentum().minus();
381 
382  G4double maxAvailMomentumSquared=
383  sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
384 
385  G4double widthOfPtSquare = 0.25; // Uzhi <Pt^2>=0.25 ???
386  G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
387 
388  G4LorentzVector Pstart(G4LorentzVector(pt,0.));
389  G4LorentzVector Pend;
390  Pend.setPx(hadron->Get4Momentum().px() - pt.x());
391  Pend.setPy(hadron->Get4Momentum().py() - pt.y());
392 
393  G4double tm1=hadron->Get4Momentum().minus() +
394  ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
395 
396  G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
397  4. * Pend.perp2() * hadron->Get4Momentum().minus()
398  / hadron->Get4Momentum().plus() ));
399 
400  G4int Sign= isProjectile ? -1 : 1;
401 
402  G4double endMinus = 0.5 * (tm1 + Sign*tm2);
403  G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
404 
405  G4double startPlus= Pstart.perp2() / startMinus;
406  G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
407 
408  Pstart.setPz(0.5*(startPlus - startMinus));
409  Pstart.setE(0.5*(startPlus + startMinus));
410 
411  Pend.setPz(0.5*(endPlus - endMinus));
412  Pend.setE(0.5*(endPlus + endMinus));
413 
414  start->Set4Momentum(Pstart);
415  end->Set4Momentum(Pend);
416 
417  #ifdef G4_FTFDEBUG
418  G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
419  G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
420  G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
421  G4cout << " sum of ends " << Pstart+Pend << G4endl;
422  G4cout << " Original " << hadron->Get4Momentum() << G4endl;
423  #endif
424 
425  return string;
426 }
427 
428 
429 // --------- private methods ----------------------
430 
431 G4double G4QGSDiffractiveExcitation::ChooseP(G4double Pmin, G4double Pmax) const
432 {
433  // choose an x between Xmin and Xmax with P(x) ~ 1/x
434  // to be improved...
435 
436  G4double range=Pmax-Pmin;
437 
438  if ( Pmin <= 0. || range <=0. )
439  {
440  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
441  throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
442  }
443 
444  G4double P;
445  /*
446  do {
447  x=Xmin + G4UniformRand() * range;
448  } while ( Xmin/x < G4UniformRand() );
449  */
450 
451  P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmin,G4UniformRand());
452 
453  //debug-hpw cout << "DiffractiveX "<<x<<G4endl;
454  return P;
455 }
456 
457 G4ThreeVector G4QGSDiffractiveExcitation::GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
458 { // @@ this method is used in FTFModel as well. Should go somewhere common!
459 
460  G4double Pt2;
461  /*
462  do {
463  pt2=widthSquare * G4Log( G4UniformRand() );
464  } while ( pt2 > maxPtSquare);
465  */
466 
467  Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.));
468 
469  G4double Pt=std::sqrt(Pt2);
470 
471  G4double phi=G4UniformRand() * twopi;
472 
473  return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
474 }
475 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
Hep3Vector boostVector() const
G4int GetPDGcode() const
Definition: G4Parton.hh:127
virtual G4ExcitedString * String(G4VSplitableHadron *aHadron, G4bool isProjectile) const
const XML_Char * target
Definition: expat.h:268
double S(double temp)
CLHEP::Hep3Vector G4ThreeVector
double x() const
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148
HepLorentzVector & rotateZ(double)
int G4int
Definition: G4Types.hh:78
static double P[]
virtual void SplitUp()=0
static constexpr double twopi
Definition: G4SIunits.hh:76
const G4ParticleDefinition * GetDefinition() const
Hep3Vector vect() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition const G4Material *G4double range
double mag() const
double py() const
bool G4bool
Definition: G4Types.hh:79
virtual G4Parton * GetNextParton()=0
double px() const
HepLorentzRotation & transform(const HepBoost &b)
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
double minus() const
const G4LorentzVector & Get4Momentum() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4double GetPDGMass() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
double mag2() const
double y() const
HepLorentzVector & rotateY(double)
static constexpr double GeV
Definition: G4SIunits.hh:217
double perp2() const
const G4ThreeVector & GetPosition() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
HepLorentzVector & transform(const HepRotation &)
void Set4Momentum(const G4LorentzVector &a4Momentum)
T sqr(const T &x)
Definition: templates.hh:145
double plus() const
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector