Geant4  10.02.p02
G4SingleDiffractiveExcitation.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: G4SingleDiffractiveExcitation.cc 93563 2015-10-26 14:46:09Z 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 
51 
53 :
54 widthOfPtSquare(-2*sqr(sigmaPt)) , minExtraMass(minextraMass),
55 minmass(x0mass)
56 {}
57 
59 {}
60 
63 {
64 
65  G4LorentzVector Pprojectile=projectile->Get4Momentum();
66  G4double Mprojectile2=sqr(projectile->GetDefinition()->GetPDGMass() + minExtraMass);
67 
68  G4LorentzVector Ptarget=target->Get4Momentum();
69  G4double Mtarget2=sqr(target->GetDefinition()->GetPDGMass() + minExtraMass);
70  // G4cout << "E proj, target :" << Pprojectile.e() << ", " <<
71  // Ptarget.e() << G4endl;
72 
73  G4bool KeepProjectile= G4UniformRand() > 0.5;
74 
75  // reset the min.mass of the non diffractive particle to its value, ( minus a bit for rounding...)
76  if ( KeepProjectile )
77  {
78  // cout << " Projectile fix" << G4endl;
79  Mprojectile2 = sqr(projectile->GetDefinition()->GetPDGMass() * (1-perCent) );
80  } else {
81  // cout << " Target fix" << G4endl;
82  Mtarget2=sqr(target->GetDefinition()->GetPDGMass() * (1-perCent) );
83  }
84 
85  // Transform momenta to cms and then rotate parallel to z axis;
86 
87  G4LorentzVector Psum;
88  Psum=Pprojectile+Ptarget;
89 
90  G4LorentzRotation toCms(-1*Psum.boostVector());
91 
92  G4LorentzVector Ptmp=toCms*Pprojectile;
93 
94  if ( Ptmp.pz() <= 0. )
95  {
96  // "String" moving backwards in CMS, abort collision !!
97  // G4cout << " abort Collision!! " << G4endl;
98  return false;
99  }
100 
101  toCms.rotateZ(-1*Ptmp.phi());
102  toCms.rotateY(-1*Ptmp.theta());
103 
104  // G4cout << "Pprojectile be4 boost " << Pprojectile << G4endl;
105  // G4cout << "Ptarget be4 boost : " << Ptarget << G4endl;
106 
107 
108 
109  G4LorentzRotation toLab(toCms.inverse());
110 
111  Pprojectile.transform(toCms);
112  Ptarget.transform(toCms);
113 
114  G4LorentzVector Qmomentum;
115  G4int whilecount=0;
116  do {
117  // Generate pt
118 
119  G4double maxPtSquare=sqr(Ptarget.pz());
120  if (whilecount++ >= 500 && (whilecount%100)==0)
121  // G4cout << "G4SingleDiffractiveExcitation::ExciteParticipants possibly looping"
122  // << ", loop count/ maxPtSquare : "
123  // << whilecount << " / " << maxPtSquare << G4endl;
124  if (whilecount > 1000 )
125  {
126  Qmomentum=G4LorentzVector(0.,0.,0.,0.);
127  // G4cout << "G4SingleDiffractiveExcitation::ExciteParticipants: Aborting loop!" << G4endl;
128  return false; // Ignore this interaction
129  }
130  Qmomentum=G4LorentzVector(GaussianPt(widthOfPtSquare,maxPtSquare),0);
131 
132 
133  // Momentum transfer
134  G4double Xmin = minmass / ( Pprojectile.e() + Ptarget.e() );
135  G4double Xmax=1.;
136  G4double Xplus =ChooseX(Xmin,Xmax);
137  G4double Xminus=ChooseX(Xmin,Xmax);
138 
139  G4double pt2=G4ThreeVector(Qmomentum.vect()).mag2();
140  G4double Qplus =-1 * pt2 / Xminus/Ptarget.minus();
141  G4double Qminus= pt2 / Xplus /Pprojectile.plus();
142 
143  if ( KeepProjectile )
144  {
145  Qminus = (sqr(projectile->GetDefinition()->GetPDGMass()) + pt2 )
146  / (Pprojectile.plus() + Qplus )
147  - Pprojectile.minus();
148  } else
149  {
150  Qplus = Ptarget.plus()
151  - (sqr(target->GetDefinition()->GetPDGMass()) + pt2 )
152  / (Ptarget.minus() - Qminus );
153  }
154 
155  Qmomentum.setPz( (Qplus-Qminus)/2 );
156  Qmomentum.setE( (Qplus+Qminus)/2 );
157 
158  // G4cout << "Qplus / Qminus " << Qplus << " / " << Qminus<<G4endl;
159  // G4cout << "pt2 " << pt2 << G4endl;
160  // G4cout << "Qmomentum " << Qmomentum << G4endl;
161  // G4cout << " Masses (P/T) : " << (Pprojectile+Qmomentum).mag() <<
162  // " / " << (Ptarget-Qmomentum).mag() << G4endl;
163 
164  } while ( (Ptarget-Qmomentum).mag2() <= Mtarget2 /* Loop checking, 26.10.2015, A.Ribon */
165  || (Pprojectile+Qmomentum).mag2() <= Mprojectile2
166  || (Ptarget-Qmomentum).e() < 0.
167  || (Pprojectile+Qmomentum).e() < 0. );
168 
169 
170  // G4double Ecms=Pprojectile.e() + Ptarget.e();
171 
172  Pprojectile += Qmomentum;
173 
174  Ptarget -= Qmomentum;
175 
176  // G4cout << "Pprojectile.e() : " << Pprojectile.e() << G4endl;
177  // G4cout << "Ptarget.e() : " << Ptarget.e() << G4endl;
178 
179  // G4cout << "end event_______________________________________________"<<G4endl;
180  //
181 
182 
183  // G4cout << "Pprojectile with Q : " << Pprojectile << G4endl;
184  // G4cout << "Ptarget with Q : " << Ptarget << G4endl;
185  // G4cout << "Projectile back: " << toLab * Pprojectile << G4endl;
186  // G4cout << "Target back: " << toLab * Ptarget << G4endl;
187 
188  // Transform back and update SplitableHadron Participant.
189  Pprojectile.transform(toLab);
190  Ptarget.transform(toLab);
191 
192  // G4cout << "G4SingleDiffractiveExcitation- Target mass " << Ptarget.mag() << G4endl;
193  // G4cout << "G4SingleDiffractiveExcitation- Projectile mass " << Pprojectile.mag() << G4endl;
194 
195  target->Set4Momentum(Ptarget);
196  projectile->Set4Momentum(Pprojectile);
197 
198 
199  return true;
200 }
201 
202 
203 
204 
205 // --------- private methods ----------------------
206 
208 {
209  // choose an x between Xmin and Xmax with P(x) ~ 1/x
210 
211  // to be improved...
212 
213  G4double range=Xmax-Xmin;
214 
215  if ( Xmin <= 0. || range <=0. )
216  {
217  G4cout << " Xmin, range : " << Xmin << " , " << range << G4endl;
218  throw G4HadronicException(__FILE__, __LINE__, "G4SingleDiffractiveExcitation::ChooseX : Invalid arguments ");
219  }
220 
221  G4double x;
222  do {
223  x=Xmin + G4UniformRand() * range;
224  } while ( Xmin/x < G4UniformRand() ); /* Loop checking, 26.10.2015, A.Ribon */
225 
226  // cout << "DiffractiveX "<<x<<G4endl;
227  return x;
228 }
229 
231 { // @@ this method is used in FTFModel as well. Should go somewhere common!
232 
233  G4double pt2;
234 
235  const G4int maxNumberOfLoops = 1000;
236  G4int loopCounter = -1;
237  do {
238  pt2=widthSquare * G4Log( G4UniformRand() );
239  } while ( ( pt2 > maxPtSquare) && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 26.10.2015, A.Ribon */
240  if ( loopCounter >= maxNumberOfLoops ) pt2 = 0.0;
241 
242  pt2=std::sqrt(pt2);
243 
244  G4double phi=G4UniformRand() * twopi;
245 
246  return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);
247 }
248 
249 
250 
251 
252 
G4bool ExciteParticipants(G4VSplitableHadron *aPartner, G4VSplitableHadron *bPartner) const
G4double ChooseX(G4double Xmin, G4double Xmax) const
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
int G4int
Definition: G4Types.hh:78
const G4ParticleDefinition * GetDefinition() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
static const double twopi
Definition: G4SIunits.hh:75
static const double perCent
Definition: G4SIunits.hh:329
const G4LorentzVector & Get4Momentum() const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double GetPDGMass() const
G4SingleDiffractiveExcitation(G4double sigmaPt=0.6 *CLHEP::GeV, G4double minExtraMass=250 *CLHEP::MeV, G4double x0mass=250 *CLHEP::MeV)
G4ThreeVector GaussianPt(G4double widthSquare, G4double maxPtSquare) const
const G4double x[NPOINTSGL]
#define G4endl
Definition: G4ios.hh:61
void Set4Momentum(const G4LorentzVector &a4Momentum)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector