ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4QGSDiffractiveExcitation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4QGSDiffractiveExcitation.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 // ---------------- G4QGSDiffractiveExcitation --------------
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 // Essential changed by V. Uzhinsky in November - December 2006
36 // in order to put it in a correspondence with original FRITIOF
37 // model. Variant of FRITIOF with nucleon de-excitation is implemented.
38 // ---------------------------------------------------------------------
39 
40 // Modified:
41 // 25-05-07 : G.Folger
42 // move from management/G4DiffractiveExcitation to to qgsm/G4QGSDiffractiveExcitation
43 //
44 
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh"
49 
51 #include "G4LorentzRotation.hh"
52 #include "G4ThreeVector.hh"
53 #include "G4ParticleDefinition.hh"
54 #include "G4VSplitableHadron.hh"
55 #include "G4ExcitedString.hh"
56 //#include "G4ios.hh"
57 
58 #include "G4Exp.hh"
59 #include "G4Log.hh"
60 #include "G4Pow.hh"
61 
62 //============================================================================
63 
64 //#define debugDoubleDiffraction
65 
66 //============================================================================
67 
69 {
70 }
71 
73 {
74 }
75 
76 
79 {
80  #ifdef debugDoubleDiffraction
81  G4cout<<G4endl<<"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<G4endl;
82  G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetParticleName()<<" "<<target->GetDefinition()->GetParticleName()<<G4endl;
83  G4cout<<"Proj 4 Mom "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
84  G4cout<<"Targ 4 Mom "<<target->Get4Momentum() <<" "<<target->Get4Momentum().mag() <<G4endl;
85  #endif
86 
87  G4LorentzVector Pprojectile=projectile->Get4Momentum();
88 
89  // -------------------- Projectile parameters -----------------------------------
90  G4bool PutOnMassShell=0;
91 
92  G4double M0projectile = Pprojectile.mag(); // Without de-excitation
93 
94  if(M0projectile < projectile->GetDefinition()->GetPDGMass())
95  {
96  PutOnMassShell=1;
97  M0projectile=projectile->GetDefinition()->GetPDGMass();
98  }
99 
100  // -------------------- Target parameters ----------------------------------------------
101  G4LorentzVector Ptarget=target->Get4Momentum();
102 
103  G4double M0target = Ptarget.mag();
104 
105  if(M0target < target->GetDefinition()->GetPDGMass())
106  {
107  PutOnMassShell=1;
108  M0target=target->GetDefinition()->GetPDGMass();
109  }
110 
111  G4LorentzVector Psum=Pprojectile+Ptarget;
112  G4double S=Psum.mag2();
113  G4double SqrtS=std::sqrt(S);
114 
115  if (SqrtS < M0projectile + M0target) {return false;} // The model cannot work for pp-interactions
116  // at Plab < 1.3 GeV/c.
117 
118  G4double Mprojectile2 = M0projectile * M0projectile;
119  G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
120 
121  // Transform momenta to cms and then rotate parallel to z axis;
122 
123  G4LorentzRotation toCms(-1*Psum.boostVector());
124 
125  G4LorentzVector Ptmp=toCms*Pprojectile;
126 
127  if ( Ptmp.pz() <= 0. )
128  {
129  // "String" moving backwards in CMS, abort collision !!
130  //G4cout << " abort Collision!! " << G4endl;
131  return false;
132  }
133 
134  toCms.rotateZ(-1*Ptmp.phi());
135  toCms.rotateY(-1*Ptmp.theta());
136 
137  G4LorentzRotation toLab(toCms.inverse());
138 
139  Pprojectile.transform(toCms);
140  Ptarget.transform(toCms);
141 
142  G4double PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
143  2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
144 
145  if (PZcms2 < 0) {return false;} // It can be in an interaction with off-shell nuclear nucleon
146 
147  G4double PZcms = std::sqrt(PZcms2);
148 
149  if (PutOnMassShell)
150  {
151  if (Pprojectile.z() > 0.)
152  {
153  Pprojectile.setPz( PZcms);
154  Ptarget.setPz( -PZcms);
155  }
156  else
157  {
158  Pprojectile.setPz(-PZcms);
159  Ptarget.setPz( PZcms);
160  };
161 
162  Pprojectile.setE(std::sqrt(Mprojectile2+sqr(Pprojectile.x())+sqr(Pprojectile.y())+PZcms2));
163  Ptarget.setE( std::sqrt( Mtarget2+sqr( Ptarget.x())+sqr( Ptarget.y())+PZcms2));
164  }
165 
166  G4double maxPtSquare = PZcms2;
167 
168  #ifdef debugDoubleDiffraction
169  G4cout << "Pprojectile after boost to CMS: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
170  G4cout << "Ptarget after boost to CMS: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
171  #endif
172 
173  G4int PrPDGcode=projectile->GetDefinition()->GetPDGEncoding();
174  G4int absPrPDGcode=std::abs(PrPDGcode);
175  G4double MinPrDiffMass(0.);
176  G4double AveragePt2(0.);
177 
178  if (M0projectile <= projectile->GetDefinition()->GetPDGMass())
179  { // Normal projectile
180  if( absPrPDGcode > 1000 ) //------Projectile is baryon --------
181  {
182  MinPrDiffMass = 1.16; // GeV
183  AveragePt2 = 0.3; // GeV^2
184  }
185  else if( absPrPDGcode == 211 || PrPDGcode == 111) //------Projectile is Pion -----------
186  {
187  MinPrDiffMass = 1.0; // GeV
188  AveragePt2 = 0.3; // GeV^2
189  }
190  else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310) //-Projectile is Kaon-
191  {
192  MinPrDiffMass = 1.1; // GeV
193  AveragePt2 = 0.3; // GeV^2
194  }
195  else //------Projectile is undefined, Nucleon assumed
196  {
197  MinPrDiffMass = 1.16; // GeV
198  AveragePt2 = 0.3; // GeV^2
199  }
200  }
201  else
202  { // Excited projectile
203  MinPrDiffMass = M0projectile + 220.0*MeV;
204  AveragePt2 = 0.3;
205  }
206 
207  MinPrDiffMass = MinPrDiffMass * GeV;
208  AveragePt2 = AveragePt2 * GeV*GeV;
209  //---------------------------------------------
210  G4double MinTrDiffMass = 1.16*GeV;
211 
212  if (SqrtS < MinPrDiffMass + MinTrDiffMass) {return false;} // The model cannot work at low energy
213 
214  G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
215  G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
216 
217  G4double Pt2;
218  G4double ProjMassT2, ProjMassT;
219  G4double TargMassT2, TargMassT;
220  G4double PMinusNew, TPlusNew;
221 
222  G4LorentzVector Qmomentum;
223  G4double Qminus, Qplus;
224 
225  G4int whilecount=0;
226  do {
227  if (whilecount++ >= 500 && (whilecount%100)==0)
228  if (whilecount > 1000 ) {
229  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
230  return false; // Ignore this interaction
231  }
232 
233  // Generate pt
234  Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
235 
236  Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
237  ProjMassT2=MinPrDiffMass2+Pt2;
238  ProjMassT =std::sqrt(ProjMassT2);
239 
240  TargMassT2=MinTrDiffMass2+Pt2;
241  TargMassT =std::sqrt(TargMassT2);
242 
243  if (SqrtS < ProjMassT + TargMassT) continue;
244 
245  PZcms2=(S*S+ProjMassT2*ProjMassT2+TargMassT2*TargMassT2-
246  2.*S*ProjMassT2-2.*S*TargMassT2-2.*ProjMassT2*TargMassT2)/4./S;
247  if (PZcms2 < 0 ) {PZcms2=0;};
248  PZcms =std::sqrt(PZcms2);
249 
250  G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
251  G4double PMinusMax=SqrtS-TargMassT;
252 
253  PMinusNew=ChooseP(PMinusMin,PMinusMax);
254  Qminus=PMinusNew-Pprojectile.minus();
255 
256  G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
257  G4double TPlusMax=SqrtS-ProjMassT;
258 
259  TPlusNew=ChooseP(TPlusMin, TPlusMax);
260  Qplus=-(TPlusNew-Ptarget.plus());
261 
262  Qmomentum.setPz( (Qplus-Qminus)/2 );
263  Qmomentum.setE( (Qplus+Qminus)/2 );
264 
265  } while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 ||
266  (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
267 
268  Pprojectile += Qmomentum;
269  Ptarget -= Qmomentum;
270 
271  // Transform back and update SplitableHadron Participant.
272  Pprojectile.transform(toLab);
273  Ptarget.transform(toLab);
274 
275  #ifdef debugDoubleDiffraction
276  G4cout << "Pprojectile after boost to Lab: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
277  G4cout << "Ptarget after boost to Lab: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
278  #endif
279 
280  target->Set4Momentum(Ptarget);
281  projectile->Set4Momentum(Pprojectile);
282 
283  return true;
284 }
285 
286 
288 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
289 {
290  hadron->SplitUp();
291  G4Parton *start= hadron->GetNextParton();
292  if ( start==NULL)
293  { G4cout << " G4QGSDiffractiveExcitation::String() Error:No start parton found"<< G4endl;
294  return NULL;
295  }
296  G4Parton *end = hadron->GetNextParton();
297  if ( end==NULL)
298  { G4cout << " G4QGSDiffractiveExcitation::String() Error:No end parton found"<< G4endl;
299  return NULL;
300  }
301 
302  G4ExcitedString * string;
303  if ( isProjectile )
304  {
305  string= new G4ExcitedString(end,start, +1);
306  } else {
307  string= new G4ExcitedString(start,end, -1);
308  }
309 
310  string->SetPosition(hadron->GetPosition());
311 
312  // momenta of string ends
313 
314  G4double maxAvailMomentumSquared=sqr(hadron->Get4Momentum().mag()/2.);
315 
316  G4double widthOfPtSquare = 0.5*sqr(GeV);
317  G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
318 
319  G4LorentzVector Pstart(G4LorentzVector(pt,0.));
320  G4LorentzVector Pend;
321  Pend.setPx(hadron->Get4Momentum().px() - pt.x());
322  Pend.setPy(hadron->Get4Momentum().py() - pt.y());
323 
324  G4double tm1=hadron->Get4Momentum().minus() +
325  ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
326 
327  G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
328  4. * Pend.perp2() * hadron->Get4Momentum().minus()
329  / hadron->Get4Momentum().plus() ));
330 
331  G4int Sign= isProjectile ? -1 : 1;
332 
333  G4double endMinus = 0.5 * (tm1 + Sign*tm2);
334  G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
335 
336  G4double startPlus= Pstart.perp2() / startMinus;
337  G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
338 
339  Pstart.setPz(0.5*(startPlus - startMinus));
340  Pstart.setE(0.5*(startPlus + startMinus));
341 
342  Pend.setPz(0.5*(endPlus - endMinus));
343  Pend.setE(0.5*(endPlus + endMinus));
344 
345  start->Set4Momentum(Pstart);
346  end->Set4Momentum(Pend);
347 
348  #ifdef debugQGSdiffExictation
349  G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
350  G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
351  G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
352  G4cout << " sum of ends " << Pstart+Pend << G4endl;
353  G4cout << " Original " << hadron->Get4Momentum() << G4endl;
354  #endif
355 
356  return string;
357 }
358 
359 
360 // --------- private methods ----------------------
361 
363 {
364  // choose an x between Xmin and Xmax with P(x) ~ 1/x
365  // to be improved...
366 
367  G4double range=Pmax-Pmin;
368 
369  if ( Pmin <= 0. || range <=0. )
370  {
371  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
372  throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
373  }
374 
375  G4double P;
376  P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmin,G4UniformRand());
377  //debug-hpw cout << "DiffractiveX "<<x<<G4endl;
378  return P;
379 }
380 
381 
383 { // @@ this method is used in FTFModel as well. Should go somewhere common!
384  G4double Pt2;
385 
386  Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.));
387 
388  G4double Pt=std::sqrt(Pt2);
389 
391 
392  return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
393 }