Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4QGSDiffractiveExcitation.cc 107867 2017-12-07 14:44:07Z 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 //============================================================================
64 
65 //#define debugDoubleDiffraction
66 
67 //============================================================================
68 
70 {
71 }
72 
74 {
75 }
76 
77 
79 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target, G4bool ) const // Uzhi Oct. 2016 , G4bool ProjectileDiffraction
80 {
81  #ifdef debugDoubleDiffraction
82  G4cout<<G4endl<<"G4QGSDiffractiveExcitation::ExciteParticipants - Double diffraction."<<G4endl;
83  G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetParticleName()<<" "<<target->GetDefinition()->GetParticleName()<<G4endl;
84  G4cout<<"Proj 4 Mom "<<projectile->Get4Momentum()<<" "<<projectile->Get4Momentum().mag()<<G4endl;
85  G4cout<<"Targ 4 Mom "<<target->Get4Momentum() <<" "<<target->Get4Momentum().mag() <<G4endl;
86  #endif
87 
88 
89 
90  G4LorentzVector Pprojectile=projectile->Get4Momentum();
91 
92  // -------------------- Projectile parameters -----------------------------------
93  G4bool PutOnMassShell=0;
94 
95  G4double M0projectile = Pprojectile.mag(); // Without de-excitation
96 
97  if(M0projectile < projectile->GetDefinition()->GetPDGMass())
98  {
99  PutOnMassShell=1;
100  M0projectile=projectile->GetDefinition()->GetPDGMass();
101  }
102 
103  // -------------------- Target parameters ----------------------------------------------
104  G4LorentzVector Ptarget=target->Get4Momentum();
105 
106  G4double M0target = Ptarget.mag();
107 
108  if(M0target < target->GetDefinition()->GetPDGMass())
109  {
110  PutOnMassShell=1;
111  M0target=target->GetDefinition()->GetPDGMass();
112  }
113 
114  G4LorentzVector Psum=Pprojectile+Ptarget;
115  G4double S=Psum.mag2();
116  G4double SqrtS=std::sqrt(S);
117 
118  if(SqrtS < M0projectile + M0target) {return false;} // The model cannot work for pp-interactions
119  // at Plab < 1.3 GeV/c. Uzhi
120 
121  G4double Mprojectile2 = M0projectile * M0projectile;
122  G4double Mtarget2 = M0target * M0target; //Ptarget.mag2(); // for AA-inter.
123 
124  // Transform momenta to cms and then rotate parallel to z axis;
125 
126  G4LorentzRotation toCms(-1*Psum.boostVector());
127 
128  G4LorentzVector Ptmp=toCms*Pprojectile;
129 
130  if ( Ptmp.pz() <= 0. )
131  {
132  // "String" moving backwards in CMS, abort collision !!
133  //G4cout << " abort Collision!! " << G4endl;
134  return false;
135  }
136 
137  toCms.rotateZ(-1*Ptmp.phi());
138  toCms.rotateY(-1*Ptmp.theta());
139 
140  G4LorentzRotation toLab(toCms.inverse());
141 
142  Pprojectile.transform(toCms);
143  Ptarget.transform(toCms);
144 
145  G4double PZcms2=(S*S+Mprojectile2*Mprojectile2+Mtarget2*Mtarget2-
146  2*S*Mprojectile2-2*S*Mtarget2-2*Mprojectile2*Mtarget2)/4./S;
147 
148  if(PZcms2 < 0) {return false;} // It can be in an interaction with off-shell nuclear nucleon
149 
150  G4double PZcms = std::sqrt(PZcms2);
151 
152  if(PutOnMassShell)
153  {
154  if(Pprojectile.z() > 0.)
155  {
156  Pprojectile.setPz( PZcms);
157  Ptarget.setPz( -PZcms);
158  }
159  else
160  {
161  Pprojectile.setPz(-PZcms);
162  Ptarget.setPz( PZcms);
163  };
164 
165  Pprojectile.setE(std::sqrt(Mprojectile2+sqr(Pprojectile.x())+sqr(Pprojectile.y())+PZcms2));
166  Ptarget.setE( std::sqrt( Mtarget2+sqr( Ptarget.x())+sqr( Ptarget.y())+PZcms2));
167  }
168 
169  G4double maxPtSquare = PZcms2;
170 
171  #ifdef debugDoubleDiffraction
172  G4cout << "Pprojectile after boost to CMS: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
173  G4cout << "Ptarget after boost to CMS: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
174  #endif
175 
176  G4int PrPDGcode=projectile->GetDefinition()->GetPDGEncoding();
177  G4int absPrPDGcode=std::abs(PrPDGcode);
178  G4double MinPrDiffMass(0.);
179  G4double AveragePt2(0.);
180 
181  if(M0projectile <= projectile->GetDefinition()->GetPDGMass())
182  { // Normal projectile
183  if( absPrPDGcode > 1000 ) //------Projectile is baryon --------
184  {
185  MinPrDiffMass = 1.16; // GeV
186  AveragePt2 = 0.3; // GeV^2
187  }
188  else if( absPrPDGcode == 211 || PrPDGcode == 111) //------Projectile is Pion -----------
189  {
190  MinPrDiffMass = 1.0; // GeV
191  AveragePt2 = 0.3; // GeV^2
192  }
193  else if( absPrPDGcode == 321 || absPrPDGcode == 130 || absPrPDGcode == 310) //-Projectile is Kaon-
194  {
195  MinPrDiffMass = 1.1; // GeV
196  AveragePt2 = 0.3; // GeV^2
197  }
198  else //------Projectile is undefined, Nucleon assumed
199  {
200  MinPrDiffMass = 1.16; // GeV
201  AveragePt2 = 0.3; // GeV^2
202  }
203  }
204  else
205  { // Excited projectile
206  MinPrDiffMass = M0projectile + 220.0*MeV;
207  AveragePt2 = 0.3;
208  }
209 
210  MinPrDiffMass = MinPrDiffMass * GeV;
211  AveragePt2 = AveragePt2 * GeV*GeV;
212 //---------------------------------------------
213  G4double MinTrDiffMass = 1.16*GeV;
214 
215  if(SqrtS < MinPrDiffMass + MinTrDiffMass) {return false;} // The model cannot work at low energy
216 
217  G4double MinPrDiffMass2 = MinPrDiffMass * MinPrDiffMass;
218  G4double MinTrDiffMass2 = MinTrDiffMass * MinTrDiffMass;
219 
220  G4double Pt2;
221  G4double ProjMassT2, ProjMassT;
222  G4double TargMassT2, TargMassT;
223  G4double PMinusNew, TPlusNew;
224 
225  G4LorentzVector Qmomentum;
226  G4double Qminus, Qplus;
227 
228  G4int whilecount=0;
229  do {
230  if (whilecount++ >= 500 && (whilecount%100)==0)
231  // G4cout << "G4QGSDiffractiveExcitation::ExciteParticipants possibly looping"
232  // << ", loop count/ maxPtSquare : "
233  // << whilecount << " / " << maxPtSquare << G4endl;
234  if (whilecount > 1000 )
235  {
236  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
237  return false; // Ignore this interaction
238  }
239 
240  // Generate pt
241  Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
242 
243  Pt2=G4ThreeVector(Qmomentum.vect()).mag2();
244  ProjMassT2=MinPrDiffMass2+Pt2;
245  ProjMassT =std::sqrt(ProjMassT2);
246 
247  TargMassT2=MinTrDiffMass2+Pt2;
248  TargMassT =std::sqrt(TargMassT2);
249 
250  if(SqrtS < ProjMassT + TargMassT) continue;
251 
252  PZcms2=(S*S+ProjMassT2*ProjMassT2+
253  TargMassT2*TargMassT2-
254  2.*S*ProjMassT2-2.*S*TargMassT2-
255  2.*ProjMassT2*TargMassT2)/4./S;
256  if(PZcms2 < 0 ) {PZcms2=0;};
257  PZcms =std::sqrt(PZcms2);
258 
259  G4double PMinusMin=std::sqrt(ProjMassT2+PZcms2)-PZcms;
260  G4double PMinusMax=SqrtS-TargMassT;
261 
262  PMinusNew=ChooseP(PMinusMin,PMinusMax);
263  Qminus=PMinusNew-Pprojectile.minus();
264 
265  G4double TPlusMin=std::sqrt(TargMassT2+PZcms2)-PZcms;
266  G4double TPlusMax=SqrtS-ProjMassT;
267 
268  TPlusNew=ChooseP(TPlusMin, TPlusMax);
269  Qplus=-(TPlusNew-Ptarget.plus());
270 
271  Qmomentum.setPz( (Qplus-Qminus)/2 );
272  Qmomentum.setE( (Qplus+Qminus)/2 );
273 
274  } while ( (Pprojectile+Qmomentum).mag2() < MinPrDiffMass2 || // Uzhi No without excitation
275  (Ptarget -Qmomentum).mag2() < MinTrDiffMass2 );
276 
277  Pprojectile += Qmomentum;
278  Ptarget -= Qmomentum;
279 
280  // Transform back and update SplitableHadron Participant.
281  Pprojectile.transform(toLab);
282  Ptarget.transform(toLab);
283 
284  #ifdef debugDoubleDiffraction
285  G4cout << "Pprojectile after boost to Lab: " << Pprojectile <<" "<<Pprojectile.mag()<<G4endl;
286  G4cout << "Ptarget after boost to Lab: " << Ptarget <<" "<<Ptarget.mag() <<G4endl;
287  #endif
288 
289  target->Set4Momentum(Ptarget);
290  projectile->Set4Momentum(Pprojectile);
291 
292  return true;
293 }
294 
295 
297 String(G4VSplitableHadron * hadron, G4bool isProjectile) const
298 {
299  hadron->SplitUp();
300  G4Parton *start= hadron->GetNextParton();
301  if ( start==NULL)
302  { G4cout << " G4QGSDiffractiveExcitation::String() Error:No start parton found"<< G4endl;
303  return NULL;
304  }
305  G4Parton *end = hadron->GetNextParton();
306  if ( end==NULL)
307  { G4cout << " G4QGSDiffractiveExcitation::String() Error:No end parton found"<< G4endl;
308  return NULL;
309  }
310 
311  G4ExcitedString * string;
312  if ( isProjectile )
313  {
314  string= new G4ExcitedString(end,start, +1);
315  } else {
316  string= new G4ExcitedString(start,end, -1);
317  }
318 
319  string->SetPosition(hadron->GetPosition());
320 
321  // momenta of string ends
322 /* // Uzhi 2016
323  G4double ptSquared= hadron->Get4Momentum().perp2();
324  G4double transverseMassSquared= hadron->Get4Momentum().plus()
325  * hadron->Get4Momentum().minus();
326 
327 
328  G4double maxAvailMomentumSquared=
329  sqr( std::sqrt(transverseMassSquared) - std::sqrt(ptSquared) );
330 */
331  G4double maxAvailMomentumSquared=sqr(hadron->Get4Momentum().mag()/2.); // Uzhi 2016
332 
333  G4double widthOfPtSquare = 0.5*sqr(GeV); //0.25; // Uzhi 2016 // Uzhi <Pt^2>=0.25 ??????????????????
334  G4ThreeVector pt=GaussianPt(widthOfPtSquare,maxAvailMomentumSquared);
335 
336  G4LorentzVector Pstart(G4LorentzVector(pt,0.));
337  G4LorentzVector Pend;
338  Pend.setPx(hadron->Get4Momentum().px() - pt.x());
339  Pend.setPy(hadron->Get4Momentum().py() - pt.y());
340 
341  G4double tm1=hadron->Get4Momentum().minus() +
342  ( Pend.perp2()-Pstart.perp2() ) / hadron->Get4Momentum().plus();
343 
344  G4double tm2= std::sqrt( std::max(0., sqr(tm1) -
345  4. * Pend.perp2() * hadron->Get4Momentum().minus()
346  / hadron->Get4Momentum().plus() ));
347 
348  G4int Sign= isProjectile ? -1 : 1;
349 
350  G4double endMinus = 0.5 * (tm1 + Sign*tm2);
351  G4double startMinus= hadron->Get4Momentum().minus() - endMinus;
352 
353  G4double startPlus= Pstart.perp2() / startMinus;
354  G4double endPlus = hadron->Get4Momentum().plus() - startPlus;
355 
356  Pstart.setPz(0.5*(startPlus - startMinus));
357  Pstart.setE(0.5*(startPlus + startMinus));
358 
359  Pend.setPz(0.5*(endPlus - endMinus));
360  Pend.setE(0.5*(endPlus + endMinus));
361 
362  start->Set4Momentum(Pstart);
363  end->Set4Momentum(Pend);
364 
365  #ifdef debugQGSdiffExictation
366  G4cout << " generated string flavors " << start->GetPDGcode() << " / " << end->GetPDGcode() << G4endl;
367  G4cout << " generated string momenta: quark " << start->Get4Momentum() << "mass : " <<start->Get4Momentum().mag()<< G4endl;
368  G4cout << " generated string momenta: Diquark " << end ->Get4Momentum() << "mass : " <<end->Get4Momentum().mag()<< G4endl;
369  G4cout << " sum of ends " << Pstart+Pend << G4endl;
370  G4cout << " Original " << hadron->Get4Momentum() << G4endl;
371  #endif
372 
373  return string;
374 }
375 
376 
377 // --------- private methods ----------------------
378 
380 {
381  // choose an x between Xmin and Xmax with P(x) ~ 1/x
382  // to be improved...
383 
384  G4double range=Pmax-Pmin;
385 
386  if ( Pmin <= 0. || range <=0. )
387  {
388  G4cout << " Pmin, range : " << Pmin << " , " << range << G4endl;
389  throw G4HadronicException(__FILE__, __LINE__, "G4QGSDiffractiveExcitation::ChooseP : Invalid arguments ");
390  }
391 
392  G4double P;
393  P=Pmin * G4Pow::GetInstance()->powA(Pmax/Pmin,G4UniformRand());
394  //debug-hpw cout << "DiffractiveX "<<x<<G4endl;
395  return P;
396 }
397 
399 { // @@ this method is used in FTFModel as well. Should go somewhere common!
400 
401  G4double Pt2;
402 
403  Pt2 = -AveragePt2 * G4Log(1. + G4UniformRand() * (G4Exp(-maxPtSquare/AveragePt2)-1.));
404 
405  G4double Pt=std::sqrt(Pt2);
406 
407  G4double phi=G4UniformRand() * twopi;
408 
409  return G4ThreeVector (Pt*std::cos(phi), Pt*std::sin(phi), 0.);
410 }
double mag2() const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4ParticleDefinition * GetDefinition() const
HepLorentzVector & transform(const HepRotation &)
CLHEP::Hep3Vector G4ThreeVector
double plus() const
const XML_Char * target
Definition: expat.h:268
static constexpr double MeV
Definition: G4SIunits.hh:214
double py() const
#define G4endl
Definition: G4ios.hh:61
double px() const
const G4String & GetParticleName() const
G4double ChooseP(G4double Pmin, G4double Pmax) const
virtual G4Parton * GetNextParton()=0
G4double GetPDGMass() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double S(double temp)
double perp2() const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
const G4ThreeVector & GetPosition() const
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
const G4ParticleDefinition const G4Material *G4double range
double minus() const
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
TMarker * pt
Definition: egs.C:25
HepLorentzVector & rotateY(double)
const G4ThreeVector const G4double const
virtual G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction=TRUE) const
G4ThreeVector GaussianPt(G4double AveragePt2, G4double maxPtSquare) const
virtual void SplitUp()=0
static double P[]
HepLorentzVector & rotateZ(double)
G4int GetPDGcode() const
Definition: G4Parton.hh:127
double mag() const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
double x() const
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
Hep3Vector vect() const
T sqr(const T &x)
Definition: templates.hh:145
CLHEP::HepLorentzVector G4LorentzVector
double y() const
const G4LorentzVector & Get4Momentum() const
Definition: G4Parton.hh:143
static constexpr double GeV
Definition: G4SIunits.hh:217
virtual G4ExcitedString * String(G4VSplitableHadron *aHadron, G4bool isProjectile) const
HepLorentzRotation & transform(const HepBoost &b)
void Set4Momentum(const G4LorentzVector &a4Momentum)
void Set4Momentum(const G4LorentzVector &aMomentum)
Definition: G4Parton.hh:148