Geant4  10.00.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 66892 2013-01-17 10:57:59Z gunter $
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 //
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 
49 {
50  SetMinEnergy( 0.0*GeV );
51  SetMaxEnergy( 100.*TeV );
52 
53  lowEnergyRecoilLimit = 100.*keV;
54  lowestEnergyLimit = 1.*MeV;
55 
83  theA = G4Alpha::Alpha();
84  theHe3 = G4He3::He3();
85 }
86 
88 {}
89 
91  const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
92 {
94  const G4HadProjectile* aParticle = &aTrack;
95  G4double ekin = aParticle->GetKineticEnergy();
96 
97  G4int A = targetNucleus.GetA_asInt();
98  G4int Z = targetNucleus.GetZ_asInt();
99 
100  if(ekin <= lowestEnergyLimit || A < 3) {
102  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
103  return &theParticleChange;
104  }
105 
106  G4double plab = aParticle->GetTotalMomentum();
107 
108  if (verboseLevel > 1)
109  G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
110  << plab/GeV << " GeV/c "
111  << " ekin(MeV) = " << ekin/MeV << " "
112  << aParticle->GetDefinition()->GetParticleName() << G4endl;
113 
114  // Scattered particle referred to axis of incident particle
115  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
116 
117  G4int N = A - Z;
118  G4int projPDG = theParticle->GetPDGEncoding();
119  if (verboseLevel > 1)
120  G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
121  << " PDGcode= " << projPDG << " on nucleus Z= " << Z
122  << " A= " << A << " N= " << N
123  << G4endl;
124 
125  G4ParticleDefinition * theDef = 0;
126 
128  G4LorentzVector lv1 = aParticle->Get4Momentum();
129  G4LorentzVector lv0(0.0,0.0,0.0,mass2);
130 
131  G4LorentzVector lv = lv0 + lv1;
132  G4ThreeVector bst = lv.boostVector();
133  lv1.boost(-bst);
134  lv0.boost(-bst);
135 
136  // Sample final particles
137  G4bool theHyperon = false;
138  G4ParticleDefinition* theRecoil = 0;
139  G4ParticleDefinition* theSecondary = 0;
140 
141  if(theParticle == theProton) {
142  theSecondary = theNeutron;
143  Z++;
144  } else if(theParticle == theNeutron) {
145  theSecondary = theProton;
146  Z--;
147  } else if(theParticle == thePiPlus) {
148  theSecondary = thePiZero;
149  Z++;
150  } else if(theParticle == thePiMinus) {
151  theSecondary = thePiZero;
152  Z--;
153  } else if(theParticle == theKPlus) {
154  if(G4UniformRand()<0.5) theSecondary = theK0S;
155  else theSecondary = theK0L;
156  Z++;
157  } else if(theParticle == theKMinus) {
158  if(G4UniformRand()<0.5) theSecondary = theK0S;
159  else theSecondary = theK0L;
160  Z--;
161  } else if(theParticle == theK0S || theParticle == theK0L) {
162  if(G4UniformRand()*A < G4double(Z)) {
163  theSecondary = theKPlus;
164  Z--;
165  } else {
166  theSecondary = theKMinus;
167  Z++;
168  }
169  } else if(theParticle == theANeutron) {
170  theSecondary = theAProton;
171  Z++;
172  } else if(theParticle == theAProton) {
173  theSecondary = theANeutron;
174  Z--;
175  } else if(theParticle == theL) {
176  G4double x = G4UniformRand();
177  if(G4UniformRand()*A < G4double(Z)) {
178  if(x < 0.2) {
179  theSecondary = theS0;
180  } else if (x < 0.4) {
181  theSecondary = theSPlus;
182  Z--;
183  } else if (x < 0.6) {
184  theSecondary = theProton;
185  theRecoil = theL;
186  theHyperon = true;
187  A--;
188  } else if (x < 0.8) {
189  theSecondary = theProton;
190  theRecoil = theS0;
191  theHyperon = true;
192  A--;
193  } else {
194  theSecondary = theNeutron;
195  theRecoil = theSPlus;
196  theHyperon = true;
197  A--;
198  }
199  } else {
200  if(x < 0.2) {
201  theSecondary = theS0;
202  } else if (x < 0.4) {
203  theSecondary = theSMinus;
204  Z++;
205  } else if (x < 0.6) {
206  theSecondary = theNeutron;
207  theRecoil = theL;
208  A--;
209  theHyperon = true;
210  } else if (x < 0.8) {
211  theSecondary = theNeutron;
212  theRecoil = theS0;
213  theHyperon = true;
214  A--;
215  } else {
216  theSecondary = theProton;
217  theRecoil = theSMinus;
218  theHyperon = true;
219  A--;
220  }
221  }
222  }
223 
224  if (Z == 1 && A == 2) theDef = theD;
225  else if (Z == 1 && A == 3) theDef = theT;
226  else if (Z == 2 && A == 3) theDef = theHe3;
227  else if (Z == 2 && A == 4) theDef = theA;
228  else {
229  theDef =
231  }
232  if(!theSecondary) { return &theParticleChange; }
233 
234  G4double m11 = theSecondary->GetPDGMass();
235  G4double m21 = theDef->GetPDGMass();
236  if(theRecoil) { m21 += theRecoil->GetPDGMass(); }
237  else { theRecoil = theDef; }
238 
239  G4double etot = lv0.e() + lv1.e();
240 
241  // kinematiacally impossible
242  if(etot < m11 + m21) {
244  theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
245  return &theParticleChange;
246  }
247 
248  G4ThreeVector p1 = lv1.vect();
249  G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
250  // G4double e2 = etot - e1;
251  G4double ptot = std::sqrt(e1*e1 - m11*m11);
252 
253  G4double tmax = 4.0*ptot*ptot;
254  G4double g2 = GeV*GeV;
255 
256  G4double t = g2*SampleT(tmax/g2, A);
257 
258  if(verboseLevel>1)
259  G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
260  << " ptot= " << ptot << G4endl;
261 
262  // Sampling in CM system
263  G4double phi = G4UniformRand()*twopi;
264  G4double cost = 1. - 2.0*t/tmax;
265  if(std::abs(cost) > 1.0) cost = 1.0;
266  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
267 
268  //if (verboseLevel > 1)
269  // G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
270 
271  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
272  v1 *= ptot;
273  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
274  G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
275 
276  nlv0.boost(bst);
277  nlv1.boost(bst);
278 
281  G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
283 
284  G4double erec = nlv0.e() - m21;
285 
286  //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl;
287 
288  if(theHyperon) {
290  aSec = new G4DynamicParticle();
291  aSec->SetDefinition(theRecoil);
292  aSec->SetKineticEnergy(0.0);
293  } else if(erec > lowEnergyRecoilLimit) {
294  aSec = new G4DynamicParticle(theRecoil, nlv0);
296  } else {
297  if(erec < 0.0) erec = 0.0;
299  }
300  return &theParticleChange;
301 }
302 
304 {
305  G4double aa, bb, cc, dd;
306  if (A <= 62.) {
307  aa = std::pow(A, 1.63);
308  bb = 14.5*std::pow(A, 0.66);
309  cc = 1.4*std::pow(A, 0.33);
310  dd = 10.;
311  } else {
312  aa = std::pow(A, 1.33);
313  bb = 60.*std::pow(A, 0.33);
314  cc = 0.4*std::pow(A, 0.40);
315  dd = 10.;
316  }
317  G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;
318  G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;
319 
320  G4double t;
321  G4double y = bb;
322  if(G4UniformRand()*(x1 + x2) < x2) y = dd;
323 
324  do {t = -std::log(G4UniformRand())/y;} while (t > tmax);
325 
326  return t;
327 }
328 
G4ParticleDefinition * theSMinus
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
static const double MeV
Definition: G4SIunits.hh:193
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:449
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:99
G4ParticleDefinition * theAXi0
static G4KaonMinus * KaonMinus()
Definition: G4KaonMinus.cc:113
void SetStatusChange(G4HadFinalStateStatus aS)
static G4AntiSigmaMinus * AntiSigmaMinus()
static G4XiZero * XiZero()
Definition: G4XiZero.cc:106
void SetMinEnergy(G4double anEnergy)
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
G4ParticleDefinition * theK0L
const G4ParticleDefinition * GetDefinition() const
static G4KaonZeroShort * KaonZeroShort()
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
bool G4bool
Definition: G4Types.hh:79
G4double GetKineticEnergy() const
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:196
G4ParticleDefinition * theXi0
static G4Neutron * Neutron()
Definition: G4Neutron.cc:104
G4ParticleDefinition * theKPlus
G4ParticleDefinition * theASMinus
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
static G4PionZero * PionZero()
Definition: G4PionZero.cc:104
static G4Deuteron * Deuteron()
Definition: G4Deuteron.cc:94
void SetKineticEnergy(G4double aEnergy)
static const G4double e1
G4ParticleDefinition * theANeutron
static G4SigmaMinus * SigmaMinus()
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()
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:197
static G4Alpha * Alpha()
Definition: G4Alpha.cc:89
G4ParticleDefinition * thePiZero
static const double keV
Definition: G4SIunits.hh:195
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)
void AddSecondary(G4DynamicParticle *aP)
static G4AntiNeutron * AntiNeutron()
G4double GetTotalMomentum() const
G4ParticleDefinition * theAntiL
CLHEP::HepLorentzVector G4LorentzVector
G4ParticleDefinition * theProton
G4ParticleDefinition * thePiMinus
G4ParticleDefinition * theHe3