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

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

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