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

Last change on this file since 1200 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
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-03-cand-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.