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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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