ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ChargeExchange.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ChargeExchange.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 //
28 //
29 // G4 Model: Charge and strangness exchange based on G4LightMedia model
30 // 28 May 2006 V.Ivanchenko
31 //
32 // Modified:
33 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
34 // 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled
35 // 12-Jun-12 A.Ribon fix warnings of shadowed variables
36 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log and G4Pow
37 //
38 
39 #include "G4ChargeExchange.hh"
40 #include "G4PhysicalConstants.hh"
41 #include "G4SystemOfUnits.hh"
42 #include "G4ParticleTable.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4IonTable.hh"
45 #include "Randomize.hh"
46 #include "G4NucleiProperties.hh"
47 
48 #include "G4Exp.hh"
49 #include "G4Log.hh"
50 #include "G4Pow.hh"
51 
52 #include "G4HadronicParameters.hh"
53 
54 
56 {
57  SetMinEnergy( 0.0*GeV );
59 
60  lowestEnergyLimit = 1.*MeV;
61 
89  theA = G4Alpha::Alpha();
90  theHe3 = G4He3::He3();
91 }
92 
94 {}
95 
97  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
98 {
100  const G4HadProjectile* aParticle = &aTrack;
101  G4double ekin = aParticle->GetKineticEnergy();
102 
103  G4int A = targetNucleus.GetA_asInt();
104  G4int Z = targetNucleus.GetZ_asInt();
105 
106  if(ekin <= lowestEnergyLimit || A < 3) {
109  return &theParticleChange;
110  }
111 
112  G4double plab = aParticle->GetTotalMomentum();
113 
114  if (verboseLevel > 1)
115  G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
116  << plab/GeV << " GeV/c "
117  << " ekin(MeV) = " << ekin/MeV << " "
118  << aParticle->GetDefinition()->GetParticleName() << G4endl;
119 
120  // Scattered particle referred to axis of incident particle
121  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
122 
123  G4int N = A - Z;
124  G4int projPDG = theParticle->GetPDGEncoding();
125  if (verboseLevel > 1)
126  G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
127  << " PDGcode= " << projPDG << " on nucleus Z= " << Z
128  << " A= " << A << " N= " << N
129  << G4endl;
130 
131  const G4ParticleDefinition* theDef = nullptr;
132 
134  G4LorentzVector lv1 = aParticle->Get4Momentum();
135  G4LorentzVector lv0(0.0,0.0,0.0,mass2);
136 
137  G4LorentzVector lv = lv0 + lv1;
138  G4ThreeVector bst = lv.boostVector();
139  lv1.boost(-bst);
140  lv0.boost(-bst);
141 
142  // Sample final particles
143  G4bool theHyperon = false;
144  const G4ParticleDefinition* theRecoil = nullptr;
145  const G4ParticleDefinition* theSecondary = nullptr;
146 
147  if(theParticle == theProton) {
148  theSecondary = theNeutron;
149  Z++;
150  } else if(theParticle == theNeutron) {
151  theSecondary = theProton;
152  Z--;
153  } else if(theParticle == thePiPlus) {
154  theSecondary = thePiZero;
155  Z++;
156  } else if(theParticle == thePiMinus) {
157  theSecondary = thePiZero;
158  Z--;
159  } else if(theParticle == theKPlus) {
160  if(G4UniformRand()<0.5) theSecondary = theK0S;
161  else theSecondary = theK0L;
162  Z++;
163  } else if(theParticle == theKMinus) {
164  if(G4UniformRand()<0.5) theSecondary = theK0S;
165  else theSecondary = theK0L;
166  Z--;
167  } else if(theParticle == theK0S || theParticle == theK0L) {
168  if(G4UniformRand()*A < G4double(Z)) {
169  theSecondary = theKPlus;
170  Z--;
171  } else {
172  theSecondary = theKMinus;
173  Z++;
174  }
175  } else if(theParticle == theANeutron) {
176  theSecondary = theAProton;
177  Z++;
178  } else if(theParticle == theAProton) {
179  theSecondary = theANeutron;
180  Z--;
181  } else if(theParticle == theL) {
183  if(G4UniformRand()*A < G4double(Z)) {
184  if(x < 0.2) {
185  theSecondary = theS0;
186  } else if (x < 0.4) {
187  theSecondary = theSPlus;
188  Z--;
189  } else if (x < 0.6) {
190  theSecondary = theProton;
191  theRecoil = theL;
192  theHyperon = true;
193  A--;
194  } else if (x < 0.8) {
195  theSecondary = theProton;
196  theRecoil = theS0;
197  theHyperon = true;
198  A--;
199  } else {
200  theSecondary = theNeutron;
201  theRecoil = theSPlus;
202  theHyperon = true;
203  A--;
204  }
205  } else {
206  if(x < 0.2) {
207  theSecondary = theS0;
208  } else if (x < 0.4) {
209  theSecondary = theSMinus;
210  Z++;
211  } else if (x < 0.6) {
212  theSecondary = theNeutron;
213  theRecoil = theL;
214  A--;
215  theHyperon = true;
216  } else if (x < 0.8) {
217  theSecondary = theNeutron;
218  theRecoil = theS0;
219  theHyperon = true;
220  A--;
221  } else {
222  theSecondary = theProton;
223  theRecoil = theSMinus;
224  theHyperon = true;
225  A--;
226  }
227  }
228  }
229 
230  if (Z == 1 && A == 2) theDef = theD;
231  else if (Z == 1 && A == 3) theDef = theT;
232  else if (Z == 2 && A == 3) theDef = theHe3;
233  else if (Z == 2 && A == 4) theDef = theA;
234  else {
235  theDef =
237  }
238  if(!theSecondary) { return &theParticleChange; }
239 
240  G4double m11 = theSecondary->GetPDGMass();
241  G4double m21 = theDef->GetPDGMass();
242  if(theRecoil) { m21 += theRecoil->GetPDGMass(); }
243  else { theRecoil = theDef; }
244 
245  G4double etot = lv0.e() + lv1.e();
246 
247  // kinematiacally impossible
248  if(etot < m11 + m21) {
251  return &theParticleChange;
252  }
253 
254  G4ThreeVector p1 = lv1.vect();
255  G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
256  // G4double e2 = etot - e1;
257  G4double ptot = std::sqrt(e1*e1 - m11*m11);
258 
259  G4double tmax = 4.0*ptot*ptot;
260  G4double g2 = GeV*GeV;
261 
262  G4double t = g2*SampleT(tmax/g2, A);
263 
264  if(verboseLevel>1) {
265  G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
266  << " ptot= " << ptot << G4endl;
267  }
268  // Sampling in CM system
270  G4double cost = 1. - 2.0*t/tmax;
271  if(std::abs(cost) > 1.0) cost = 1.0;
272  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
273 
274  //if (verboseLevel > 1)
275  // G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
276 
277  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
278  v1 *= ptot;
279  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
280  G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
281 
282  nlv0.boost(bst);
283  nlv1.boost(bst);
284 
287  G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
289 
290  G4double erec = std::max(nlv0.e() - m21, 0.0);
291 
292  //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl;
293 
294  if(theHyperon) {
296  aSec = new G4DynamicParticle();
297  aSec->SetDefinition(theRecoil);
298  aSec->SetKineticEnergy(0.0);
299  } else if(erec > GetRecoilEnergyThreshold()) {
300  aSec = new G4DynamicParticle(theRecoil, nlv0);
302  } else {
304  }
305  return &theParticleChange;
306 }
307 
309 {
310  G4double aa, bb, cc, dd;
311  G4Pow* g4pow = G4Pow::GetInstance();
312  if (A <= 62.) {
313  aa = g4pow->powZ(A, 1.63);
314  bb = 14.5*g4pow->powZ(A, 0.66);
315  cc = 1.4*g4pow->powZ(A, 0.33);
316  dd = 10.;
317  } else {
318  aa = g4pow->powZ(A, 1.33);
319  bb = 60.*g4pow->powZ(A, 0.33);
320  cc = 0.4*g4pow->powZ(A, 0.40);
321  dd = 10.;
322  }
323  G4double x1 = (1.0 - G4Exp(-tmax*bb))*aa/bb;
324  G4double x2 = (1.0 - G4Exp(-tmax*dd))*cc/dd;
325 
326  G4double t;
327  G4double y = bb;
328  if(G4UniformRand()*(x1 + x2) < x2) y = dd;
329 
330  const G4int maxNumberOfLoops = 10000;
331  G4int loopCounter = 0;
332  do {
333  t = -G4Log(G4UniformRand())/y;
334  } while ( (t > tmax) &&
335  ++loopCounter < maxNumberOfLoops ); /* Loop checking, 10.08.2015, A.Ribon */
336  if ( loopCounter >= maxNumberOfLoops ) {
337  t = 0.0;
338  }
339  return t;
340 }
341