source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4UHadronElasticProcess.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 12.9 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// $Id: G4UHadronElasticProcess.cc,v 1.39 2008/10/22 08:16:40 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29// Geant4 Hadron Elastic Scattering Process -- header file
30//
31// Created 21 April 2006 V.Ivanchenko
32//
33// Modified:
34// 24.04.06 V.Ivanchenko add neutron scattering on hydrogen from CHIPS
35// 07.06.06 V.Ivanchenko fix problem of rotation of final state
36// 25.07.06 V.Ivanchenko add 19 MeV low energy for CHIPS
37// 26.09.06 V.Ivanchenko add lowestEnergy
38// 20.10.06 V.Ivanchenko initialise lowestEnergy=0 for neitrals, eV for charged
39// 23.01.07 V.Ivanchnko add cross section interfaces with Z and A
40// 02.05.07 V.Ivanchnko add He3
41//
42
43#include "G4UHadronElasticProcess.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
57G4UHadronElasticProcess::G4UHadronElasticProcess(const G4String& pName, G4double)
58 : G4HadronicProcess(pName), lowestEnergy(0.0), first(true)
59{
60 SetProcessSubType(fHadronElastic);
61 AddDataSet(new G4HadronElasticDataSet);
62 theProton = G4Proton::Proton();
63 theNeutron = G4Neutron::Neutron();
64 thEnergy = 19.0*MeV;
65 verboseLevel= 1;
66 qCManager = 0;
67}
68
69G4UHadronElasticProcess::~G4UHadronElasticProcess()
70{
71}
72
73void G4UHadronElasticProcess::SetQElasticCrossSection(G4VQCrossSection* p)
74{
75 qCManager = p;
76}
77
78void G4UHadronElasticProcess::
79BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
80{
81 if(first) {
82 first = false;
83 if(!qCManager) qCManager = G4QElasticCrossSection::GetPointer();
84 theParticle = &aParticleType;
85 pPDG = theParticle->GetPDGEncoding();
86
87 store = G4HadronicProcess::GetCrossSectionDataStore();
88
89 // defined lowest threshold for the projectile
90 if(theParticle->GetPDGCharge() != 0.0) lowestEnergy = eV;
91
92 // if(verboseLevel>1 ||
93 // (verboseLevel==1 && theParticle == theNeutron)) {
94 if(verboseLevel>1 && theParticle == theNeutron) {
95 // G4cout << G4endl;
96 G4cout << "G4UHadronElasticProcess for "
97 << theParticle->GetParticleName()
98 << " PDGcode= " << pPDG
99 << " Elow(MeV)= " << thEnergy/MeV
100 << " Elowest(eV)= " << lowestEnergy/eV
101 << G4endl;
102 }
103 }
104 G4HadronicProcess::BuildPhysicsTable(aParticleType);
105 //store->BuildPhysicsTable(aParticleType);
106}
107
108G4double G4UHadronElasticProcess::GetMeanFreePath(const G4Track& track,
109 G4double,
110 G4ForceCondition* cond)
111{
112 *cond = NotForced;
113 const G4DynamicParticle* dp = track.GetDynamicParticle();
114 cross = 0.0;
115 G4double x = DBL_MAX;
116
117 // Compute cross sesctions
118 const G4Material* material = track.GetMaterial();
119 const G4ElementVector* theElementVector = material->GetElementVector();
120 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
121 G4double temp = material->GetTemperature();
122 G4int nelm = material->GetNumberOfElements();
123
124#ifdef G4VERBOSE
125 if(verboseLevel>1)
126 G4cout << "G4UHadronElasticProcess get mfp for "
127 << theParticle->GetParticleName()
128 << " p(GeV)= " << dp->GetTotalMomentum()/GeV
129 << " in " << material->GetName()
130 << G4endl;
131#endif
132
133 for (G4int i=0; i<nelm; i++) {
134 const G4Element* elm = (*theElementVector)[i];
135 G4double x = GetMicroscopicCrossSection(dp, elm, temp);
136 cross += theAtomNumDensityVector[i]*x;
137 xsec[i] = cross;
138 }
139
140#ifdef G4VERBOSE
141 if(verboseLevel>1)
142 G4cout << "G4UHadronElasticProcess cross(1/mm)= " << cross
143 << " E(MeV)= " << dp->GetKineticEnergy()
144 << " " << theParticle->GetParticleName()
145 << " in " << material->GetName()
146 << G4endl;
147#endif
148
149 if(cross > DBL_MIN) x = 1./cross;
150 return x;
151}
152
153G4double G4UHadronElasticProcess::GetMicroscopicCrossSection(
154 const G4DynamicParticle* dp,
155 const G4Element* elm,
156 G4double temp)
157{
158 // gives the microscopic cross section in GEANT4 internal units
159 G4int iz = G4int(elm->GetZ());
160 G4double x = 0.0;
161
162 // CHIPS cross sections
163 if(iz <= 2 && dp->GetKineticEnergy() > thEnergy &&
164 (theParticle == theProton || theParticle == theNeutron)) {
165
166 G4double momentum = dp->GetTotalMomentum();
167 G4IsotopeVector* isv = elm->GetIsotopeVector();
168 G4int ni = 0;
169 if(isv) ni = isv->size();
170
171 x = 0.0;
172 if(ni == 0) {
173 G4int N = G4int(elm->GetN()+0.5) - iz;
174 x = qCManager->GetCrossSection(false,momentum,iz,N,pPDG);
175 xsecH[0] = x;
176#ifdef G4VERBOSE
177 if(verboseLevel>1)
178 G4cout << "G4UHadronElasticProcess compute CHIPS CS for Z= " << iz
179 << " N= " << N << " pdg= " << pPDG
180 << " mom(GeV)= " << momentum/GeV
181 << " " << qCManager << G4endl;
182#endif
183 } else {
184 G4double* ab = elm->GetRelativeAbundanceVector();
185 for(G4int j=0; j<ni; j++) {
186 G4int N = (*isv)[j]->GetN() - iz;
187 if(iz == 1) {
188 if(N > 1) N = 1;
189 } else {
190 N = 2;
191 }
192#ifdef G4VERBOSE
193 if(verboseLevel>1)
194 G4cout << "G4UHadronElasticProcess compute CHIPS CS for Z= " << iz
195 << " N= " << N << " pdg= " << pPDG
196 << " mom(GeV)= " << momentum/GeV
197 << " " << qCManager << G4endl;
198#endif
199 G4double y = ab[j]*qCManager->GetCrossSection(false,momentum,iz,N,pPDG);
200 x += y;
201 xsecH[j] = x;
202 }
203 }
204
205 // GHAD cross section
206 } else {
207#ifdef G4VERBOSE
208 if(verboseLevel>1)
209 G4cout << "G4UHadronElasticProcess compute GHAD CS for element "
210 << elm->GetName()
211 << G4endl;
212#endif
213 x = store->GetCrossSection(dp, elm, temp);
214 }
215 // NaN finder
216 if(!(x < 0.0 || x >= 0.0)) {
217 if (verboseLevel > 1) {
218 G4cout << "G4UHadronElasticProcess:WARNING: Z= " << iz
219 << " pdg= " << pPDG
220 << " mom(GeV)= " << dp->GetTotalMomentum()/GeV
221 << " cross= " << x
222 << " set to zero"
223 << G4endl;
224 }
225 x = 0.0;
226 }
227
228#ifdef G4VERBOSE
229 if(verboseLevel>1)
230 G4cout << "G4UHadronElasticProcess cross(mb)= " << x/millibarn
231 << " E(MeV)= " << dp->GetKineticEnergy()
232 << " " << theParticle->GetParticleName()
233 << " in Z= " << iz
234 << G4endl;
235#endif
236
237 return x;
238}
239
240G4VParticleChange* G4UHadronElasticProcess::PostStepDoIt(
241 const G4Track& track,
242 const G4Step& step)
243{
244 G4ForceCondition cn;
245 aParticleChange.Initialize(track);
246 G4double kineticEnergy = track.GetKineticEnergy();
247 if(kineticEnergy <= lowestEnergy)
248 return G4VDiscreteProcess::PostStepDoIt(track,step);
249
250 G4double mfp = GetMeanFreePath(track, 0.0, &cn);
251 if(mfp == DBL_MAX)
252 return G4VDiscreteProcess::PostStepDoIt(track,step);
253
254 G4Material* material = track.GetMaterial();
255
256 // Select element
257 const G4ElementVector* theElementVector = material->GetElementVector();
258 G4Element* elm = (*theElementVector)[0];
259 G4int nelm = material->GetNumberOfElements() - 1;
260 if (nelm > 0) {
261 G4double x = G4UniformRand()*cross;
262 G4int i = -1;
263 do {i++;} while (x > xsec[i] && i < nelm);
264 elm = (*theElementVector)[i];
265 }
266 G4double Z = elm->GetZ();
267 G4double A = G4double(G4int(elm->GetN()+0.5));
268 G4int iz = G4int(Z);
269
270 // Select isotope
271 G4IsotopeVector* isv = elm->GetIsotopeVector();
272 G4int ni = 0;
273 if(isv) ni = isv->size();
274
275 if(ni == 1) {
276 A = G4double((*isv)[0]->GetN());
277 } else if(ni > 1) {
278
279 G4double* ab = elm->GetRelativeAbundanceVector();
280 G4int j = -1;
281 ni--;
282 // Special treatment of hydrogen and helium for CHIPS
283 if(iz <= 2 && kineticEnergy > thEnergy &&
284 (theParticle == theProton || theParticle == theNeutron)) {
285 G4double x = G4UniformRand()*xsecH[ni];
286 do {j++;} while (x > xsecH[j] && j < ni);
287
288 // GHAD cross sections
289 } else {
290 G4double y = G4UniformRand();
291 do {
292 j++;
293 y -= ab[j];
294 } while (y > 0.0 && j < ni);
295 }
296 A = G4double((*isv)[j]->GetN());
297 }
298
299 G4HadronicInteraction* hadi =
300 ChooseHadronicInteraction( kineticEnergy, material, elm);
301
302 // Initialize the hadronic projectile from the track
303 // G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
304 G4HadProjectile thePro(track);
305 if(verboseLevel>1)
306 G4cout << "G4UHadronElasticProcess::PostStepDoIt for "
307 << theParticle->GetParticleName()
308 << " Target Z= " << Z
309 << " A= " << A << G4endl;
310 targetNucleus.SetParameters(A, Z);
311
312 aParticleChange.Initialize(track);
313 G4HadFinalState* result = hadi->ApplyYourself(thePro, targetNucleus);
314 G4ThreeVector indir = track.GetMomentumDirection();
315 G4ThreeVector outdir = (result->GetMomentumChange()).rotateUz(indir);
316
317 if(verboseLevel>1)
318 G4cout << "Efin= " << result->GetEnergyChange()
319 << " de= " << result->GetLocalEnergyDeposit()
320 << " nsec= " << result->GetNumberOfSecondaries()
321 << " dir= " << outdir
322 << G4endl;
323
324 aParticleChange.ProposeEnergy(result->GetEnergyChange());
325 aParticleChange.ProposeMomentumDirection(outdir);
326 if(result->GetNumberOfSecondaries() > 0) {
327 aParticleChange.SetNumberOfSecondaries(1);
328 G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
329 G4ThreeVector pdir = p->GetMomentumDirection();
330 // G4cout << "recoil " << pdir << G4endl;
331 pdir = pdir.rotateUz(indir);
332 // G4cout << "recoil rotated " << pdir << G4endl;
333 p->SetMomentumDirection(pdir);
334 aParticleChange.AddSecondary(p);
335 } else {
336 aParticleChange.SetNumberOfSecondaries(0);
337 aParticleChange.ProposeLocalEnergyDeposit(result->GetLocalEnergyDeposit());
338 }
339 result->Clear();
340
341 return G4VDiscreteProcess::PostStepDoIt(track,step);
342}
343
344G4bool G4UHadronElasticProcess::
345IsApplicable(const G4ParticleDefinition& aParticleType)
346{
347 return (aParticleType == *(G4PionPlus::PionPlus()) ||
348 aParticleType == *(G4PionMinus::PionMinus()) ||
349 aParticleType == *(G4KaonPlus::KaonPlus()) ||
350 aParticleType == *(G4KaonZeroShort::KaonZeroShort()) ||
351 aParticleType == *(G4KaonZeroLong::KaonZeroLong()) ||
352 aParticleType == *(G4KaonMinus::KaonMinus()) ||
353 aParticleType == *(G4Proton::Proton()) ||
354 aParticleType == *(G4AntiProton::AntiProton()) ||
355 aParticleType == *(G4Neutron::Neutron()) ||
356 aParticleType == *(G4AntiNeutron::AntiNeutron()) ||
357 aParticleType == *(G4Lambda::Lambda()) ||
358 aParticleType == *(G4AntiLambda::AntiLambda()) ||
359 aParticleType == *(G4SigmaPlus::SigmaPlus()) ||
360 aParticleType == *(G4SigmaZero::SigmaZero()) ||
361 aParticleType == *(G4SigmaMinus::SigmaMinus()) ||
362 aParticleType == *(G4AntiSigmaPlus::AntiSigmaPlus()) ||
363 aParticleType == *(G4AntiSigmaZero::AntiSigmaZero()) ||
364 aParticleType == *(G4AntiSigmaMinus::AntiSigmaMinus()) ||
365 aParticleType == *(G4XiZero::XiZero()) ||
366 aParticleType == *(G4XiMinus::XiMinus()) ||
367 aParticleType == *(G4AntiXiZero::AntiXiZero()) ||
368 aParticleType == *(G4AntiXiMinus::AntiXiMinus()) ||
369 aParticleType == *(G4Deuteron::Deuteron()) ||
370 aParticleType == *(G4Triton::Triton()) ||
371 aParticleType == *(G4He3::He3()) ||
372 aParticleType == *(G4Alpha::Alpha()) ||
373 aParticleType == *(G4OmegaMinus::OmegaMinus()) ||
374 aParticleType == *(G4AntiOmegaMinus::AntiOmegaMinus()));
375}
376
377void G4UHadronElasticProcess::
378DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
379{
380 store->DumpPhysicsTable(aParticleType);
381}
382
Note: See TracBrowser for help on using the repository browser.