source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4ChargeExchangeProcess.cc@ 900

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

import all except CVS
File size: 11.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: G4ChargeExchangeProcess.cc,v 1.9 2007/01/30 10:23:26 vnivanch Exp $
28// GEANT4 tag $Name: $
29//
30//
31// Geant4 Hadron Elastic Scattering Process -- header file
32//
33// Created 21 April 2006 V.Ivanchenko
34//
35// Modified:
36// 24-Apr-06 V.Ivanchenko add neutron scattering on hydrogen from CHIPS
37// 07-Jun-06 V.Ivanchenko fix problem of rotation of final state
38// 25-Jul-06 V.Ivanchenko add 19 MeV low energy for CHIPS
39// 23-Jan-07 V.Ivanchenko add cross section interfaces with Z and A
40// and do not use CHIPS for cross sections
41//
42
43#include "G4ChargeExchangeProcess.hh"
44#include "globals.hh"
45#include "G4CrossSectionDataStore.hh"
46#include "G4HadronElasticDataSet.hh"
47#include "G4VQCrossSection.hh"
48#include "G4QElasticCrossSection.hh"
49#include "G4QCHIPSWorld.hh"
50#include "G4Element.hh"
51#include "G4ElementVector.hh"
52#include "G4IsotopeVector.hh"
53#include "G4Neutron.hh"
54#include "G4Proton.hh"
55#include "G4HadronElastic.hh"
56#include "G4PhysicsLinearVector.hh"
57
58G4ChargeExchangeProcess::G4ChargeExchangeProcess(const G4String& procName)
59 : G4HadronicProcess(procName), first(true)
60{
61 thEnergy = 19.*MeV;
62 verboseLevel= 1;
63 qCManager = 0;
64 AddDataSet(new G4HadronElasticDataSet);
65 theProton = G4Proton::Proton();
66 theNeutron = G4Neutron::Neutron();
67 theAProton = G4AntiProton::AntiProton();
68 theANeutron = G4AntiNeutron::AntiNeutron();
69 thePiPlus = G4PionPlus::PionPlus();
70 thePiMinus = G4PionMinus::PionMinus();
71 thePiZero = G4PionZero::PionZero();
72 theKPlus = G4KaonPlus::KaonPlus();
73 theKMinus = G4KaonMinus::KaonMinus();
74 theK0S = G4KaonZeroShort::KaonZeroShort();
75 theK0L = G4KaonZeroLong::KaonZeroLong();
76 theL = G4Lambda::Lambda();
77 theAntiL = G4AntiLambda::AntiLambda();
78 theSPlus = G4SigmaPlus::SigmaPlus();
79 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
80 theSMinus = G4SigmaMinus::SigmaMinus();
81 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
82 theS0 = G4SigmaZero::SigmaZero();
83 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
84 theXiMinus = G4XiMinus::XiMinus();
85 theXi0 = G4XiZero::XiZero();
86 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
87 theAXi0 = G4AntiXiZero::AntiXiZero();
88 theOmega = G4OmegaMinus::OmegaMinus();
89 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
90 theD = G4Deuteron::Deuteron();
91 theT = G4Triton::Triton();
92 theA = G4Alpha::Alpha();
93 theA = G4He3::He3();
94}
95
96G4ChargeExchangeProcess::~G4ChargeExchangeProcess()
97{
98 delete factors;
99}
100
101void G4ChargeExchangeProcess::SetQElasticCrossSection(G4VQCrossSection* p)
102{
103 qCManager = p;
104}
105
106void G4ChargeExchangeProcess::
107BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
108{
109 if(first) {
110 first = false;
111 theParticle = &aParticleType;
112 pPDG = theParticle->GetPDGEncoding();
113
114 store = G4HadronicProcess::GetCrossSectionDataStore();
115
116 const size_t n = 10;
117 if(theParticle == thePiPlus || theParticle == thePiMinus ||
118 theParticle == theKPlus || theParticle == theKMinus ||
119 theParticle == theK0S || theParticle == theK0L) {
120
121 G4double F[n] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.1,0.09,0.07};
122 factors = new G4PhysicsLinearVector(0.0,1.8*GeV,n);
123 for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
124
125 } else {
126
127 G4double F[n] = {0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0.0};
128 factors = new G4PhysicsLinearVector(0.0,3.6*GeV,n);
129 for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
130 }
131
132 if(verboseLevel>1)
133 G4cout << "G4ChargeExchangeProcess for "
134 << theParticle->GetParticleName()
135 << G4endl;
136 }
137 store->BuildPhysicsTable(aParticleType);
138}
139
140G4double G4ChargeExchangeProcess::GetMeanFreePath(const G4Track& track,
141 G4double,
142 G4ForceCondition* cond)
143{
144 *cond = NotForced;
145 const G4DynamicParticle* dp = track.GetDynamicParticle();
146 const G4Material* material = track.GetMaterial();
147 cross = 0.0;
148 G4double x = DBL_MAX;
149
150 // The process is effective only above the threshold
151 if(dp->GetKineticEnergy() < thEnergy) return x;
152
153 // Compute cross sesctions
154 const G4ElementVector* theElementVector = material->GetElementVector();
155 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
156 G4double temp = material->GetTemperature();
157 G4int nelm = material->GetNumberOfElements();
158 if(verboseLevel>1)
159 G4cout << "G4ChargeExchangeProcess get mfp for "
160 << theParticle->GetParticleName()
161 << " p(GeV)= " << dp->GetTotalMomentum()/GeV
162 << " in " << material->GetName()
163 << G4endl;
164 for (G4int i=0; i<nelm; i++) {
165 const G4Element* elm = (*theElementVector)[i];
166 G4double x = GetMicroscopicCrossSection(dp, elm, temp);
167 cross += theAtomNumDensityVector[i]*x;
168 xsec[i] = cross;
169 }
170 if(verboseLevel>1)
171 G4cout << "G4ChargeExchangeProcess cross(1/mm)= " << cross
172 << " E(MeV)= " << dp->GetKineticEnergy()
173 << " " << theParticle->GetParticleName()
174 << " in " << material->GetName()
175 << G4endl;
176 if(cross > DBL_MIN) x = 1./cross;
177
178 return x;
179}
180
181G4double G4ChargeExchangeProcess::GetMicroscopicCrossSection(
182 const G4DynamicParticle* dp,
183 const G4Element* elm,
184 G4double temp)
185{
186 // gives the microscopic cross section in GEANT4 internal units
187 G4double Z = elm->GetZ();
188 G4int iz = G4int(Z);
189 G4double x = 0.0;
190 if(iz == 1) return x;
191
192 if(verboseLevel>1)
193 G4cout << "G4ChargeExchangeProcess compute GHAD CS for element "
194 << elm->GetName()
195 << G4endl;
196 x = store->GetCrossSection(dp, elm, temp);
197
198 // NaN finder
199 if(!(x < 0.0 || x >= 0.0)) {
200 if (verboseLevel > -1) {
201 G4cout << "G4ChargeExchangeProcess WARNING: Z= " << iz
202 << " pdg= " << pPDG
203 << " mom(GeV)= " << dp->GetTotalMomentum()/GeV
204 << " cross= " << x
205 << " set to zero"
206 << G4endl;
207 }
208 x = 0.0;
209 }
210
211 if(verboseLevel>1)
212 G4cout << "G4ChargeExchangeProcess cross(mb)= " << x/millibarn
213 << " E(MeV)= " << dp->GetKineticEnergy()
214 << " " << theParticle->GetParticleName()
215 << " in Z= " << iz
216 << G4endl;
217 G4bool b;
218 G4double A = elm->GetN();
219 x *= factors->GetValue(dp->GetTotalMomentum(), b)/std::pow(A, 0.42);
220 if(theParticle == thePiPlus || theParticle == theProton ||
221 theParticle == theKPlus || theParticle == theANeutron)
222 x *= (1.0 - Z/A);
223
224 else if(theParticle == thePiMinus || theParticle == theNeutron ||
225 theParticle == theKMinus || theParticle == theAProton)
226 x *= Z/A;
227
228 return x;
229}
230
231G4VParticleChange* G4ChargeExchangeProcess::PostStepDoIt(
232 const G4Track& track,
233 const G4Step& step)
234{
235 G4ForceCondition* cn = 0;
236 aParticleChange.Initialize(track);
237 G4double mfp = GetMeanFreePath(track, 0.0, cn);
238 if(mfp == DBL_MAX) return G4VDiscreteProcess::PostStepDoIt(track,step);
239
240 G4double kineticEnergy = track.GetKineticEnergy();
241 G4Material* material = track.GetMaterial();
242
243 // Select element
244 const G4ElementVector* theElementVector = material->GetElementVector();
245 G4Element* elm = (*theElementVector)[0];
246 G4int nelm = material->GetNumberOfElements() - 1;
247 if (nelm > 0) {
248 G4double x = G4UniformRand()*cross;
249 G4int i = -1;
250 do {i++;} while (x > xsec[i] && i < nelm);
251 elm = (*theElementVector)[i];
252 }
253 G4double Z = elm->GetZ();
254 G4double A = G4double(G4int(elm->GetN()+0.5));
255
256 // Select isotope
257 G4IsotopeVector* isv = elm->GetIsotopeVector();
258 G4int ni = 0;
259 if(isv) ni = isv->size();
260
261 if(ni == 1) {
262 A = G4double((*isv)[0]->GetN());
263 } else if(ni > 1) {
264
265 G4double* ab = elm->GetRelativeAbundanceVector();
266 G4int j = -1;
267 ni--;
268 G4double y = G4UniformRand();
269 do {
270 j++;
271 y -= ab[j];
272 } while (y > 0.0 && j < ni);
273 A = G4double((*isv)[j]->GetN());
274 }
275 G4HadronicInteraction* hadi =
276 ChooseHadronicInteraction( kineticEnergy, material, elm);
277
278 // Initialize the hadronic projectile from the track
279 G4HadProjectile thePro(track);
280 if(verboseLevel>1)
281 G4cout << "G4ChargeExchangeProcess::PostStepDoIt for "
282 << theParticle->GetParticleName()
283 << " Target Z= " << Z
284 << " A= " << A << G4endl;
285 targetNucleus.SetParameters(A, Z);
286
287 aParticleChange.Initialize(track);
288 G4HadFinalState* result = hadi->ApplyYourself(thePro, targetNucleus);
289 G4ThreeVector indir = track.GetMomentumDirection();
290 G4int nsec = result->GetNumberOfSecondaries();
291
292 if(verboseLevel>1)
293 G4cout << "Efin= " << result->GetEnergyChange()
294 << " de= " << result->GetLocalEnergyDeposit()
295 << " nsec= " << nsec
296 << G4endl;
297
298
299 if(nsec > 0) {
300 aParticleChange.ProposeEnergy(0.0);
301 aParticleChange.ProposeTrackStatus(fStopAndKill);
302 aParticleChange.ProposeLocalEnergyDeposit(result->GetLocalEnergyDeposit());
303 aParticleChange.SetNumberOfSecondaries(nsec);
304 for(G4int j=0; j<nsec; j++) {
305 G4DynamicParticle* p = result->GetSecondary(j)->GetParticle();
306 G4ThreeVector pdir = p->GetMomentumDirection();
307 // G4cout << "recoil " << pdir << G4endl;
308 pdir = pdir.rotateUz(indir);
309 // G4cout << "recoil rotated " << pdir << G4endl;
310 p->SetMomentumDirection(pdir);
311 aParticleChange.AddSecondary(p);
312 }
313 }
314 result->Clear();
315
316 return G4VDiscreteProcess::PostStepDoIt(track,step);
317}
318
319G4bool G4ChargeExchangeProcess::
320IsApplicable(const G4ParticleDefinition& aParticleType)
321{
322 const G4ParticleDefinition* p = &aParticleType;
323 return (p == thePiPlus || p == thePiMinus ||
324 p == theProton || p == theNeutron ||
325 p == theAProton|| p == theANeutron||
326 p == theKPlus || p == theKMinus ||
327 p == theK0S || p == theK0L ||
328 p == theL);
329}
330
331void G4ChargeExchangeProcess::
332DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
333{
334 store->DumpPhysicsTable(aParticleType);
335}
Note: See TracBrowser for help on using the repository browser.