source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchange.cc @ 1337

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago

tag geant4.9.4 beta 1 + modifs locales
File size: 10.2 KB
Line 
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,v 1.16 2009/09/22 16:21:46 vnivanch Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30//
31// G4 Model: Charge and strangness exchange based on G4LightMedia model
32//           28 May 2006 V.Ivanchenko
33//
34// Modified:
35// 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
36// 25-Jul-06 V.Ivanchenko add 19 MeV low energy, below which S-wave is sampled
37//
38
39#include "G4ChargeExchange.hh"
40#include "G4ParticleTable.hh"
41#include "G4ParticleDefinition.hh"
42#include "G4IonTable.hh"
43#include "Randomize.hh"
44#include "G4NucleiProperties.hh"
45
46G4ChargeExchange::G4ChargeExchange() : G4HadronicInteraction("Charge Exchange")
47{
48  SetMinEnergy( 0.0*GeV );
49  SetMaxEnergy( 100.*TeV );
50
51  lowEnergyRecoilLimit = 100.*keV; 
52  lowestEnergyLimit    = 1.*MeV; 
53
54  theProton   = G4Proton::Proton();
55  theNeutron  = G4Neutron::Neutron();
56  theAProton  = G4AntiProton::AntiProton();
57  theANeutron = G4AntiNeutron::AntiNeutron();
58  thePiPlus   = G4PionPlus::PionPlus();
59  thePiMinus  = G4PionMinus::PionMinus();
60  thePiZero   = G4PionZero::PionZero();
61  theKPlus    = G4KaonPlus::KaonPlus();
62  theKMinus   = G4KaonMinus::KaonMinus();
63  theK0S      = G4KaonZeroShort::KaonZeroShort();
64  theK0L      = G4KaonZeroLong::KaonZeroLong();
65  theL        = G4Lambda::Lambda();
66  theAntiL    = G4AntiLambda::AntiLambda();
67  theSPlus    = G4SigmaPlus::SigmaPlus();
68  theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
69  theSMinus   = G4SigmaMinus::SigmaMinus();
70  theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
71  theS0       = G4SigmaZero::SigmaZero();
72  theAS0      = G4AntiSigmaZero::AntiSigmaZero();
73  theXiMinus  = G4XiMinus::XiMinus();
74  theXi0      = G4XiZero::XiZero();
75  theAXiMinus = G4AntiXiMinus::AntiXiMinus();
76  theAXi0     = G4AntiXiZero::AntiXiZero();
77  theOmega    = G4OmegaMinus::OmegaMinus();
78  theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
79  theD        = G4Deuteron::Deuteron();
80  theT        = G4Triton::Triton();
81  theA        = G4Alpha::Alpha();
82  theA        = G4He3::He3();
83}
84
85G4ChargeExchange::~G4ChargeExchange()
86{}
87
88G4HadFinalState* G4ChargeExchange::ApplyYourself(
89                 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
90{
91  theParticleChange.Clear();
92  const G4HadProjectile* aParticle = &aTrack;
93  G4double ekin = aParticle->GetKineticEnergy();
94
95  G4double aTarget = targetNucleus.GetN();
96  G4double zTarget = targetNucleus.GetZ();
97
98  G4int Z = static_cast<G4int>(zTarget+0.5);
99  G4int A = static_cast<G4int>(aTarget+0.5);
100
101  if(ekin <= lowestEnergyLimit || A < 3) {
102    theParticleChange.SetEnergyChange(ekin);
103    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
104    return &theParticleChange;
105  }
106
107  G4double plab = aParticle->GetTotalMomentum();
108
109  if (verboseLevel > 1)
110    G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
111           << plab/GeV << " GeV/c "
112           << " ekin(MeV) = " << ekin/MeV << "  "
113           << aParticle->GetDefinition()->GetParticleName() << G4endl;
114
115  // Scattered particle referred to axis of incident particle
116  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
117
118  G4int N = A - Z;
119  G4int projPDG = theParticle->GetPDGEncoding();
120  if (verboseLevel > 1)
121    G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
122           << " PDGcode= " << projPDG << " on nucleus Z= " << Z
123           << " A= " << A << " N= " << N
124           << G4endl;
125
126  G4ParticleDefinition * theDef = 0;
127
128  G4double m2 = G4NucleiProperties::GetNuclearMass((G4double)A, (G4double)Z);
129  G4LorentzVector lv1 = aParticle->Get4Momentum();
130  G4LorentzVector lv0(0.0,0.0,0.0,m2);
131
132  G4LorentzVector lv  = lv0 + lv1;
133  G4ThreeVector bst = lv.boostVector();
134  lv1.boost(-bst);
135  lv0.boost(-bst);
136
137  // Sample final particles
138  G4bool theHyperon = false;
139  G4ParticleDefinition* theRecoil = 0;
140  G4ParticleDefinition* theSecondary = 0;
141
142  if(theParticle == theProton) {
143    theSecondary = theNeutron;
144    Z++;
145  } else if(theParticle == theNeutron) {
146    theSecondary = theProton;
147    Z--;
148  } else if(theParticle == thePiPlus) {
149    theSecondary = thePiZero;
150    Z++;
151  } else if(theParticle == thePiMinus) {
152    theSecondary = thePiZero;
153    Z--;
154  } else if(theParticle == theKPlus) {
155    if(G4UniformRand()<0.5) theSecondary = theK0S;
156    else  theSecondary = theK0L;
157    Z++;
158  } else if(theParticle == theKMinus) {
159    if(G4UniformRand()<0.5) theSecondary = theK0S;
160    else  theSecondary = theK0L;
161    Z--;
162  } else if(theParticle == theK0S || theParticle == theK0L) {
163    if(G4UniformRand()*aTarget < zTarget) {
164      theSecondary = theKPlus;
165      Z--;
166    } else {
167      theSecondary = theKMinus;
168      Z++;
169    }
170  } else if(theParticle == theANeutron) {
171    theSecondary = theAProton;
172    Z++;
173  } else if(theParticle == theAProton) {
174    theSecondary = theANeutron;
175    Z--;
176  } else if(theParticle == theL) {
177    G4double x = G4UniformRand();
178    if(G4UniformRand()*aTarget < zTarget) {
179      if(x < 0.2) {
180        theSecondary = theS0;
181      } else if (x < 0.4) {
182        theSecondary = theSPlus;
183        Z--;
184      } else if (x < 0.6) {
185        theSecondary = theProton;
186        theRecoil = theL;
187        theHyperon = true;
188        A--;
189      } else if (x < 0.8) {
190        theSecondary = theProton;
191        theRecoil = theS0;
192        theHyperon = true;
193        A--;
194      } else {
195        theSecondary = theNeutron;
196        theRecoil = theSPlus;
197        theHyperon = true;
198        A--;
199      }
200    } else {
201      if(x < 0.2) {
202        theSecondary = theS0;
203      } else if (x < 0.4) {
204        theSecondary = theSMinus;
205        Z++;
206      } else if (x < 0.6) {
207        theSecondary = theNeutron;
208        theRecoil = theL;
209        A--;
210        theHyperon = true;
211      } else if (x < 0.8) {
212        theSecondary = theNeutron;
213        theRecoil = theS0;
214        theHyperon = true;
215        A--;
216      } else {
217        theSecondary = theProton;
218        theRecoil = theSMinus;
219        theHyperon = true;
220        A--;
221      }
222    }
223  }
224
225  if (Z == 1 && A == 2) theDef = theD;
226  else if (Z == 1 && A == 3) theDef = theT;
227  else if (Z == 2 && A == 3) theDef = theHe3;
228  else if (Z == 2 && A == 4) theDef = theA;
229  else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
230
231  G4double m11 = theSecondary->GetPDGMass();
232  G4double m21 = theDef->GetPDGMass();
233  if(theRecoil)  m21 += theRecoil->GetPDGMass();
234  else           theRecoil = theDef;
235
236  G4double etot = lv0.e() + lv1.e();
237
238  // kinematiacally impossible
239  if(etot < m11 + m21) {
240    theParticleChange.SetEnergyChange(ekin);
241    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
242    return &theParticleChange;
243  }
244
245  G4ThreeVector p1 = lv1.vect();
246  G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
247  //  G4double e2 = etot - e1;
248  G4double ptot = std::sqrt(e1*e1 - m11*m11);
249
250  G4double tmax = 4.0*ptot*ptot;
251  G4double g2 = GeV*GeV; 
252
253  G4double t = g2*SampleT(tmax/g2,aTarget);
254
255  if(verboseLevel>1)
256    G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
257           << " ptot= " << ptot << G4endl;
258
259  // Sampling in CM system
260  G4double phi  = G4UniformRand()*twopi;
261  G4double cost = 1. - 2.0*t/tmax;
262  if(std::abs(cost) > 1.0) cost = 1.0;
263  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
264
265  //if (verboseLevel > 1)
266  //  G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
267
268  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
269  v1 *= ptot;
270  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
271  G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
272
273  nlv0.boost(bst);
274  nlv1.boost(bst);
275
276  theParticleChange.SetStatusChange(stopAndKill);
277  theParticleChange.SetEnergyChange(0.0);
278  G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
279  theParticleChange.AddSecondary(aSec);
280
281  G4double erec =  nlv0.e() - m21;
282
283  //G4cout << "erec= " <<erec << " Esec= " << aSec->GetKineticEnergy() << G4endl; 
284
285  if(theHyperon) {
286    theParticleChange.SetLocalEnergyDeposit(erec);
287    aSec = new G4DynamicParticle();
288    aSec->SetDefinition(theRecoil);
289    aSec->SetKineticEnergy(0.0);
290  } else if(erec > lowEnergyRecoilLimit) {
291    aSec = new G4DynamicParticle(theRecoil, nlv0);
292    theParticleChange.AddSecondary(aSec);
293  } else {
294    if(erec < 0.0) erec = 0.0;
295    theParticleChange.SetLocalEnergyDeposit(erec);
296  }
297  return &theParticleChange;
298}
299
300G4double G4ChargeExchange::SampleT(G4double tmax, G4double A)
301{
302  G4double aa, bb, cc, dd;
303  if (A <= 62.) {
304    aa = std::pow(A, 1.63);
305    bb = 14.5*std::pow(A, 0.66);
306    cc = 1.4*std::pow(A, 0.33);
307    dd = 10.;
308  } else {
309    aa = std::pow(A, 1.33);
310    bb = 60.*std::pow(A, 0.33);
311    cc = 0.4*std::pow(A, 0.40);
312    dd = 10.;
313  }
314  G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;
315  G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;
316 
317  G4double t;
318  G4double y = bb;
319  if(G4UniformRand()*(x1 + x2) < x2) y = dd;
320
321  do {t = -std::log(G4UniformRand())/y;} while (t > tmax);
322
323  return t;
324}
325
Note: See TracBrowser for help on using the repository browser.