Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/parton_string/qgsm/src/G4QuarkExchange.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // ------------------------------------------------------------
 28 //      GEANT 4 class implemetation file
 29 //
 30 //      ---------------- G4QuarkExchange --------------
 31 //             by V. Uzhinsky, October 2016.
 32 //       QuarkExchange is used by strings models.
 33 //      Take a projectile and a target.
 34 //Simulate Q exchange with excitation of projectile or target.
 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 2020
 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 June 2020 : suppression of strange quark pair prodution,
 58                                      //                  i.e. u:d:s=1:1:0.04 . Need to be tuned!
 59 }
 60 
 61 G4QuarkExchange::~G4QuarkExchange(){}
 62 
 63 G4bool G4QuarkExchange::
 64 ExciteParticipants(G4VSplitableHadron *projectile, G4VSplitableHadron *target) const
 65 {
 66   #ifdef debugQuarkExchange
 67     G4cout<<G4endl<<"G4QuarkExchange::ExciteParticipants"<<G4endl;
 68   #endif
 69 
 70   G4LorentzVector Pprojectile = projectile->Get4Momentum();
 71   G4double Mprojectile        = projectile->GetDefinition()->GetPDGMass();
 72   G4double Mprojectile2       = sqr(Mprojectile);
 73 
 74   G4LorentzVector Ptarget = target->Get4Momentum();
 75   G4double Mtarget        = target->GetDefinition()->GetPDGMass();
 76   G4double Mtarget2       = sqr(Mtarget);
 77 
 78   #ifdef debugQuarkExchange
 79     G4cout<<"Proj Targ "<<projectile->GetDefinition()->GetPDGEncoding()<<" "<<target->GetDefinition()->GetPDGEncoding()<<G4endl;
 80     G4cout<<"Proj. 4-Mom "<<Pprojectile<<" "<<Pprojectile.mag()<<G4endl
 81           <<"Targ. 4-Mom "<<Ptarget    <<" "<<Ptarget.mag()   <<G4endl;
 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-Mtarget "<<SqrtS<<" "<<Mprojectile<<" "<<Mtarget
 90           <<" "<<SqrtS-Mprojectile-Mtarget<<G4endl;
 91   #endif
 92   if (SqrtS-Mprojectile-Mtarget <= 250.0*MeV) {
 93     #ifdef debugQuarkExchange
 94       G4cerr<<"Energy is too small for quark exchange!"<<G4endl;
 95       G4cerr<<"Projectile: "<<projectile->GetDefinition()->GetPDGEncoding()<<" "
 96             <<Pprojectile<<" "<<Pprojectile.mag()<<G4endl;
 97       G4cerr<<"Target:     "<<target->GetDefinition()->GetPDGEncoding()<<" "
 98             <<Ptarget<<" "<<Ptarget.mag()<<G4endl; 
 99       G4cerr<<"sqrt(S) = "<<SqrtS<<" Mp + Mt = "<<Pprojectile.mag()+Ptarget.mag()<<G4endl;
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, abort collision !!
111     //         G4cout << " abort Collision!! " << G4endl;
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 " << Pprojectile << G4endl;
125     G4cout << "Ptarget      in CMS " << Ptarget     << G4endl;
126   #endif
127   G4double maxPtSquare=sqr(Ptarget.pz());
128 
129   G4double ProjectileMinDiffrMass = Pprojectile.mag()/GeV;
130   G4double TargetMinDiffrMass     = Ptarget.mag()/GeV;
131 
132   G4double AveragePt2(0.);
133 
134   G4int    PDGcode=projectile->GetDefinition()->GetPDGEncoding();
135   G4int absPDGcode=std::abs(PDGcode); 
136 
137   G4bool ProjectileDiffraction = true;
138 
139   // Also for heavy hadrons, assume 50% probability of projectile diffraction.
140   if ( absPDGcode > 1000 )                          { ProjectileDiffraction = G4UniformRand() <= 0.5; }
141   if ( (absPDGcode == 211) || (absPDGcode == 111) ) { ProjectileDiffraction = G4UniformRand() <= 0.66; }
142   if ( (absPDGcode == 321) || (absPDGcode == 311)  || 
143        (   PDGcode == 130) || (   PDGcode == 310) ) { ProjectileDiffraction = G4UniformRand() <= 0.5; }
144   if ( absPDGcode > 400  &&  absPDGcode < 600 )     { ProjectileDiffraction = G4UniformRand() <= 0.5; }
145 
146   //G4cout<<"ProjectileDiffr "<<ProjectileDiffraction<<G4endl;
147 
148   if ( ProjectileDiffraction ) {
149     if ( absPDGcode > 1000 )                            //------Projectile is baryon --------
150     {
151       if ( absPDGcode > 4000 && absPDGcode < 6000 )  // Projectile is a charm or bottom baryon
152       {
153         ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25;  // GeV
154         AveragePt2 = 0.3;                                                                      // GeV^2
155       }
156       else
157       {
158         ProjectileMinDiffrMass = 1.16;              // GeV
159         AveragePt2 = 0.3;                           // GeV^2
160       }
161     }
162     else if( absPDGcode == 211 || absPDGcode ==  111) //------Projectile is Pion -----------
163     {
164       ProjectileMinDiffrMass = 1.0;               // GeV
165       AveragePt2 = 0.3;                           // GeV^2
166     }
167     else if( absPDGcode == 321 || absPDGcode == 130 || absPDGcode == 310) //Projectile is Kaon
168     {
169       ProjectileMinDiffrMass = 1.1;               // GeV
170       AveragePt2 = 0.3;                           // GeV^2
171     }
172     else if( absPDGcode == 22)                        //------Projectile is Gamma -----------
173     {
174       ProjectileMinDiffrMass = 0.25;             // GeV
175       AveragePt2 = 0.36;                         // GeV^2
176     }
177     else if( absPDGcode > 400 && absPDGcode < 600)  // Projectile is a charm or bottom meson
178     {
179       ProjectileMinDiffrMass = projectile->GetDefinition()->GetPDGMass()/CLHEP::GeV + 0.25;  // GeV
180       AveragePt2 = 0.3;                                                                      // GeV^2
181     }
182     else                                             //------Projectile is undefined, Nucleon assumed
183     {
184       ProjectileMinDiffrMass = 1.1;              // GeV
185       AveragePt2 = 0.3;                          // GeV^2
186     };
187 
188     ProjectileMinDiffrMass = ProjectileMinDiffrMass * GeV;
189     Mprojectile2=sqr(ProjectileMinDiffrMass);
190 
191     if (G4UniformRand() <= 0.5) TargetMinDiffrMass += 0.22; 
192     TargetMinDiffrMass *= GeV;
193     Mtarget2 = sqr( TargetMinDiffrMass) ;
194   }
195   else
196   {
197     if (G4UniformRand() <= 0.5) ProjectileMinDiffrMass += 0.22; 
198     ProjectileMinDiffrMass *=GeV;
199     Mprojectile2=sqr(ProjectileMinDiffrMass);
200 
201     TargetMinDiffrMass = 1.16*GeV;                     // For target nucleon
202     Mtarget2 = sqr( TargetMinDiffrMass) ;
203     AveragePt2 = 0.3;                                  // GeV^2
204   }   // end of if ( ProjectileDiffraction )
205 
206   AveragePt2 = AveragePt2 * GeV*GeV;
207 
208   if ( SqrtS - (ProjectileMinDiffrMass+TargetMinDiffrMass) < 220.0*MeV ) return false;
209 
210   //----------------------- 
211   G4double Pt2, PZcms, PZcms2;
212   G4double ProjMassT2, ProjMassT;
213   G4double TargMassT2, TargMassT;
214   G4double PMinusMin, PMinusMax,  sqrtPMinusMin, sqrtPMinusMax;
215   G4double TPlusMin, TPlusMax,    sqrtTPlusMin,  sqrtTPlusMax;
216   G4double PMinusNew, PPlusNew, TPlusNew(0.), TMinusNew;
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 interaction
230     }
231     
232     // Generate pt
233     Qmomentum=G4LorentzVector(GaussianPt(AveragePt2,maxPtSquare),0);
234 
235     Pt2 = G4ThreeVector( Qmomentum.vect() ).mag2();
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  TargMassT  SqrtS  S  ProjectileDiffraction"<<G4endl;
243       G4cout<<whilecount<<" "<<Pt2<<" "<<ProjMassT<<" "<<TargMassT<<" "<<SqrtS<<" "<<S<<" "<<ProjectileDiffraction<<G4endl;
244     #endif
245 
246     if ( SqrtS < ProjMassT + TargMassT + 220.0*MeV ) continue;
247 
248     PZcms2 = ( S*S + ProjMassT2*ProjMassT2 + TargMassT2*TargMassT2
249                - 2.0*S*ProjMassT2 - 2.0*S*TargMassT2 - 2.0*ProjMassT2*TargMassT2 ) / 4.0 / S;
250 
251     if ( PZcms2 < 0 ) continue;
252 
253     PZcms = std::sqrt( PZcms2 );
254 
255     if ( ProjectileDiffraction )
256     { // The projectile will fragment, the target will saved.
257       PMinusMin = std::sqrt( ProjMassT2 + PZcms2 ) - PZcms;
258       PMinusMax = SqrtS - TargMassT;
259       sqrtPMinusMin = std::sqrt(PMinusMin); sqrtPMinusMax = std::sqrt(PMinusMax);
260 
261       if ( absPDGcode > 1000 ) 
262       { 
263   PMinusNew = PMinusMax * (1.0 - (1.0 - PMinusMin/PMinusMax) 
264                                      * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
265       } 
266       else if ( (absPDGcode == 211) || (absPDGcode == 111) ) 
267       {
268         while (true)
269         {
270           x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusMin)*G4UniformRand();
271     y=G4UniformRand();
272           if ( y < 1.0-0.7 * x/sqrtPMinusMax ) break;
273   }
274   PMinusNew = sqr(x);
275       } 
276       else if ( (absPDGcode == 321) || (absPDGcode == 311)  || 
277           (   PDGcode == 130) || (   PDGcode == 310) ) 
278       {  // For K-mesons it must be found !!! Uzhi
279         while (true)
280         {
281     x=sqrtPMinusMax-(sqrtPMinusMax-sqrtPMinusMin)*G4UniformRand();
282     y=G4UniformRand();
283           if ( y < 1.0-0.7 * x/sqrtPMinusMax ) break;
284         }
285   PMinusNew = sqr(x);
286       } 
287       else
288       { 
289   PMinusNew = PMinusMax * (1.0 - (1.0 - PMinusMin/PMinusMax) 
290                                  * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
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 projectile will saved.
302       TPlusMin = std::sqrt( TargMassT2 + PZcms2 ) - PZcms;
303       TPlusMax = SqrtS - ProjMassT;
304       sqrtTPlusMin = std::sqrt(TPlusMin); sqrtTPlusMax = std::sqrt(TPlusMax);
305 
306       if ( absPDGcode > 1000 ) 
307       { 
308         TPlusNew = TPlusMax * (1.0 - (1.0 - TPlusMin/TPlusMax) 
309                                     * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
310       } 
311       else if ( (absPDGcode == 211) || (absPDGcode == 111) ) 
312       {
313         while (true)
314         {
315           x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPlusMin)*G4UniformRand();
316           y=G4UniformRand();
317           if ( y < 1.0-0.7 * x/sqrtTPlusMax ) break;
318         }
319   TPlusNew = sqr(x);
320       } 
321       else if ( (absPDGcode == 321) || (absPDGcode == 311)  || 
322           (   PDGcode == 130) || (   PDGcode == 310) ) 
323       { // For K-mesons it must be found !!! Uzhi
324   while (true)
325   {
326           x=sqrtTPlusMax-(sqrtTPlusMax-sqrtTPlusMin)*G4UniformRand();
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 - TPlusMin/TPlusMax) 
334                                     * G4Pow::GetInstance()->powA(G4UniformRand(),0.3333) );
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 (Pprojectile + Qmomentum).mag2()  Mprojectile2"<<G4endl;
349       G4cout<<ProjectileDiffraction<<" "<<( Pprojectile + Qmomentum ).mag2()<<" "<< Mprojectile2<<G4endl;
350       G4cout<<"TargetDiffraction     (Ptarget     - Qmomentum).mag2()  Mtarget2"<<G4endl;
351       G4cout<<!ProjectileDiffraction<<" "<<( Ptarget    - Qmomentum ).mag2()<<" "<< Mtarget2<<G4endl;
352     #endif
353 
354   } while ( ( ProjectileDiffraction&&( Pprojectile + Qmomentum ).mag2() <  Mprojectile2 ) ||
355             (!ProjectileDiffraction&&( Ptarget     - Qmomentum ).mag2() <  Mtarget2       )   ); 
356     // Repeat the sampling because there was not any excitation
357 
358   Pprojectile += Qmomentum;
359 
360   Ptarget     -= Qmomentum;
361 
362   // Transform back and update SplitableHadron Participant.
363   Pprojectile.transform(toLab);
364   Ptarget.transform(toLab);
365 
366   #ifdef debugQuarkExchange
367     G4cout << "Pprojectile  in Lab. " << Pprojectile << G4endl;
368     G4cout << "Ptarget      in Lab. " << Ptarget     << G4endl;
369     G4cout << "G4QuarkExchange: Projectile mass  " <<  Pprojectile.mag() << G4endl;
370     G4cout << "G4QuarkExchange: Target mass      " <<  Ptarget.mag() << G4endl;
371   #endif
372 
373   target->Set4Momentum(Ptarget);
374   projectile->Set4Momentum(Pprojectile);
375 
376   //=================================== Quark exchange ================================
377   projectile->SplitUp();
378   target->SplitUp();
379 
380   G4Parton* PrQuark = nullptr;
381   G4Parton* TrQuark = nullptr;
382 
383   if( projectile->GetDefinition()->GetBaryonNumber() >= 0 ) {  // Uzhi June 2020
384     // Quark exchange ----
385     PrQuark = projectile->GetNextParton(); 
386     TrQuark = target->GetNextParton(); 
387     G4ParticleDefinition * Tmp = PrQuark->GetDefinition();
388     PrQuark->SetDefinition(TrQuark->GetDefinition());
389     TrQuark->SetDefinition(Tmp);
390     return true;
391   }
392 
393   // Quark exchage for projectile anti-baryon (annihilation and new Q pair creation ---
394   // This part added by Uzhi June 2020
395   PrQuark = projectile->GetNextAntiParton();
396   TrQuark = target->GetNextParton(); 
397   if( -PrQuark->GetDefinition()->GetPDGEncoding() == TrQuark->GetDefinition()->GetPDGEncoding() ){
398     G4int QuarkCode = 1 + (int)(G4UniformRand()/StrangeSuppress);
399     G4ParticleDefinition* AQpr = G4ParticleTable::GetParticleTable()->FindParticle(-QuarkCode);
400     G4ParticleDefinition*  Qtr = G4ParticleTable::GetParticleTable()->FindParticle( QuarkCode);
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(G4double widthSquare, G4double maxPtSquare) const
414 {            //  @@ this method is used in FTFModel as well. Should go somewhere common!
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) && ++loopCounter < maxNumberOfLoops );  /* Loop checking, 07.08.2015, A.Ribon */
423   if ( loopCounter >= maxNumberOfLoops ) {
424     pt2 = 0.99*maxPtSquare;  // Just an acceptable value, without any physics consideration. 
425   }
426 
427   pt2=std::sqrt(pt2);
428 
429   G4double phi=G4UniformRand() * twopi;
430 
431   return G4ThreeVector (pt2*std::cos(phi), pt2*std::sin(phi), 0.);    
432 }
433 
434