Geant4  10.02.p01
G4ChargeExchange.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: G4ChargeExchange.cc 91897 2015-08-10 09:55:12Z gcosmo $
28 //
29 //
30 // G4 Model: Charge and strangness exchange based on G4LightMedia model
31 // 28 May 2006 V.Ivanchenko
32 //
33 // Modified:
34 // 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
35 // 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled
36 // 12-Jun-12 A.Ribon fix warnings of shadowed variables
37 // 06-Aug-15 A.Ribon migrating to G4Exp, G4Log and G4Pow
38 //
39 
40 #include "G4ChargeExchange.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4ParticleTable.hh"
44 #include "G4ParticleDefinition.hh"
45 #include "G4IonTable.hh"
46 #include "Randomize.hh"
47 #include "G4NucleiProperties.hh"
48 
49 #include "G4Exp.hh"
50 #include "G4Log.hh"
51 #include "G4Pow.hh"
52 
53 
55 {
56  SetMinEnergy( 0.0*GeV );
57  SetMaxEnergy( 100.*TeV );
58 
59  lowEnergyRecoilLimit = 100.*keV;
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) {
108  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
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  G4ParticleDefinition * theDef = 0;
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  G4ParticleDefinition* theRecoil = 0;
145  G4ParticleDefinition* theSecondary = 0;
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) {
250  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
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
269  G4double phi = G4UniformRand()*twopi;
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 = nlv0.e() - m21;
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 > lowEnergyRecoilLimit) {
300  aSec = new G4DynamicParticle(theRecoil, nlv0);
302  } else {
303  if(erec < 0.0) erec = 0.0;
305  }
306  return &theParticleChange;
307 }
308 
310 {
311  G4double aa, bb, cc, dd;
312  if (A <= 62.) {
313  aa = G4Pow::GetInstance()->powA(A, 1.63);
314  bb = 14.5*G4Pow::GetInstance()->powA(A, 0.66);
315  cc = 1.4*G4Pow::GetInstance()->powA(A, 0.33);
316  dd = 10.;
317  } else {
318  aa = G4Pow::GetInstance()->powA(A, 1.33);
319  bb = 60.*G4Pow::GetInstance()->powA(A, 0.33);
320  cc = 0.4*G4Pow::GetInstance()->powA(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 
340  return t;
341 }
342 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4ParticleDefinition * theSMinus
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
static const double MeV
Definition: G4SIunits.hh:211
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4ParticleDefinition * theL
G4ParticleDefinition * theOmega
virtual G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
G4ParticleDefinition * theASPlus
static G4AntiOmegaMinus * AntiOmegaMinus()
G4ParticleDefinition * theXiMinus
CLHEP::Hep3Vector G4ThreeVector
G4double lowEnergyRecoilLimit
G4ParticleDefinition * theNeutron
G4ParticleDefinition * theAOmega
G4ParticleDefinition * theS0
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:491
static G4OmegaMinus * OmegaMinus()
G4ParticleDefinition * theAProton
static G4KaonZeroLong * KaonZeroLong()
G4ParticleDefinition * theAS0
G4ParticleDefinition * theT
G4ParticleDefinition * theA
G4ParticleDefinition * theKMinus
static G4AntiSigmaPlus * AntiSigmaPlus()
G4ParticleDefinition * theK0S
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
static G4SigmaZero * SigmaZero()
Definition: G4SigmaZero.cc:102
G4ParticleDefinition * theAXi0
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
void SetStatusChange(G4HadFinalStateStatus aS)
static ulg bb
Definition: csz_inflate.cc:344
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4XiZero * XiZero()
Definition: G4XiZero.cc:106
void SetMinEnergy(G4double anEnergy)
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4ParticleDefinition * theK0L
double A(double temperature)
const G4ParticleDefinition * GetDefinition() const
static G4KaonZeroShort * KaonZeroShort()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
static const double twopi
Definition: G4SIunits.hh:75
static G4XiMinus * XiMinus()
Definition: G4XiMinus.cc:106
static G4AntiXiMinus * AntiXiMinus()
static G4Triton * Triton()
Definition: G4Triton.cc:95
static G4Proton * Proton()
Definition: G4Proton.cc:93
static G4PionPlus * PionPlus()
Definition: G4PionPlus.cc:98
static const double GeV
Definition: G4SIunits.hh:214
G4ParticleDefinition * theXi0
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4ParticleDefinition * theKPlus
G4ParticleDefinition * theASMinus
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
Definition: G4PionZero.cc:108
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
void SetKineticEnergy(G4double aEnergy)
static const G4double e1
G4ParticleDefinition * theANeutron
static G4SigmaMinus * SigmaMinus()
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
void SetEnergyChange(G4double anEnergy)
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double SampleT(G4double p, G4double A)
static G4AntiLambda * AntiLambda()
G4ParticleDefinition * theD
static G4PionMinus * PionMinus()
Definition: G4PionMinus.cc:98
virtual ~G4ChargeExchange()
static G4AntiSigmaZero * AntiSigmaZero()
const G4double x[NPOINTSGL]
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
G4ParticleDefinition * thePiPlus
void SetMaxEnergy(const G4double anEnergy)
void SetLocalEnergyDeposit(G4double aE)
static G4AntiXiZero * AntiXiZero()
#define G4endl
Definition: G4ios.hh:61
static const double TeV
Definition: G4SIunits.hh:215
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4ParticleDefinition * thePiZero
static const double keV
Definition: G4SIunits.hh:213
static G4SigmaPlus * SigmaPlus()
Definition: G4SigmaPlus.cc:108
G4ParticleDefinition * theAXiMinus
static G4Lambda * Lambda()
Definition: G4Lambda.cc:108
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
static G4KaonPlus * KaonPlus()
Definition: G4KaonPlus.cc:113
static G4He3 * He3()
Definition: G4He3.cc:94
G4ParticleDefinition * theSPlus
void SetMomentumChange(const G4ThreeVector &aV)
static G4AntiNeutron * AntiNeutron()
G4double GetTotalMomentum() const
G4ParticleDefinition * theAntiL
CLHEP::HepLorentzVector G4LorentzVector
G4ParticleDefinition * theProton
G4ParticleDefinition * thePiMinus
G4ParticleDefinition * theHe3