ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4IonsKoxCrossSection.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4IonsKoxCrossSection.cc
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 // 18-Sep-2003 First version is written by T. Koi
27 // 10-Nov-2003 Bug fix at Cal. ke_per_n and D T. Koi
28 // 12-Nov-2003 Add energy check at lower side T. Koi
29 // 26-Dec-2006 Add isotope dependence D. Wright
30 // 14-Mar-2011 Moved constructor, destructor and virtual methods to source by V.Ivanchenko
31 // 19-Aug-2011 V.Ivanchenko move to new design and make x-section per element
32 
33 #include "G4IonsKoxCrossSection.hh"
34 #include "G4PhysicalConstants.hh"
35 #include "G4SystemOfUnits.hh"
36 #include "G4DynamicParticle.hh"
37 #include "G4NucleiProperties.hh"
38 #include "G4HadTmpUtil.hh"
39 #include "G4NistManager.hh"
40 
42  : G4VCrossSectionDataSet("IonsKox"),
43 // lowerLimit ( 10*MeV ),
44  r0 ( 1.1*fermi ), rc ( 1.3*fermi )
45 {}
46 
48 {}
49 
50 void
52 {
53  outFile << "G4IonsKoxCrossSection calculates the total reaction cross\n"
54  << "section for nucleus-nucleus scattering using the Kox\n"
55  << "parameterization. It is valid for projectiles and targets\n"
56  << "of all Z, at projectile energies up to 10 GeV/n. If the\n"
57  << "projectile energy is less than 10 MeV/n, a zero cross section\n"
58  << "is returned.\n";
59 }
60 
62  G4int, const G4Material*)
63 {
64  return (1 <= aDP->GetDefinition()->GetBaryonNumber());
65 }
66 
67 G4double
69  const G4DynamicParticle* aParticle, G4int ZZ, const G4Material*)
70 {
71  G4double xsection = 0.0;
72 
73  G4int Ap = aParticle->GetDefinition()->GetBaryonNumber();
74  G4int Zp = G4int(aParticle->GetDefinition()->GetPDGCharge() / eplus + 0.5);
75  G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
76 
77  // Apply energy check, if less than lower limit then 0 value is returned
78  // if ( ke_per_N < lowerLimit ) return xsection;
79 
80  G4int At = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(ZZ));
81  G4int Zt = ZZ;
82 
83  G4double one_third = 1.0 / 3.0;
84 
85  G4double cubicrAt = G4Pow::GetInstance()->powA ( G4double(At) , G4double(one_third) );
86  G4double cubicrAp = G4Pow::GetInstance()->powA ( G4double(Ap) , G4double(one_third) );
87 
88  // rc divide fermi
89  G4double Bc = Zt * Zp / ( (rc/fermi) * (cubicrAp+cubicrAt) );
90 
91  G4double targ_mass = G4NucleiProperties::GetNuclearMass(At, Zt);
92  G4double proj_mass = aParticle->GetMass();
93  G4double proj_momentum = aParticle->GetMomentum().mag();
94 
95  G4double Ecm = calEcm ( proj_mass , targ_mass , proj_momentum );
96  if( Ecm <= Bc) return xsection;
97 
98  G4double Rvol = r0 * ( cubicrAp + cubicrAt );
99 
100 // G4double ke_per_N = aParticle->GetKineticEnergy() / Ap;
101  G4double c = calCeValue ( ke_per_N / MeV );
102 
103  G4double a = 1.85;
104  G4double Rsurf = r0 * (a*cubicrAp * cubicrAt/(cubicrAp + cubicrAt) - c);
105  G4double D = 5.0 * ( At - 2 * Zt ) * Zp / ( Ap * At );
106  Rsurf = Rsurf + D * fermi; // multiply D by fermi
107 
108  G4double Rint = Rvol + Rsurf;
109  xsection = pi * Rint * Rint * ( 1 - Bc / ( Ecm / MeV ) );
110 
111  return xsection;
112 }
113 
114 G4double
116 {
117  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
118  G4double Ecm = std::sqrt ( mp * mp + mt * mt + 2 * Elab * mt );
119  G4double Pcm = Plab * mt / Ecm;
120  G4double KEcm = std::sqrt ( Pcm * Pcm + mp * mp ) - mp;
121  return KEcm;
122 }
123 
124 
126 {
127  // Calculate c value
128  // This value is indepenent from projectile and target particle
129  // ke is projectile kinetic energy per nucleon in the Lab system with MeV unit
130  // fitting function is made by T. Koi
131  // There are no data below 30 MeV/n in Kox et al.,
132 
133  G4double Ce;
134  G4double log10_ke = std::log10 ( ke );
135  if (log10_ke > 1.5)
136  {
137  Ce = - 10.0 / G4Pow::GetInstance()->powA ( G4double(log10_ke) , G4double(5) ) + 2.0;
138  }
139  else
140  {
141  Ce = (-10.0/G4Pow::GetInstance()->powA(G4double(1.5), G4double(5) ) + 2.0) /
143 
144  }
145  return Ce;
146 }
147