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

Last change on this file since 962 was 962, checked in by garnier, 15 years ago

update processes
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-ref-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.