ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SingleDiffractiveExcitation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SingleDiffractiveExcitation.cc
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 // ------------------------------------------------------------
28 // GEANT 4 class implemetation file
29 //
30 // ---------------- G4SingleDiffractiveExcitation --------------
31 // by Gunter Folger, October 1998.
32 // diffractive Excitation used by strings models
33 // Take a projectile and a target
34 // excite the projectile and target
35 // ------------------------------------------------------------
36 
37 #include "G4SingleDiffractiveExcitation.hh"
38 #include "globals.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "Randomize.hh"
42 #include "G4LorentzRotation.hh"
43 #include "G4ThreeVector.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4VSplitableHadron.hh"
46 #include "G4ExcitedString.hh"
47 
48 #include "G4Log.hh"
49 #include "G4Pow.hh"
50 
51 //#define debugSingleDiffraction
52 
53 G4SingleDiffractiveExcitation::G4SingleDiffractiveExcitation(){}
54 
55 G4SingleDiffractiveExcitation::~G4SingleDiffractiveExcitation(){}
56 
57 G4bool G4SingleDiffractiveExcitation::
58 ExciteParticipants( G4VSplitableHadron *projectile, G4VSplitableHadron *target,
59  G4bool ProjectileDiffraction ) const
60 {
61  #ifdef debugSingleDiffraction
62  G4cout<<G4endl<<"G4SingleDiffractiveExcitation::ExciteParticipants"<<G4endl;
63  #endif
64 
65  G4LorentzVector Pprojectile=projectile->Get4Momentum();
66  G4double Mprojectile = projectile->GetDefinition()->GetPDGMass();
67  G4double Mprojectile2=sqr(projectile->GetDefinition()->GetPDGMass());
68 
69  G4LorentzVector Ptarget=target->Get4Momentum();
70  G4double Mtarget = target->GetDefinition()->GetPDGMass();
71  G4double Mtarget2=sqr(target->GetDefinition()->GetPDGMass());
72 
73  #ifdef debugSingleDiffraction
74  G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl;
75  G4cout<<"Pr Tr 4-Mom "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl
76  <<" "<<Ptarget <<" "<<Ptarget.mag() <<G4endl;
77  #endif
78 
79  G4LorentzVector Psum=Pprojectile+Ptarget;
80  G4double SqrtS=Psum.mag();
81  G4double S =Psum.mag2();
82 
83  #ifdef debugSingleDiffraction
84  G4cout<<"SqrtS-Mprojectile-Mtarget "<<SqrtS<<" "<<Mprojectile<<" "<<Mtarget
85  <<" "<<SqrtS-Mprojectile-Mtarget<<G4endl;
86  #endif
87  if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) {
88  #ifdef debugSingleDiffraction
89  G4cerr<<"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<" "
90  <<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
91  G4cerr<<"Target: "<<target->GetDefinition()->GetPDGEncoding()<<" "
92  <<Ptarget<<" "<<Ptarget.mag()<<G4endl;
93  G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "<<Pprojectile.mag()+Ptarget.mag()<<G4endl;
94  #endif
95  return true;
96  }
97 
98  G4LorentzRotation toCms(-1*Psum.boostVector());
99 
100  G4LorentzVector Ptmp=toCms*Pprojectile;
101 
102  if ( Ptmp.pz() <= 0. )
103  {
104  // "String" moving backwards in CMS, abort collision !!
105  // G4cout << " abort Collision!! " << G4endl;
106  return false;
107  }
108 
109  toCms.rotateZ(-1*Ptmp.phi());
110  toCms.rotateY(-1*Ptmp.theta());
111 
112  G4LorentzRotation toLab(toCms.inverse());
113 
114  Pprojectile.transform(toCms);
115  Ptarget.transform(toCms);
116  #ifdef debugSingleDiffraction
117  G4cout << "Pprojectile in CMS " << Pprojectile << G4endl;
118  G4cout << "Ptarget in CMS " << Ptarget << G4endl;
119  #endif
120  G4double maxPtSquare=sqr(Ptarget.pz());
121 
122  G4double ProjectileMinDiffrMass(0.), TargetMinDiffrMass(0.);
123  G4double AveragePt2(0.);
124  G4int absPDGcode=std::abs(projectile->GetDefinition()->GetPDGEncoding());
125 
126  if ( ProjectileDiffraction ) {
127  if ( absPDGcode > 1000 ) //------Projectile is baryon --------
128  {
129  ProjectileMinDiffrMass = 1.16; // GeV
130  AveragePt2 = 0.3; // GeV^2
131  }
132  else if( absPDGcode == 211 || absPDGcode == 111) //------Projectile is Pion -----------
133  {
134  ProjectileMinDiffrMass = 1.0; // GeV
135  AveragePt2 = 0.3; // GeV^2
136  }
137  else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon
138  {
139  ProjectileMinDiffrMass = 1.1; // GeV
140  AveragePt2 = 0.3; // GeV^2
141  }
142  else if( absPDGcode == 22) //------Projectile is Gamma -----------
143  {
144  ProjectileMinDiffrMass = 0.25; // GeV
145  AveragePt2 = 0.36; // GeV^2
146  }
147  else //------Projectile is undefined, Nucleon assumed
148  {
149  ProjectileMinDiffrMass = 1.1; // GeV
150  AveragePt2 = 0.3; // GeV^2
151  };
152 
153  ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV;
154  Mprojectile2=sqr(ProjectileMinDiffrMass);
155  }
156  else
157  {
158  TargetMinDiffrMass = 1.16*GeV; // For target nucleon
159  Mtarget2 = sqr( TargetMinDiffrMass) ;
160  AveragePt2 = 0.3; // GeV^2
161  } // end of if ( ProjectileDiffraction )
162 
163  AveragePt2 = AveragePt2 * GeV*GeV;
164 
165  G4double Pt2, PZcms, PZcms2;
166  G4double ProjMassT2, ProjMassT;
167  G4double TargMassT2, TargMassT;
168  G4double PMinusMin, PMinusMax;
169  G4double TPlusMin, TPlusMax;
170  G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
171 
172  G4LorentzVector Qmomentum;
173  G4double Qminus, Qplus;
174 
175  G4int whilecount=0;
176  do {
177  whilecount++;
178 
179  if (whilecount > 1000 )
180  {
181  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
182  return false; // Ignore this interaction
183  }
184 
185  // Generate pt
186  Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
187 
188  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
189 
190  ProjMassT2 = Mprojectile2 + Pt2;
191  ProjMassT = std::sqrt( ProjMassT2 );
192  TargMassT2 = Mtarget2 + Pt2;
193  TargMassT = std::sqrt( TargMassT2 );
194 
195  #ifdef debugSingleDiffraction
196  G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl;
197  #endif
198  if ( SqrtS < ProjMassT + TargMassT ) continue;
199 
200  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
201  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
202 
203  if ( PZcms2 < 0 ) continue;
204 
205  PZcms = std::sqrt( PZcms2 );
206 
207  if ( ProjectileDiffraction )
208  { // The projectile will fragment, the target will saved.
209  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
210  PMinusMax = SqrtS - TargMassT;
211 
212  PMinusNew = ChooseX( PMinusMin, PMinusMax );
213  TMinusNew = SqrtS - PMinusNew;
214 
215  Qminus = Ptarget.minus() - TMinusNew;
216  TPlusNew = TargMassT2 / TMinusNew;
217  Qplus = Ptarget.plus() - TPlusNew;
218 
219  } else { // The target will fragment, the projectile will saved.
220  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
221  TPlusMax = SqrtS - ProjMassT;
222 
223  TPlusNew = ChooseX( TPlusMin, TPlusMax );
224  PPlusNew = SqrtS - TPlusNew;
225 
226  Qplus = PPlusNew - Pprojectile.plus();
227  PMinusNew = ProjMassT2 / PPlusNew;
228  Qminus = PMinusNew - Pprojectile.minus();
229  }
230 
231  Qmomentum.setPz( (Qplus - Qminus)/2 );
232  Qmomentum.setE( (Qplus + Qminus)/2 );
233 
234  #ifdef debugSingleDiffraction
235  G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl;
236  G4cout<<!ProjectileDiffraction<<" "<<( Ptarget - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl;
237  #endif
238 
239  } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) ||
240  (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) );
241  // Repeat the sampling because there was not any excitation
242 
243  Pprojectile += Qmomentum;
244 
245  Ptarget -= Qmomentum;
246 
247  // Transform back and update SplitableHadron Participant.
248  Pprojectile.transform(toLab);
249  Ptarget.transform(toLab);
250 
251  #ifdef debugSingleDiffraction
252  G4cout << "Pprojectile in Lab. " << Pprojectile << G4endl;
253  G4cout << "Ptarget in Lab. " << Ptarget << G4endl;
254  G4cout << "G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() << G4endl;
255  G4cout << "G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() << G4endl;
256  #endif
257 
258  target->Set4Momentum(Ptarget);
259  projectile->Set4Momentum(Pprojectile);
260 
261  return true;
262 }
263 
264 // --------- private methods ----------------------
265 
266 G4double G4SingleDiffractiveExcitation::ChooseX(G4double Xmin, G4double Xmax) const
267 {
268  // choose an x between Xmin and Xmax with P(x) ~ 1/x
269  G4double range=Xmax-Xmin;
270 
271  if ( Xmin <= 0. || range <=0. )
272  {
273  G4cout << " Xmin, range : " << Xmin << " , " << range << G4endl;
274  throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation::ChooseX : Invalid arguments ");
275  }
276 
277  G4double x = Xmin*G4Pow::GetInstance()->powA(Xmax/Xmin, G4UniformRand() );
278  // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin) - G4UniformRand() * (1.0/std::sqrt(Xmin) - 1.0/std::sqrt(Xmax)));
279  return x;
280 }
281 
282 
283 G4ThreeVector G4SingleDiffractiveExcitation::GaussianPt(G4double widthSquare, G4double maxPtSquare) const
284 { // @@ this method is used in FTFModel as well. Should go somewhere common!
285 
286  G4double pt2;
287 
288  const G4int maxNumberOfLoops = 1000;
289  G4int loopCounter = 0;
290  do {
291  pt2=-widthSquare * G4Log( G4UniformRand() );
292  } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
293  if ( loopCounter >= maxNumberOfLoops ) {
294  pt2 = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration.
295  }
296 
297  pt2=std::sqrt(pt2);
298 
299  G4double phi=G4UniformRand() * twopi;
300 
301  return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);
302 }
303