Geant4  10.01.p03
G4TripathiCrossSection.cc
Go to the documentation of this file.
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 // Implementation of formulas in analogy to NASA technical paper 3621 by
27 // Tripathi, et al.
28 //
29 // 26-Dec-2006 Isotope dependence added by D. Wright
30 // 19-Aug-2011 V.Ivanchenko move to new design and make x-section per element
31 //
32 
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4DynamicParticle.hh"
37 #include "G4ParticleTable.hh"
38 #include "G4IonTable.hh"
39 #include "G4HadTmpUtil.hh"
40 #include "G4Proton.hh"
41 #include "G4NistManager.hh"
42 
44  : G4VCrossSectionDataSet("Tripathi")
45 {}
46 
48 {}
49 
50 G4bool
52  G4int, const G4Material*)
53 {
54  G4bool result = false;
55  if ( (aPart->GetDefinition()->GetBaryonNumber()>2.5) &&
56  ( aPart->GetKineticEnergy()/aPart->GetDefinition()->GetBaryonNumber()<1*GeV) ) {
57  result = true;
58  }
59  return result;
60 }
61 
64  const G4Material*)
65 {
66  G4double result = 0.;
67  G4double targetAtomicNumber = G4NistManager::Instance()->GetAtomicMassAmu(ZZ);
68  G4double nTargetProtons = ZZ;
69 
70  G4double kineticEnergy = aPart->GetKineticEnergy()/MeV;
71  G4double nProjProtons = aPart->GetDefinition()->GetPDGCharge();
72  G4double projectileAtomicNumber =
73  aPart->GetDefinition()->GetBaryonNumber();
74 
75  static const G4double nuleonRadius=1.1E-15;
76  static const G4double myNuleonRadius=1.36E-15;
77 
78  // needs target mass
79  G4double targetMass =
81  ->GetIonMass(G4lrint(nTargetProtons), G4lrint(targetAtomicNumber));
82  G4LorentzVector pTarget(0,0,0,targetMass);
83  G4LorentzVector pProjectile(aPart->Get4Momentum());
84  pTarget = pTarget+pProjectile;
85  G4double E_cm = (pTarget.mag()-targetMass-pProjectile.m())/MeV;
86  if(E_cm <= DBL_MIN) { return result; }
87  // done
88  G4double r_rms_p = 0.6 * myNuleonRadius *
89  std::pow(projectileAtomicNumber, 1./3.);
90  G4double r_rms_t = 0.6 * myNuleonRadius *
91  std::pow(targetAtomicNumber, 1./3.);
92 
93  // done
94  G4double r_p = 1.29*r_rms_p/nuleonRadius ;
95  G4double r_t = 1.29*r_rms_t/nuleonRadius;
96 
97  // done
98  G4double Radius = r_p + r_t +
99  1.2*(std::pow(targetAtomicNumber, 1./3.) +
100  std::pow(projectileAtomicNumber, 1./3.))/std::pow(E_cm, 1./3.);
101 
102  //done
103  G4double B = 1.44*nProjProtons*nTargetProtons/Radius;
104  if(E_cm <= B) return result;
105  // done
106  G4double Energy = kineticEnergy/projectileAtomicNumber;
107 
108  // done
109  //
110  // Note that this correction to G4TripathiCrossSection is just to accurately
111  // reflect Tripathi's algorithm. However, if you're using alpha
112  // particles/protons consider using the more accurate
113  // G4TripathiLightCrossSection, which Tripathi developed specifically for
114  // light systems.
115  //
116 
117  G4double D;
118  if (nProjProtons==1 && projectileAtomicNumber==1)
119  {
120  D = 2.05;
121  }
122  else if (nProjProtons==2 && projectileAtomicNumber==4)
123  {
124  D = 2.77-(8.0E-3*targetAtomicNumber)+
125  (1.8E-5*targetAtomicNumber*targetAtomicNumber)
126  - 0.8/(1+std::exp((250.-Energy)/75.));
127  }
128  else
129  {
130  //
131  // This is the original value used in the G4TripathiCrossSection
132  // implementation, and was used for all projectile/target conditions.
133  // I'm not touching this, although judging from Tripathi's paper, this is
134  // valid for cases where the nucleon density changes little with A.
135  //
136  D = 1.75;
137  }
138  // done
139  G4double C_E = D * (1-std::exp(-Energy/40.)) -
140  0.292*std::exp(-Energy/792.)*std::cos(0.229*std::pow(Energy, 0.453));
141 
142  // done
143  G4double S = std::pow(projectileAtomicNumber, 1./3.)*
144  std::pow(targetAtomicNumber, 1./3.)/
145  (std::pow(projectileAtomicNumber, 1./3.) +
146  std::pow(targetAtomicNumber, 1./3.));
147 
148  // done
149  G4double deltaE = 1.85*S + 0.16*S/std::pow(E_cm,1./3.) - C_E +
150  0.91*(targetAtomicNumber-2.*nTargetProtons)*nProjProtons/
151  (targetAtomicNumber*projectileAtomicNumber);
152 
153  // done
154  result = pi * nuleonRadius*nuleonRadius *
155  std::pow(( std::pow(targetAtomicNumber, 1./3.) +
156  std::pow(projectileAtomicNumber, 1./3.) + deltaE),2.) *
157  (1-B/E_cm);
158 
159  if(result < 0.) { result = 0.; }
160  return result*m2;
161 
162 }
163 
static const double MeV
Definition: G4SIunits.hh:193
G4double GetKineticEnergy() const
const G4double pi
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
G4IonTable * GetIonTable() const
static const double m2
Definition: G4SIunits.hh:111
bool G4bool
Definition: G4Types.hh:79
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
Definition: G4IonTable.cc:1249
static const double GeV
Definition: G4SIunits.hh:196
G4LorentzVector Get4Momentum() const
static G4ParticleTable * GetParticleTable()
int G4lrint(double ad)
Definition: templates.hh:163
#define DBL_MIN
Definition: templates.hh:75
G4double GetAtomicMassAmu(const G4String &symb) const
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
virtual G4double GetElementCrossSection(const G4DynamicParticle *, G4int Z, const G4Material *)
virtual G4bool IsElementApplicable(const G4DynamicParticle *aPart, G4int Z, const G4Material *)
CLHEP::HepLorentzVector G4LorentzVector