Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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 // $Id: G4SingleDiffractiveExcitation.cc 108191 2018-01-18 16:10:21Z gcosmo $
28 // ------------------------------------------------------------
29 // GEANT 4 class implemetation file
30 //
31 // ---------------- G4SingleDiffractiveExcitation --------------
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 // ------------------------------------------------------------
37 
39 #include "globals.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "Randomize.hh"
43 #include "G4LorentzRotation.hh"
44 #include "G4ThreeVector.hh"
45 #include "G4ParticleDefinition.hh"
46 #include "G4VSplitableHadron.hh"
47 #include "G4ExcitedString.hh"
48 
49 #include "G4Log.hh"
50 #include "G4Pow.hh"
51 
52 //#define debugSingleDiffraction
53 
55 
57 
59 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target, 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 //----------------------- Uzhi Oct. 2016 Start
123  G4double ProjectileMinDiffrMass(0.), TargetMinDiffrMass(0.);
124  G4double AveragePt2(0.);
125  G4int absPDGcode=std::abs(projectile->GetDefinition()->GetPDGEncoding());
126 
127  if ( ProjectileDiffraction ) {
128  if( absPDGcode > 1000 ) //------Projectile is baryon --------
129  {
130  ProjectileMinDiffrMass = 1.16; // GeV
131  AveragePt2 = 0.3; // GeV^2
132  }
133  else if( absPDGcode == 211 || absPDGcode == 111) //------Projectile is Pion -----------
134  {
135  ProjectileMinDiffrMass = 1.0; // GeV
136  AveragePt2 = 0.3; // GeV^2
137  }
138  else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon
139  {
140  ProjectileMinDiffrMass = 1.1; // GeV
141  AveragePt2 = 0.3; // GeV^2
142  }
143  else if( absPDGcode == 22) //------Projectile is Gamma -----------
144  {
145  ProjectileMinDiffrMass = 0.25; // GeV
146  AveragePt2 = 0.36; // GeV^2
147  }
148  else //------Projectile is undefined, Nucleon assumed
149  {
150  ProjectileMinDiffrMass = 1.1; // GeV
151  AveragePt2 = 0.3; // GeV^2
152  };
153 
154  ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV;
155  Mprojectile2=sqr(ProjectileMinDiffrMass);
156  }
157  else
158  {
159  TargetMinDiffrMass = 1.16*GeV; // For target nucleon
160  Mtarget2 = sqr( TargetMinDiffrMass) ;
161  AveragePt2 = 0.3; // GeV^2
162  } // end of if ( ProjectileDiffraction )
163 
164  AveragePt2 = AveragePt2 * GeV*GeV; // Uzhi 6 Oct. 2016
165 //----------------------- Uzhi Oct. 2016 End
166 
167  G4double Pt2, PZcms, PZcms2;
168  G4double ProjMassT2, ProjMassT;
169  G4double TargMassT2, TargMassT;
170  G4double PMinusMin, PMinusMax;
171  //G4double PPlusMin , PPlusMax;
172  G4double TPlusMin, TPlusMax;
173  G4double PMinusNew, PPlusNew, TPlusNew, TMinusNew;
174 
175  G4LorentzVector Qmomentum;
176  G4double Qminus, Qplus;
177 
178  G4int whilecount=0;
179  do {
180  whilecount++;
181 
182  if (whilecount > 1000 )
183  {
184  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
185  return false; // Ignore this interaction
186  }
187 
188  // Generate pt
189  Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
190 
191  Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
192 
193  ProjMassT2 = Mprojectile2 + Pt2;
194  ProjMassT = std::sqrt( ProjMassT2 );
195  TargMassT2 = Mtarget2 + Pt2;
196  TargMassT = std::sqrt( TargMassT2 );
197 
198  #ifdef debugSingleDiffraction
199  G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl;
200  #endif
201  if ( SqrtS < ProjMassT + TargMassT ) continue;
202 
203  PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
204  - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
205 
206  if ( PZcms2 < 0 ) continue;
207 
208  PZcms = std::sqrt( PZcms2 );
209 
210  if ( ProjectileDiffraction )
211  { // The projectile will fragment, the target will saved.
212  PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
213  PMinusMax = SqrtS - TargMassT;
214 
215  PMinusNew = ChooseX( PMinusMin, PMinusMax );
216  TMinusNew = SqrtS - PMinusNew;
217 
218  Qminus = Ptarget.minus() - TMinusNew;
219  TPlusNew = TargMassT2 / TMinusNew;
220  Qplus = Ptarget.plus() - TPlusNew;
221 
222  } else {// The target will fragment, the projectile will saved.
223  TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
224  TPlusMax = SqrtS - ProjMassT;
225 
226  TPlusNew = ChooseX( TPlusMin, TPlusMax ); // TPlusMax; //
227  PPlusNew = SqrtS - TPlusNew;
228 
229  Qplus = PPlusNew - Pprojectile.plus();
230  PMinusNew = ProjMassT2 / PPlusNew;
231  Qminus = PMinusNew - Pprojectile.minus();
232  }
233 
234 
235  Qmomentum.setPz( (Qplus - Qminus)/2 );
236  Qmomentum.setE( (Qplus + Qminus)/2 );
237 
238  #ifdef debugSingleDiffraction
239  G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl;
240  G4cout<<!ProjectileDiffraction<<" "<<( Ptarget - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl;
241  #endif
242 
243  } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() < Mprojectile2 ) ||
244  (!ProjectileDiffraction&&( Ptarget - Qmomentum ).mag2() < Mtarget2 ) );
245  // Repeat the sampling because there was not any excitation
246 
247 
248  Pprojectile += Qmomentum;
249 
250  Ptarget -= Qmomentum;
251 
252  // Transform back and update SplitableHadron Participant.
253  Pprojectile.transform(toLab);
254  Ptarget.transform(toLab);
255 
256  #ifdef debugSingleDiffraction
257  G4cout << "Pprojectile in Lab. " << Pprojectile << G4endl;
258  G4cout << "Ptarget in Lab. " << Ptarget << G4endl;
259  G4cout << "G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() << G4endl;
260  G4cout << "G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() << G4endl;
261  #endif
262 
263 //G4int Uzhi; G4cin>>Uzhi;
264 
265  target->Set4Momentum(Ptarget);
266  projectile->Set4Momentum(Pprojectile);
267 
268  return true;
269 }
270 
271 // --------- private methods ----------------------
272 
274 {
275  // choose an x between Xmin and Xmax with P(x) ~ 1/x
276  G4double range=Xmax-Xmin;
277 
278  if ( Xmin <= 0. || range <=0. )
279  {
280  G4cout << " Xmin, range : " << Xmin << " , " << range << G4endl;
281  throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation::ChooseX : Invalid arguments ");
282  }
283 
284  G4double x = Xmin*G4Pow::GetInstance()->powA(Xmax/Xmin, G4UniformRand() );
285 // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin) - G4UniformRand() * (1.0/std::sqrt(Xmin) - 1.0/std::sqrt(Xmax)));
286  return x;
287 }
288 
290 { // @@ this method is used in FTFModel as well. Should go somewhere common!
291 
292  G4double pt2;
293 
294  const G4int maxNumberOfLoops = 1000;
295  G4int loopCounter = 0;
296  do {
297  pt2=-widthSquare * G4Log( G4UniformRand() );
298  } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 07.08.2015, A.Ribon */
299  if ( loopCounter >= maxNumberOfLoops ) {
300  pt2 = 0.99*maxPtSquare; // Just an acceptable value, without any physics consideration.
301  }
302 
303  pt2=std::sqrt(pt2);
304 
305  G4double phi=G4UniformRand() * twopi;
306 
307  return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);
308 }
309 
310 
311 
312 
313 
Float_t x
Definition: compare.C:6
double mag2() const
const G4ParticleDefinition * GetDefinition() const
HepLorentzVector & transform(const HepRotation &)
CLHEP::Hep3Vector G4ThreeVector
const XML_Char * target
Definition: expat.h:268
static constexpr double MeV
Definition: G4SIunits.hh:214
#define G4endl
Definition: G4ios.hh:61
G4double GetPDGMass() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double S(double temp)
int Xmax
G4double ChooseX(G4double Xmin, G4double Xmax) const
double G4double
Definition: G4Types.hh:76
bool G4bool
Definition: G4Types.hh:79
static G4Pow * GetInstance()
Definition: G4Pow.cc:57
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:242
const G4ParticleDefinition const G4Material *G4double range
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner, G4bool ProjectileDiffraction) const
#define G4UniformRand()
Definition: Randomize.hh:53
static constexpr double twopi
Definition: G4SIunits.hh:76
HepLorentzVector & rotateY(double)
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const
HepLorentzVector & rotateZ(double)
double mag() const
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
const G4LorentzVector & Get4Momentum() const
Hep3Vector vect() const
T sqr(const T &x)
Definition: templates.hh:145
CLHEP::HepLorentzVector G4LorentzVector
int Xmin
static constexpr double GeV
Definition: G4SIunits.hh:217
HepLorentzRotation & transform(const HepBoost &b)
void Set4Momentum(const G4LorentzVector &a4Momentum)