| 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 // ---------------- G4QuarkExchange -----
31 // by V. Uzhinsky, October 2016.
32 // QuarkExchange is used by strings mode
33 // Take a projectile and a target.
34 //Simulate Q exchange with excitation of proje
35 // -------------------------------------------
36
37 #include "G4QuarkExchange.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 "G4ParticleTable.hh" // Uzhi June 20
49
50 #include "G4Log.hh"
51 #include "G4Pow.hh"
52
53 //#define debugQuarkExchange
54
55 G4QuarkExchange::G4QuarkExchange()
56 {
57 StrangeSuppress = (1.0-0.04)/2.0; // Uzhi J
58 //
59 }
60
61 G4QuarkExchange::~G4QuarkExchange(){}
62
63 G4bool G4QuarkExchange::
64 ExciteParticipants(G4VSplitableHadron *project
65 {
66 #ifdef debugQuarkExchange
67 G4cout<<G4endl<<"G4QuarkExchange::ExcitePa
68 #endif
69
70 G4LorentzVector Pprojectile = projectile->Ge
71 G4double Mprojectile = projectile->Ge
72 G4double Mprojectile2 = sqr(Mprojectil
73
74 G4LorentzVector Ptarget = target->Get4Moment
75 G4double Mtarget = target->GetDefinit
76 G4double Mtarget2 = sqr(Mtarget);
77
78 #ifdef debugQuarkExchange
79 G4cout<<"Proj Targ "<<projectile->GetDefin
80 G4cout<<"Proj. 4-Mom "<<Pprojectile<<" "<<
81 <<"Targ. 4-Mom "<<Ptarget <<" "<<
82 #endif
83
84 G4LorentzVector Psum=Pprojectile+Ptarget;
85 G4double SqrtS=Psum.mag();
86 G4double S =Psum.mag2();
87
88 #ifdef debugQuarkExchange
89 G4cout<<"SS Mpr Mtr SqrtS-Mprojectile-Mtar
90 <<" "<<SqrtS-Mprojectile-Mtarget<<G4
91 #endif
92 if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV)
93 #ifdef debugQuarkExchange
94 G4cerr<<"Energy is too small for quark e
95 G4cerr<<"Projectile: "<<projectile->GetD
96 <<Pprojectile<<" "<<Pprojectile.ma
97 G4cerr<<"Target: "<<target->GetDefin
98 <<Ptarget<<" "<<Ptarget.mag()<<G4e
99 G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt =
100 #endif
101 return true;
102 }
103
104 G4LorentzRotation toCms(-1*Psum.boostVector(
105
106 G4LorentzVector Ptmp=toCms*Pprojectile;
107
108 if ( Ptmp.pz() <= 0. )
109 {
110 // "String" moving backwards in CMS, abor
111 // G4cout << " abort Collision!! "
112 return false;
113 }
114
115 toCms.rotateZ(-1*Ptmp.phi());
116 toCms.rotateY(-1*Ptmp.theta());
117
118 G4LorentzRotation toLab(toCms.inverse());
119
120 Pprojectile.transform(toCms);
121 Ptarget.transform(toCms);
122
123 #ifdef debugQuarkExchange
124 G4cout << "Pprojectile in CMS " << Pproje
125 G4cout << "Ptarget in CMS " << Ptarge
126 #endif
127 G4double maxPtSquare=sqr(Ptarget.pz());
128
129 G4double ProjectileMinDiffrMass = Pprojectil
130 G4double TargetMinDiffrMass = Ptarget.ma
131
132 G4double AveragePt2(0.);
133
134 G4int PDGcode=projectile->GetDefinition()
135 G4int absPDGcode=std::abs(PDGcode);
136
137 G4bool ProjectileDiffraction = true;
138
139 // Also for heavy hadrons, assume 50% probab
140 if ( absPDGcode > 1000 )
141 if ( (absPDGcode == 211) || (absPDGcode == 1
142 if ( (absPDGcode == 321) || (absPDGcode == 3
143 ( PDGcode == 130) || ( PDGcode == 3
144 if ( absPDGcode > 400 && absPDGcode < 600
145
146 //G4cout<<"ProjectileDiffr "<<ProjectileDiff
147
148 if ( ProjectileDiffraction ) {
149 if ( absPDGcode > 1000 )
150 {
151 if ( absPDGcode > 4000 && absPDGcode < 6
152 {
153 ProjectileMinDiffrMass = projectile->G
154 AveragePt2 = 0.3;
155 }
156 else
157 {
158 ProjectileMinDiffrMass = 1.16;
159 AveragePt2 = 0.3;
160 }
161 }
162 else if( absPDGcode == 211 || absPDGcode =
163 {
164 ProjectileMinDiffrMass = 1.0;
165 AveragePt2 = 0.3;
166 }
167 else if( absPDGcode == 321 || absPDGcode =
168 {
169 ProjectileMinDiffrMass = 1.1;
170 AveragePt2 = 0.3;
171 }
172 else if( absPDGcode == 22)
173 {
174 ProjectileMinDiffrMass = 0.25;
175 AveragePt2 = 0.36;
176 }
177 else if( absPDGcode > 400 && absPDGcode <
178 {
179 ProjectileMinDiffrMass = projectile->Get
180 AveragePt2 = 0.3;
181 }
182 else
183 {
184 ProjectileMinDiffrMass = 1.1;
185 AveragePt2 = 0.3;
186 };
187
188 ProjectileMinDiffrMass = ProjectileMinDiff
189 Mprojectile2=sqr(ProjectileMinDiffrMass);
190
191 if (G4UniformRand() <= 0.5) TargetMinDiffr
192 TargetMinDiffrMass *= GeV;
193 Mtarget2 = sqr( TargetMinDiffrMass) ;
194 }
195 else
196 {
197 if (G4UniformRand() <= 0.5) ProjectileMinD
198 ProjectileMinDiffrMass *=GeV;
199 Mprojectile2=sqr(ProjectileMinDiffrMass);
200
201 TargetMinDiffrMass = 1.16*GeV;
202 Mtarget2 = sqr( TargetMinDiffrMass) ;
203 AveragePt2 = 0.3;
204 } // end of if ( ProjectileDiffraction )
205
206 AveragePt2 = AveragePt2 * GeV*GeV;
207
208 if ( SqrtS - (ProjectileMinDiffrMass+TargetM
209
210 //-----------------------
211 G4double Pt2, PZcms, PZcms2;
212 G4double ProjMassT2, ProjMassT;
213 G4double TargMassT2, TargMassT;
214 G4double PMinusMin, PMinusMax, sqrtPMinusMi
215 G4double TPlusMin, TPlusMax, sqrtTPlusMin
216 G4double PMinusNew, PPlusNew, TPlusNew(0.),
217
218 G4LorentzVector Qmomentum;
219 G4double Qminus, Qplus;
220
221 G4double x(0.), y(0.);
222 G4int whilecount=0;
223 do {
224 whilecount++;
225
226 if (whilecount > 1000 )
227 {
228 Qmomentum=G4LorentzVector(0.,0.,0.,0.);
229 return false; // Ignore this intera
230 }
231
232 // Generate pt
233 Qmomentum=G4LorentzVector(GaussianPt(Avera
234
235 Pt2 = G4ThreeVector( Qmomentum.vect() ).ma
236 ProjMassT2 = Mprojectile2 + Pt2;
237 ProjMassT = std::sqrt( ProjMassT2 );
238 TargMassT2 = Mtarget2 + Pt2;
239 TargMassT = std::sqrt( TargMassT2 );
240
241 #ifdef debugQuarkExchange
242 G4cout<<"whilecount Pt2 ProjMassT Tar
243 G4cout<<whilecount<<" "<<Pt2<<" "<<ProjM
244 #endif
245
246 if ( SqrtS < ProjMassT + TargMassT + 220.0
247
248 PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + T
249 - 2.0*S*ProjMassT2 - 2.0*S*Targ
250
251 if ( PZcms2 < 0 ) continue;
252
253 PZcms = std::sqrt( PZcms2 );
254
255 if ( ProjectileDiffraction )
256 { // The projectile will fragment, the tar
257 PMinusMin = std::sqrt( ProjMassT2 + PZcm
258 PMinusMax = SqrtS - TargMassT;
259 sqrtPMinusMin = std::sqrt(PMinusMin); sq
260
261 if ( absPDGcode > 1000 )
262 {
263 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinus
264 * G4Pow::
265 }
266 else if ( (absPDGcode == 211) || (absPDG
267 {
268 while (true)
269 {
270 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtP
271 y=G4UniformRand();
272 if ( y < 1.0-0.7 * x/sqrtPMinusMax )
273 }
274 PMinusNew = sqr(x);
275 }
276 else if ( (absPDGcode == 321) || (absPDG
277 ( PDGcode == 130) || ( PDGcode =
278 { // For K-mesons it must be found !!!
279 while (true)
280 {
281 x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusM
282 y=G4UniformRand();
283 if ( y < 1.0-0.7 * x/sqrtPMinusMax )
284 }
285 PMinusNew = sqr(x);
286 }
287 else
288 {
289 PMinusNew = PMinusMax * (1.0 - (1.0 - PMinus
290 * G4Pow::GetI
291 };
292
293 TMinusNew = SqrtS - PMinusNew;
294
295 Qminus = Ptarget.minus() - TMinusNew;
296 TPlusNew = TargMassT2 / TMinusNew;
297 Qplus = Ptarget.plus() - TPlusNew;
298
299 }
300 else
301 { // The target will fragment, the project
302 TPlusMin = std::sqrt( TargMassT2 + PZcms
303 TPlusMax = SqrtS - ProjMassT;
304 sqrtTPlusMin = std::sqrt(TPlusMin); sqrt
305
306 if ( absPDGcode > 1000 )
307 {
308 TPlusNew = TPlusMax * (1.0 - (1.0 - TP
309 * G4Pow::G
310 }
311 else if ( (absPDGcode == 211) || (absPDG
312 {
313 while (true)
314 {
315 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPl
316 y=G4UniformRand();
317 if ( y < 1.0-0.7 * x/sqrtTPlusMax )
318 }
319 TPlusNew = sqr(x);
320 }
321 else if ( (absPDGcode == 321) || (absPDG
322 ( PDGcode == 130) || ( PDGcode =
323 { // For K-mesons it must be found !!! U
324 while (true)
325 {
326 x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPl
327 y=G4UniformRand();
328 if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break;
329 }
330 }
331 else
332 {
333 TPlusNew = TPlusMax * (1.0 - (1.0 - TP
334 * G4Pow::G
335 };
336
337 PPlusNew = SqrtS - TPlusNew;
338
339 Qplus = PPlusNew - Pprojectile.plus();
340 PMinusNew = ProjMassT2 / PPlusNew;
341 Qminus = PMinusNew - Pprojectile.minus()
342 }
343
344 Qmomentum.setPz( (Qplus - Qminus)/2 );
345 Qmomentum.setE( (Qplus + Qminus)/2 );
346
347 #ifdef debugQuarkExchange
348 G4cout<<"ProjectileDiffraction (Pproject
349 G4cout<<ProjectileDiffraction<<" "<<( Pp
350 G4cout<<"TargetDiffraction (Ptarget
351 G4cout<<!ProjectileDiffraction<<" "<<( P
352 #endif
353
354 } while ( ( ProjectileDiffraction&&( Pprojec
355 (!ProjectileDiffraction&&( Ptarget
356 // Repeat the sampling because there was n
357
358 Pprojectile += Qmomentum;
359
360 Ptarget -= Qmomentum;
361
362 // Transform back and update SplitableHadron
363 Pprojectile.transform(toLab);
364 Ptarget.transform(toLab);
365
366 #ifdef debugQuarkExchange
367 G4cout << "Pprojectile in Lab. " << Pproj
368 G4cout << "Ptarget in Lab. " << Ptarg
369 G4cout << "G4QuarkExchange: Projectile mas
370 G4cout << "G4QuarkExchange: Target mass
371 #endif
372
373 target->Set4Momentum(Ptarget);
374 projectile->Set4Momentum(Pprojectile);
375
376 //=================================== Quark
377 projectile->SplitUp();
378 target->SplitUp();
379
380 G4Parton* PrQuark = nullptr;
381 G4Parton* TrQuark = nullptr;
382
383 if( projectile->GetDefinition()->GetBaryonNu
384 // Quark exchange ----
385 PrQuark = projectile->GetNextParton();
386 TrQuark = target->GetNextParton();
387 G4ParticleDefinition * Tmp = PrQuark->GetD
388 PrQuark->SetDefinition(TrQuark->GetDefinit
389 TrQuark->SetDefinition(Tmp);
390 return true;
391 }
392
393 // Quark exchage for projectile anti-baryon
394 // This part added by Uzhi June 2020
395 PrQuark = projectile->GetNextAntiParton();
396 TrQuark = target->GetNextParton();
397 if( -PrQuark->GetDefinition()->GetPDGEncodin
398 G4int QuarkCode = 1 + (int)(G4UniformRand(
399 G4ParticleDefinition* AQpr = G4ParticleTab
400 G4ParticleDefinition* Qtr = G4ParticleTab
401 if( (AQpr != nullptr) && (Qtr != nullptr)
402 PrQuark->SetDefinition(AQpr);
403 TrQuark->SetDefinition( Qtr);
404 }
405 }
406
407 return true;
408 }
409
410
411 // --------- private methods -----------------
412
413 G4ThreeVector G4QuarkExchange::GaussianPt(G4do
414 { // @@ this method is used in FTF
415
416 G4double pt2;
417
418 const G4int maxNumberOfLoops = 1000;
419 G4int loopCounter = 0;
420 do {
421 pt2=-widthSquare * G4Log( G4UniformRand()
422 } while ( ( pt2 > maxPtSquare) && ++loopCoun
423 if ( loopCounter >= maxNumberOfLoops ) {
424 pt2 = 0.99*maxPtSquare; // Just an accept
425 }
426
427 pt2=std::sqrt(pt2);
428
429 G4double phi=G4UniformRand() * twopi;
430
431 return G4ThreeVector (pt2*std::cos(phi), pt2
432 }
433
434