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