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

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 11.4 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.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.