| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 // GEANT 4 class implemetation file
29 //
30 // ---------------- G4SingleDiffractiveEx
31 // by Gunter Folger, October 1998.
32 // diffractive Excitation used by strings
33 // Take a projectile and a target
34 // excite the projectile and target
35 // -------------------------------------------
36
37 #include "G4SingleDiffractiveExcitation.hh"
38 #include "globals.hh"
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "Randomize.hh"
42 #include "G4LorentzRotation.hh"
43 #include "G4ThreeVector.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4VSplitableHadron.hh"
46 #include "G4ExcitedString.hh"
47
48 #include "G4Log.hh"
49 #include "G4Pow.hh"
50
51 //#define debugSingleDiffraction
52
53 G4SingleDiffractiveExcitation::G4SingleDiffrac
54
55 G4SingleDiffractiveExcitation::~G4SingleDiffra
56
57 G4bool G4SingleDiffractiveExcitation::
58 ExciteParticipants( G4VSplitableHadron *projec
59 G4bool ProjectileDiffracti
60 {
61 #ifdef debugSingleDiffraction
62 G4cout<<G4endl<<"G4SingleDiffractiveExcita
63 #endif
64
65 G4LorentzVector Pprojectile=projectile->Get4
66 G4double Mprojectile = projectile->GetDef
67 G4double Mprojectile2=sqr(projectile->GetDef
68
69 G4LorentzVector Ptarget=target->Get4Momentum
70 G4double Mtarget = target->GetDefinition(
71 G4double Mtarget2=sqr(target->GetDefinition(
72
73 #ifdef debugSingleDiffraction
74 G4cout<<"Proj Targ "<<projectile->GetDefin
75 G4cout<<"Pr Tr 4-Mom "<<Pprojectile<<" "<<
76 <<" "<<Ptarget <<" "<<
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 "<<Sqrt
85 <<" "<<SqrtS-Mprojectile-Mtarget<<G4
86 #endif
87 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV)
88 #ifdef debugSingleDiffraction
89 G4cerr<<"Projectile: "<<projectile->GetDef
90 <<Pprojectile<<" "<<Pprojectile.mag(
91 G4cerr<<"Target: "<<target->GetDefinit
92 <<Ptarget<<" "<<Ptarget.mag()<<G4end
93 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "
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, abor
105 // G4cout << " abort Collision!! "
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 " << Pproje
118 G4cout << "Ptarget in CMS " << Ptarge
119 #endif
120 G4double maxPtSquare=sqr(Ptarget.pz());
121
122 G4double ProjectileMinDiffrMass(0.), TargetM
123 G4double AveragePt2(0.);
124 G4int absPDGcode=std::abs(projectile->GetDef
125
126 if ( ProjectileDiffraction ) {
127 if ( absPDGcode > 1000 )
128 {
129 if ( absPDGcode > 4000 && absPDGcode < 6
130 {
131 ProjectileMinDiffrMass = projectile->G
132 AveragePt2 = 0.3;
133 }
134 else
135 {
136 ProjectileMinDiffrMass = 1.16;
137 AveragePt2 = 0.3;
138 }
139 }
140 else if( absPDGcode == 211 || absPDGcode =
141 {
142 ProjectileMinDiffrMass = 1.0;
143 AveragePt2 = 0.3;
144 }
145 else if( absPDGcode == 321 || absPDGcode =
146 {
147 ProjectileMinDiffrMass = 1.1;
148 AveragePt2 = 0.3;
149 }
150 else if( absPDGcode == 22)
151 {
152 ProjectileMinDiffrMass = 0.25;
153 AveragePt2 = 0.36;
154 }
155 else if( absPDGcode > 400 && absPDGcode <
156 {
157 ProjectileMinDiffrMass = projectile->Get
158 AveragePt2 = 0.3;
159 }
160 else
161 {
162 ProjectileMinDiffrMass = 1.1;
163 AveragePt2 = 0.3;
164 };
165
166 ProjectileMinDiffrMass = ProjectileMinDiff
167 Mprojectile2=sqr(ProjectileMinDiffrMass);
168 }
169 else
170 {
171 TargetMinDiffrMass = 1.16*GeV;
172 Mtarget2 = sqr( TargetMinDiffrMass) ;
173 AveragePt2 = 0.3;
174 } // end of if ( ProjectileDiffraction )
175
176 AveragePt2 = AveragePt2 * GeV*GeV;
177
178 G4double Pt2, PZcms, PZcms2;
179 G4double ProjMassT2, ProjMassT;
180 G4double TargMassT2, TargMassT;
181 G4double PMinusMin, PMinusMax;
182 G4double TPlusMin, TPlusMax;
183 G4double PMinusNew, PPlusNew, TPlusNew, TMin
184
185 G4LorentzVector Qmomentum;
186 G4double Qminus, Qplus;
187
188 G4int whilecount=0;
189 do {
190 whilecount++;
191
192 if (whilecount > 1000 )
193 {
194 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
195 return false; // Ignore this intera
196 }
197
198 // Generate pt
199 Qmomentum=G4LorentzVector(GaussianPt(Avera
200
201 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma
202
203 ProjMassT2 = Mprojectile2 + Pt2;
204 ProjMassT = std::sqrt( ProjMassT2 );
205 TargMassT2 = Mtarget2 + Pt2;
206 TargMassT = std::sqrt( TargMassT2 );
207
208 #ifdef debugSingleDiffraction
209 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjM
210 #endif
211 if ( SqrtS < ProjMassT + TargMassT ) conti
212
213 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + T
214 - 2.0*S*ProjMassT2 - 2.0*S*TargM
215
216 if ( PZcms2 < 0 ) continue;
217
218 PZcms = std::sqrt( PZcms2 );
219
220 if ( ProjectileDiffraction )
221 { // The projectile will fragment, t
222 PMinusMin = std::sqrt( ProjMassT2 + PZcm
223 PMinusMax = SqrtS - TargMassT;
224
225 PMinusNew = ChooseX( PMinusMin, PMinusMa
226 TMinusNew = SqrtS - PMinusNew;
227
228 Qminus = Ptarget.minus() - TMinusNew;
229 TPlusNew = TargMassT2 / TMinusNew;
230 Qplus = Ptarget.plus() - TPlusNew;
231
232 } else { // The target will fragment, the
233 TPlusMin = std::sqrt( TargMassT2 + PZcms
234 TPlusMax = SqrtS - ProjMassT;
235
236 TPlusNew = ChooseX( TPlusMin, TPlusMax )
237 PPlusNew = SqrtS - TPlusNew;
238
239 Qplus = PPlusNew - Pprojectile.plus();
240 PMinusNew = ProjMassT2 / PPlusNew;
241 Qminus = PMinusNew - Pprojectile.minus()
242 }
243
244 Qmomentum.setPz( (Qplus - Qminus)/2 );
245 Qmomentum.setE( (Qplus + Qminus)/2 );
246
247 #ifdef debugSingleDiffraction
248 G4cout<<ProjectileDiffraction<<" "<<( Pp
249 G4cout<<!ProjectileDiffraction<<" "<<( P
250 #endif
251
252 } while ( ( ProjectileDiffraction&&( Pprojec
253 (!ProjectileDiffraction&&( Ptarget
254 // Repeat the sampling because there was n
255
256 Pprojectile += Qmomentum;
257
258 Ptarget -= Qmomentum;
259
260 // Transform back and update SplitableHadron
261 Pprojectile.transform(toLab);
262 Ptarget.transform(toLab);
263
264 #ifdef debugSingleDiffraction
265 G4cout << "Pprojectile in Lab. " << Pproj
266 G4cout << "Ptarget in Lab. " << Ptarg
267 G4cout << "G4SingleDiffractiveExcitation-
268 G4cout << "G4SingleDiffractiveExcitation-
269 #endif
270
271 target->Set4Momentum(Ptarget);
272 projectile->Set4Momentum(Pprojectile);
273
274 return true;
275 }
276
277 // --------- private methods -----------------
278
279 G4double G4SingleDiffractiveExcitation::Choose
280 {
281 // choose an x between Xmin and Xmax with P(
282 G4double range=Xmax-Xmin;
283
284 if ( Xmin <= 0. || range <=0. )
285 {
286 G4cout << " Xmin, range : " << Xmin << " ,
287 throw G4HadronicException(__FILE__, __LINE
288 }
289
290 G4double x = Xmin*G4Pow::GetInstance()->powA
291 // G4double x = 1.0/sqr(1.0/std::sqrt(Xmin)
292 return x;
293 }
294
295
296 G4ThreeVector G4SingleDiffractiveExcitation::G
297 { // @@ this method is used in FTF
298
299 G4double pt2;
300
301 const G4int maxNumberOfLoops = 1000;
302 G4int loopCounter = 0;
303 do {
304 pt2=-widthSquare * G4Log( G4UniformRand()
305 } while ( ( pt2 > maxPtSquare) && ++loopCoun
306 if ( loopCounter >= maxNumberOfLoops ) {
307 pt2 = 0.99*maxPtSquare; // Just an accept
308 }
309
310 pt2=std::sqrt(pt2);
311
312 G4double phi=G4UniformRand() * twopi;
313
314 return G4ThreeVector (pt2*std::cos(phi), pt2
315 }
316
317