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

Last change on this file since 847 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 11.4 KB
RevLine 
[819]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.11 2007/05/25 17:46:52 dennis Exp $
28// GEANT4 tag $Name:  $
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 "G4QElasticCrossSection.hh"
44#include "G4VQCrossSection.hh"
45#include "G4ElasticHadrNucleusHE.hh"
46#include "Randomize.hh"
47#include "G4HadronElastic.hh"
48
49
50G4ChargeExchange::G4ChargeExchange(G4HadronElastic* hel, G4double elim,
51                                   G4double ehigh)
52: G4HadronicInteraction("G4ChargeExchange"),
53  fElastic(hel),
54  native(false),
55  ekinlim(elim),
56  ekinhigh(ehigh)
57{
58  SetMinEnergy( 0.0*GeV );
59  SetMaxEnergy( 100.*TeV );
60  ekinlow = 19.0*MeV;
61
62  verboseLevel= 0;
63  if(!fElastic) {
64    native = true;
65    fElastic = new G4HadronElastic();
66  }
67  qCManager   = fElastic->GetCS();
68  hElastic    = fElastic->GetHElastic();
69
70  theProton   = G4Proton::Proton();
71  theNeutron  = G4Neutron::Neutron();
72  theAProton  = G4AntiProton::AntiProton();
73  theANeutron = G4AntiNeutron::AntiNeutron();
74  thePiPlus   = G4PionPlus::PionPlus();
75  thePiMinus  = G4PionMinus::PionMinus();
76  thePiZero   = G4PionZero::PionZero();
77  theKPlus    = G4KaonPlus::KaonPlus();
78  theKMinus   = G4KaonMinus::KaonMinus();
79  theK0S      = G4KaonZeroShort::KaonZeroShort();
80  theK0L      = G4KaonZeroLong::KaonZeroLong();
81  theL        = G4Lambda::Lambda();
82  theAntiL    = G4AntiLambda::AntiLambda();
83  theSPlus    = G4SigmaPlus::SigmaPlus();
84  theASPlus   = G4AntiSigmaPlus::AntiSigmaPlus();
85  theSMinus   = G4SigmaMinus::SigmaMinus();
86  theASMinus  = G4AntiSigmaMinus::AntiSigmaMinus();
87  theS0       = G4SigmaZero::SigmaZero();
88  theAS0      = G4AntiSigmaZero::AntiSigmaZero();
89  theXiMinus  = G4XiMinus::XiMinus();
90  theXi0      = G4XiZero::XiZero();
91  theAXiMinus = G4AntiXiMinus::AntiXiMinus();
92  theAXi0     = G4AntiXiZero::AntiXiZero();
93  theOmega    = G4OmegaMinus::OmegaMinus();
94  theAOmega   = G4AntiOmegaMinus::AntiOmegaMinus();
95  theD        = G4Deuteron::Deuteron();
96  theT        = G4Triton::Triton();
97  theA        = G4Alpha::Alpha();
98  theA        = G4He3::He3();
99}
100
101G4ChargeExchange::~G4ChargeExchange()
102{
103  if(native) delete fElastic;
104}
105
106G4HadFinalState* G4ChargeExchange::ApplyYourself(
107                 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
108{
109  theParticleChange.Clear();
110  const G4HadProjectile* aParticle = &aTrack;
111  G4double ekin = aParticle->GetKineticEnergy();
112
113  G4double aTarget = targetNucleus.GetN();
114  G4double zTarget = targetNucleus.GetZ();
115
116  G4int Z = static_cast<G4int>(zTarget+0.5);
117  G4int A = static_cast<G4int>(aTarget+0.5);
118
119  if(ekin == 0.0 || A < 3) {
120    theParticleChange.SetEnergyChange(ekin);
121    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
122    return &theParticleChange;
123  }
124
125  G4double plab = aParticle->GetTotalMomentum();
126
127  if (verboseLevel > 1)
128    G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
129           << plab/GeV << " GeV/c "
130           << " ekin(MeV) = " << ekin/MeV << "  "
131           << aParticle->GetDefinition()->GetParticleName() << G4endl;
132
133  // Scattered particle referred to axis of incident particle
134  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
135  G4double m1 = theParticle->GetPDGMass();
136
137  G4int N = A - Z;
138  G4int projPDG = theParticle->GetPDGEncoding();
139  if (verboseLevel > 1)
140    G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
141           << " PDGcode= " << projPDG << " on nucleus Z= " << Z
142           << " A= " << A << " N= " << N
143           << G4endl;
144
145  G4ParticleDefinition * theDef = 0;
146
147  if (Z == 1 && A == 3) theDef = theT;
148  else if (Z == 2 && A == 3) theDef = theHe3;
149  else if (Z == 2 && A == 4) theDef = theA;
150  else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
151
152  G4double m2 = theDef->GetPDGMass();
153  G4LorentzVector lv1 = aParticle->Get4Momentum();
154  G4LorentzVector lv0(0.0,0.0,0.0,m2);
155
156  G4LorentzVector lv  = lv0 + lv1;
157  G4ThreeVector bst = lv.boostVector();
158  lv1.boost(-bst);
159  lv0.boost(-bst);
160
161  // Sample final particles
162  G4bool theHyperon = false;
163  G4ParticleDefinition* theRecoil = 0;
164  G4ParticleDefinition* theSecondary = 0;
165
166  if(theParticle == theProton) {
167    theSecondary = theNeutron;
168    Z++;
169  } else if(theParticle == theNeutron) {
170    theSecondary = theProton;
171    Z--;
172  } else if(theParticle == thePiPlus) {
173    theSecondary = thePiZero;
174    Z++;
175  } else if(theParticle == thePiMinus) {
176    theSecondary = thePiZero;
177    Z--;
178  } else if(theParticle == theKPlus) {
179    if(G4UniformRand()<0.5) theSecondary = theK0S;
180    else  theSecondary = theK0L;
181    Z++;
182  } else if(theParticle == theKMinus) {
183    if(G4UniformRand()<0.5) theSecondary = theK0S;
184    else  theSecondary = theK0L;
185    Z--;
186  } else if(theParticle == theK0S || theParticle == theK0L) {
187    if(G4UniformRand()*aTarget < zTarget) {
188      theSecondary = theKPlus;
189      Z--;
190    } else {
191      theSecondary = theKMinus;
192      Z++;
193    }
194  } else if(theParticle == theANeutron) {
195    theSecondary = theAProton;
196    Z++;
197  } else if(theParticle == theAProton) {
198    theSecondary = theANeutron;
199    Z--;
200  } else if(theParticle == theL) {
201    G4double x = G4UniformRand();
202    if(G4UniformRand()*aTarget < zTarget) {
203      if(x < 0.2) {
204        theSecondary = theS0;
205      } else if (x < 0.4) {
206        theSecondary = theSPlus;
207        Z--;
208      } else if (x < 0.6) {
209        theSecondary = theProton;
210        theRecoil = theL;
211        theHyperon = true;
212        A--;
213      } else if (x < 0.8) {
214        theSecondary = theProton;
215        theRecoil = theS0;
216        theHyperon = true;
217        A--;
218      } else {
219        theSecondary = theNeutron;
220        theRecoil = theSPlus;
221        theHyperon = true;
222        A--;
223      }
224    } else {
225      if(x < 0.2) {
226        theSecondary = theS0;
227      } else if (x < 0.4) {
228        theSecondary = theSMinus;
229        Z++;
230      } else if (x < 0.6) {
231        theSecondary = theNeutron;
232        theRecoil = theL;
233        A--;
234        theHyperon = true;
235      } else if (x < 0.8) {
236        theSecondary = theNeutron;
237        theRecoil = theS0;
238        theHyperon = true;
239        A--;
240      } else {
241        theSecondary = theProton;
242        theRecoil = theSMinus;
243        theHyperon = true;
244        A--;
245      }
246    }
247  }
248
249  if (Z == 1 && A == 2) theDef = theD;
250  else if (Z == 1 && A == 3) theDef = theT;
251  else if (Z == 2 && A == 3) theDef = theHe3;
252  else if (Z == 2 && A == 4) theDef = theA;
253  else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
254
255  G4double m11 = theSecondary->GetPDGMass();
256  G4double m21 = theDef->GetPDGMass();
257  if(theRecoil)  m21 += theRecoil->GetPDGMass();
258  else           theRecoil = theDef;
259
260  G4double etot = lv0.e() + lv1.e();
261  if(etot < m11 + m21) return &theParticleChange;
262
263  G4ThreeVector p1 = lv1.vect();
264  G4double e1 = 0.5*etot*(1.0 + (m21*m21 - m11*m11)/(etot*etot));
265  //  G4double e2 = etot - e1;
266  G4double ptot = std::sqrt(e1*e1 - m11*m11);
267
268  G4double tmax = 4.0*ptot*ptot;
269  G4double t = 0.0;
270
271  // Choose generator
272  G4ElasticGenerator gtype = fLElastic;
273  if ((theParticle == theProton || theParticle == theNeutron) && 
274      Z <= 2 && ekin >= ekinlow) {
275    gtype = fQElastic;
276  } else {
277    if(ekin >= ekinlow)       gtype = fSWave;
278    else if(ekin >= ekinhigh) gtype = fHElastic;
279  }
280
281  // Sample t
282  if(gtype == fQElastic) {
283    if (verboseLevel > 1)
284      G4cout << "G4ChargeExchange: Z= " << Z << " N= "
285             << N << " pdg= " <<  projPDG
286             << " mom(GeV)= " << plab/GeV << "  " << qCManager << G4endl;
287    if(Z == 1 && N == 2) N = 1;
288    else if (Z == 2 && N == 1) N = 2;
289    G4double cs = qCManager->GetCrossSection(false,plab,Z,N,projPDG);
290    if(cs > 0.0) t = qCManager->GetExchangeT(Z,N,projPDG);
291    else gtype = fSWave;
292  }
293
294  if(gtype == fHElastic) {
295    t = hElastic->SampleT(theParticle,plab,Z,A);
296    if(t > tmax) gtype = fSWave;
297  }
298
299  if(gtype == fLElastic) {
300    t = GeV*GeV*fElastic->SampleT(ptot,m1,m2,aTarget);
301    if(t > tmax) gtype = fSWave;
302  }
303
304  // NaN finder
305  if(!(t < 0.0 || t >= 0.0)) {
306    if (verboseLevel > -1) {
307      G4cout << "G4ChargeExchange:WARNING: Z= " << Z << " N= " 
308             << N << " pdg= " <<  projPDG
309             << " mom(GeV)= " << plab/GeV
310             << " the model type " << gtype;
311      if(gtype ==  fQElastic) G4cout << " CHIPS ";
312      else if(gtype ==  fLElastic) G4cout << " LElastic ";
313      else if(gtype ==  fHElastic) G4cout << " HElastic ";
314      G4cout << " t= " << t
315             << " S-wave will be sampled" 
316             << G4endl; 
317    }
318    gtype = fSWave;
319  }
320
321  if(gtype == fSWave) t = G4UniformRand()*tmax;
322
323  if(verboseLevel>1)
324    G4cout <<"type= " << gtype <<" t= " << t << " tmax= " << tmax
325           << " ptot= " << ptot << G4endl;
326
327  // Sampling in CM system
328  G4double phi  = G4UniformRand()*twopi;
329  G4double cost = 1. - 2.0*t/tmax;
330  if(std::abs(cost) > 1.0) cost = -1.0 + 2.0*G4UniformRand();
331  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
332
333  if (verboseLevel > 1)
334    G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
335
336  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
337  v1 *= ptot;
338  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
339  G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
340
341  nlv0.boost(bst);
342  nlv1.boost(bst);
343
344  theParticleChange.SetStatusChange(stopAndKill);
345  G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
346  theParticleChange.AddSecondary(aSec);
347
348  G4double erec =  nlv0.e() - m21;
349  if(theHyperon) {
350    theParticleChange.SetLocalEnergyDeposit(erec);
351    aSec = new G4DynamicParticle();
352    aSec->SetDefinition(theRecoil);
353    aSec->SetKineticEnergy(0.0);
354  } else if(erec > ekinlim) {
355    aSec = new G4DynamicParticle(theRecoil, nlv0);
356    theParticleChange.AddSecondary(aSec);
357  } else {
358    theParticleChange.SetLocalEnergyDeposit(erec);
359  }
360  return &theParticleChange;
361}
Note: See TracBrowser for help on using the repository browser.