Geant4  10.01.p02
G4TripathiLightCrossSection.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 // * *
21 // * Parts of this code which have been developed by QinetiQ Ltd *
22 // * under contract to the European Space Agency (ESA) are the *
23 // * intellectual property of ESA. Rights to use, copy, modify and *
24 // * redistribute this software for general public use are granted *
25 // * in compliance with any licensing, distribution and development *
26 // * policy adopted by the Geant4 Collaboration. This code has been *
27 // * written by QinetiQ Ltd for the European Space Agency, under ESA *
28 // * contract 17191/03/NL/LvH (Aurora Programme). *
29 // * *
30 // * By using, copying, modifying or distributing the software (or *
31 // * any work based on the software) you agree to acknowledge its *
32 // * use in resulting scientific publications, and indicate your *
33 // * acceptance of all terms of the Geant4 Software license. *
34 // ********************************************************************
35 //
36 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
37 //
38 // MODULE: G4TripathiLightCrossSection.cc
39 //
40 // Version: B.1
41 // Date: 15/04/04
42 // Author: P R Truscott
43 // Organisation: QinetiQ Ltd, UK
44 // Customer: ESA/ESTEC, NOORDWIJK
45 // Contract: 17191/03/NL/LvH
46 //
47 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
48 //
49 // CHANGE HISTORY
50 // --------------
51 //
52 // 6 October 2003, P R Truscott, QinetiQ Ltd, UK
53 // Created.
54 //
55 // 15 March 2004, P R Truscott, QinetiQ Ltd, UK
56 // Beta release
57 //
58 // 24 November 2010 J. M. Quesada bug fixed in X_m
59 // (according to eq. 14 in
60 // R.K. Tripathi et al. Nucl. Instr. and Meth. in Phys. Res. B 155 (1999) 349-356)
61 //
62 // 19 Aug 2011 V.Ivanchenko move to new design and make x-section per element
63 //
64 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
66 //
68 #include "G4PhysicalConstants.hh"
69 #include "G4SystemOfUnits.hh"
70 #include "G4DynamicParticle.hh"
71 #include "G4WilsonRadius.hh"
72 #include "G4NucleiProperties.hh"
73 #include "G4HadTmpUtil.hh"
74 #include "G4NistManager.hh"
75 #include "G4Pow.hh"
76 
78  : G4VCrossSectionDataSet("TripathiLightIons")
79 {
80  // Constructor only needs to instantiate the object which provides functions
81  // to calculate the nuclear radius, and some other constants used to
82  // calculate cross-sections.
83 
85  r_0 = 1.1 * fermi;
86 
87  // The following variable is set to true if
88  // G4TripathiLightCrossSection::GetCrossSection is going to be called from
89  // within G4TripathiLightCrossSection::GetCrossSection to check whether the
90  // cross-section is behaviing anomalously in the low-energy region.
91 
92  lowEnergyCheck = false;
93 }
95 //
97 {
98  //
99  // Destructor just needs to delete the pointer to the G4WilsonRadius object.
100  //
101  delete theWilsonRadius;
102 }
104 //
105 G4bool
107  G4int ZT, const G4Material*)
108 {
109  G4bool result = false;
110  G4int AT = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(ZT));
111  G4int ZP = G4lrint(theProjectile->GetDefinition()->GetPDGCharge()/eplus);
112  G4int AP = theProjectile->GetDefinition()->GetBaryonNumber();
113  if (theProjectile->GetKineticEnergy()/AP < 10.0*GeV &&
114  ((AT==1 && ZT==1) || (AP==1 && ZP==1) ||
115  (AT==1 && ZT==0) || (AP==1 && ZP==0) ||
116  (AT==2 && ZT==1) || (AP==2 && ZP==1) ||
117  (AT==3 && ZT==2) || (AP==3 && ZP==2) ||
118  (AT==4 && ZT==2) || (AP==4 && ZP==2))) { result = true; }
119  return result;
120 }
121 
123 //
124 G4double
126  G4int ZT, const G4Material*)
127 {
128  // Initialise the result.
129  G4double result = 0.0;
130 
131  // Get details of the projectile and target (nucleon number, atomic number,
132  // kinetic enery and energy/nucleon.
133 
135  G4int AT = G4lrint(xAT);
136  G4double EA = theProjectile->GetKineticEnergy()/MeV;
137  G4int AP = theProjectile->GetDefinition()->GetBaryonNumber();
138  G4double xAP= G4double(AP);
139  G4double ZP = G4lrint(theProjectile->GetDefinition()->GetPDGCharge()/eplus);
140  G4double E = EA / xAP;
141 
142  G4Pow* g4pow = G4Pow::GetInstance();
143 
144  G4double AT13 = g4pow->Z13(AT);
145  G4double AP13 = g4pow->Z13(AP);
146 
147  // Determine target mass and energy within the centre-of-mass frame.
148 
150  G4LorentzVector pT(0.0, 0.0, 0.0, mT);
151  G4LorentzVector pP(theProjectile->Get4Momentum());
152  pT += pP;
153  G4double E_cm = (pT.mag()-mT-pP.m())/MeV;
154 
155  //G4cout << G4endl;
156  //G4cout << "### EA= " << EA << " ZT= " << ZT << " AT= " << AT
157  // << " ZP= " << ZP << " AP= " << AP << " E_cm= " << E_cm
158  // << " Elim= " << (0.8 + 0.04*ZT)*xAP << G4endl;
159 
160  if (E_cm <= 0.0) { return 0.; }
161  if (E_cm <= (0.8 + 0.04*ZT)*xAP && !lowEnergyCheck) { return 0.; }
162 
163  G4double E_cm13 = g4pow->A13(E_cm);
164 
165  // Determine nuclear radii. Note that the r_p and r_T are defined differently
166  // from Wilson et al.
167 
170 
171  G4double r_p = 1.29*r_rms_p;
172  G4double r_t = 1.29*r_rms_t;
173 
174  G4double Radius = (r_p + r_t)/fermi + 1.2*(AT13 + AP13)/E_cm13;
175 
176  G4double B = 1.44 * ZP * ZT / Radius;
177 
178  // Now determine other parameters associated with the parametric
179  // formula, depending upon the projectile and target.
180 
181  G4double T1 = 0.0;
182  G4double D = 0.0;
183  G4double G = 0.0;
184 
185  if ((AT==1 && ZT==1) || (AP==1 && ZP==1)) {
186  T1 = 23.0;
187  D = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
188 
189  } else if ((AT==1 && ZT==0) || (AP==1 && ZP==0)) {
190  T1 = 18.0;
191  D = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
192 
193  } else if ((AT==2 && ZT==1) || (AP==2 && ZP==1)) {
194  T1 = 23.0;
195  D = 1.65 + 0.1/(1+std::exp((500.0-E)/200.0));
196 
197  } else if ((AT==3 && ZT==2) || (AP==3 && ZP==2)) {
198  T1 = 40.0;
199  D = 1.55;
200 
201  } else if (AP==4 && ZP==2) {
202  if (AT==4 && ZT==2) {T1 = 40.0; G = 300.0;}
203  else if (ZT==4) {T1 = 25.0; G = 300.0;}
204  else if (ZT==7) {T1 = 40.0; G = 500.0;}
205  else if (ZT==13) {T1 = 25.0; G = 300.0;}
206  else if (ZT==26) {T1 = 40.0; G = 300.0;}
207  else {T1 = 40.0; G = 75.0;}
208  D = 2.77 - 8.0E-3*AT + 1.8E-5*AT*AT-0.8/(1.0+std::exp((250.0-E)/G));
209  }
210  else if (AT==4 && ZT==2) {
211  if (AP==4 && ZP==2) {T1 = 40.0; G = 300.0;}
212  else if (ZP==4) {T1 = 25.0; G = 300.0;}
213  else if (ZP==7) {T1 = 40.0; G = 500.0;}
214  else if (ZP==13) {T1 = 25.0; G = 300.0;}
215  else if (ZP==26) {T1 = 40.0; G = 300.0;}
216  else {T1 = 40.0; G = 75.0;}
217  D = 2.77 - 8.0E-3*AP + 1.8E-5*AP*AP-0.8/(1.0+std::exp((250.0-E)/G));
218  }
219 
220  // C_E, S, deltaE, X1, S_L and X_m correspond directly with the original
221  // formulae of Tripathi et al in his report.
222  //G4cout << "E= " << E << " T1= " << T1 << " AP= " << AP << " ZP= " << ZP
223  // << " AT= " << AT << " ZT= " << ZT << G4endl;
224  G4double C_E = D*(1.0-std::exp(-E/T1)) -
225  0.292*std::exp(-E/792.0)*std::cos(0.229*std::pow(E,0.453));
226 
227  G4double S = AP13*AT13/(AP13 + AT13);
228 
229  G4double deltaE = 0.0;
230  G4double X1 = 0.0;
231  if (AT >= AP)
232  {
233  deltaE = 1.85*S + 0.16*S/E_cm13 - C_E + 0.91*(AT-2*ZT)*ZP/(xAT*xAP);
234  X1 = 2.83 - 3.1E-2*AT + 1.7E-4*AT*AT;
235  }
236  else
237  {
238  deltaE = 1.85*S + 0.16*S/E_cm13 - C_E + 0.91*(AP-2*ZP)*ZT/(xAT*xAP);
239  X1 = 2.83 - 3.1E-2*AP + 1.7E-4*AP*AP;
240  }
241  G4double S_L = 1.2 + 1.6*(1.0-std::exp(-E/15.0));
242  //JMQ 241110 bug fixed
243  G4double X_m = 1.0 - X1*std::exp(-E/(X1*S_L));
244 
245  //G4cout << "deltaE= " << deltaE << " X1= " << X1 << " S_L= " << S_L << " X_m= " << X_m << G4endl;
246 
247  // R_c is also highly dependent upon the A and Z of the projectile and
248  // target.
249 
250  G4double R_c = 1.0;
251  if (AP==1 && ZP==1)
252  {
253  if (AT==2 && ZT==1) R_c = 13.5;
254  else if (AT==3 && ZT==2) R_c = 21.0;
255  else if (AT==4 && ZT==2) R_c = 27.0;
256  else if (ZT==3) R_c = 2.2;
257  }
258  else if (AT==1 && ZT==1)
259  {
260  if (AP==2 && ZP==1) R_c = 13.5;
261  else if (AP==3 && ZP==2) R_c = 21.0;
262  else if (AP==4 && ZP==2) R_c = 27.0;
263  else if (ZP==3) R_c = 2.2;
264  }
265  else if (AP==2 && ZP==1)
266  {
267  if (AT==2 && ZT==1) R_c = 13.5;
268  else if (AT==4 && ZT==2) R_c = 13.5;
269  else if (AT==12 && ZT==6) R_c = 6.0;
270  }
271  else if (AT==2 && ZT==1)
272  {
273  if (AP==2 && ZP==1) R_c = 13.5;
274  else if (AP==4 && ZP==2) R_c = 13.5;
275  else if (AP==12 && ZP==6) R_c = 6.0;
276  }
277  else if ((AP==4 && ZP==2 && (ZT==73 || ZT==79)) ||
278  (AT==4 && ZT==2 && (ZP==73 || ZP==79))) R_c = 0.6;
279 
280  // Find the total cross-section. Check that it's value is positive, and if
281  // the energy is less that 10 MeV/nuc, find out if the cross-section is
282  // increasing with decreasing energy. If so this is a sign that the function
283  // is behaving badly at low energies, and the cross-section should be
284  // set to zero.
285 
286  G4double xr = r_0*(AT13 + AP13 + deltaE);
287  result = pi * xr * xr * (1.0 - R_c*B/E_cm) * X_m;
288  //G4cout << " result= " << result << " E= " << E << " check= "<< lowEnergyCheck << G4endl;
289  if (result < 0.0) {
290  result = 0.0;
291 
292  } else if (!lowEnergyCheck && E < 6.0) {
293  G4double f = 0.95;
294  G4DynamicParticle slowerProjectile = *theProjectile;
295  slowerProjectile.SetKineticEnergy(f * EA * MeV);
296 
297  G4bool savelowenergy = lowEnergyCheck;
298  SetLowEnergyCheck(true);
299  G4double resultp = GetElementCrossSection(&slowerProjectile, ZT);
300  SetLowEnergyCheck(savelowenergy);
301  //G4cout << " newres= " << resultp << " f*EA= " << f*EA << G4endl;
302  if (resultp > result) { result = 0.0; }
303  }
304 
305  return result;
306 }
307 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static const double MeV
Definition: G4SIunits.hh:193
static G4double GetNuclearMass(const G4double A, const G4double Z)
G4double GetWilsonRMSRadius(G4double A)
G4double GetKineticEnergy() const
Definition: G4Pow.hh:56
const G4double pi
G4ParticleDefinition * GetDefinition() const
int G4int
Definition: G4Types.hh:78
static G4NistManager * Instance()
virtual G4double GetElementCrossSection(const G4DynamicParticle *theProjectile, G4int Z, const G4Material *mat=0)
virtual G4bool IsElementApplicable(const G4DynamicParticle *theProjectile, G4int Z, const G4Material *)
G4double Z13(G4int Z) const
Definition: G4Pow.hh:127
bool G4bool
Definition: G4Types.hh:79
static const double GeV
Definition: G4SIunits.hh:196
void SetKineticEnergy(G4double aEnergy)
G4LorentzVector Get4Momentum() const
G4double A13(G4double A) const
Definition: G4Pow.hh:132
int G4lrint(double ad)
Definition: templates.hh:163
G4double GetAtomicMassAmu(const G4String &symb) const
double G4double
Definition: G4Types.hh:76
static const double eplus
Definition: G4SIunits.hh:178
G4double GetPDGCharge() const
static const double fermi
Definition: G4SIunits.hh:93
CLHEP::HepLorentzVector G4LorentzVector