source: trunk/source/processes/hadronic/cross_sections/src/G4TripathiLightCrossSection.cc@ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 12.3 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// * *
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// J. M. Quesada 24 November 2010 bug fixed in X_m
59//(according to eq. 14 in R.K. Tripathi et al. Nucl. Instr. and Meth. in Phys. Res. B 155 (1999) 349-356)
60//
61// %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62///////////////////////////////////////////////////////////////////////////////
63//
64#include "G4TripathiLightCrossSection.hh"
65#include "G4WilsonRadius.hh"
66#include "G4ParticleTable.hh"
67#include "G4IonTable.hh"
68#include "G4HadTmpUtil.hh"
69
70
71///////////////////////////////////////////////////////////////////////////////
72//
73G4TripathiLightCrossSection::G4TripathiLightCrossSection ()
74{
75//
76//
77// Constructor only needs to instantiate the object which provides functions
78// to calculate the nuclear radius, and some other constants used to
79// calculate cross-sections.
80//
81 theWilsonRadius = new G4WilsonRadius();
82 r_0 = 1.1 * fermi;
83 third = 1.0/3.0;
84//
85//
86// The following variable is set to true if
87// G4TripathiLightCrossSection::GetCrossSection is going to be called from
88// within G4TripathiLightCrossSection::GetCrossSection to check whether the
89// cross-section is behaviing anomalously in the low-energy region.
90//
91 lowEnergyCheck = false;
92}
93///////////////////////////////////////////////////////////////////////////////
94//
95G4TripathiLightCrossSection::~G4TripathiLightCrossSection ()
96{
97//
98//
99// Destructor just needs to delete the pointer to the G4WilsonRadius object.
100//
101 delete theWilsonRadius;
102}
103///////////////////////////////////////////////////////////////////////////////
104//
105G4bool G4TripathiLightCrossSection::IsApplicable
106 (const G4DynamicParticle* theProjectile, const G4Element* theTarget)
107{
108 G4int Z = G4lrint(theTarget->GetZ());
109 G4int A = G4lrint(theTarget->GetN());
110 return IsIsoApplicable(theProjectile, Z, A);
111}
112
113
114G4bool
115G4TripathiLightCrossSection::IsIsoApplicable(const G4DynamicParticle* theProjectile,
116 G4int ZZ, G4int AA)
117{
118 G4bool result = false;
119 const G4double AT = AA;
120 const G4double ZT = ZZ;
121 const G4double ZP = theProjectile->GetDefinition()->GetPDGCharge();
122 const G4double AP = theProjectile->GetDefinition()->GetBaryonNumber();
123 if (theProjectile->GetKineticEnergy()/
124 theProjectile->GetDefinition()->GetBaryonNumber()<10.0*GeV &&
125 ((AT==1 && ZT==1) || (AP==1 && ZP==1) ||
126 (AT==1 && ZT==0) || (AP==1 && ZP==0) ||
127 (AT==2 && ZT==1) || (AP==2 && ZP==1) ||
128 (AT==3 && ZT==2) || (AP==3 && ZP==2) ||
129 (AT==4 && ZT==2) || (AP==4 && ZP==2))) result = true;
130 return result;
131}
132
133///////////////////////////////////////////////////////////////////////////////
134//
135G4double
136G4TripathiLightCrossSection::GetZandACrossSection(const G4DynamicParticle* theProjectile,
137 G4int ZZ, G4int AA, G4double /*theTemperature*/)
138{
139//
140// Initialise the result.
141 G4double result = 0.0;
142//
143//
144// Get details of the projectile and target (nucleon number, atomic number,
145// kinetic enery and energy/nucleon.
146//
147 const G4double AT = AA;
148 const G4double ZT = ZZ;
149 const G4double EA = theProjectile->GetKineticEnergy()/MeV;
150 const G4double AP = theProjectile->GetDefinition()->GetBaryonNumber();
151 const G4double ZP = theProjectile->GetDefinition()->GetPDGCharge();
152 G4double E = EA / AP;
153//
154//
155// Determine target mass and energy within the centre-of-mass frame.
156//
157 G4double mT = G4ParticleTable::GetParticleTable()
158 ->GetIonTable()
159 ->GetIonMass(static_cast<G4int>(ZT), static_cast<G4int>(AT));
160 G4LorentzVector pT(0.0, 0.0, 0.0, mT);
161 G4LorentzVector pP(theProjectile->Get4Momentum());
162 pT = pT + pP;
163 G4double E_cm = (pT.mag()-mT-pP.m())/MeV;
164//
165//
166// Determine nuclear radii. Note that the r_p and r_T are defined differently
167// from Wilson et al.
168//
169 G4WilsonRadius theWilsonNuclearRadius;
170 G4double r_rms_p = theWilsonRadius->GetWilsonRMSRadius(AP);
171 G4double r_rms_t = theWilsonRadius->GetWilsonRMSRadius(AT);
172
173 G4double r_p = 1.29*r_rms_p;
174 G4double r_t = 1.29*r_rms_t;
175
176 G4double Radius = (r_p + r_t)/fermi + 1.2*(std::pow(AT, third) + std::pow(AP, third))/
177 std::pow(E_cm, third);
178
179 G4double B = 1.44 * ZP * ZT / Radius;
180//
181//
182// Now determine other parameters associated with the parametric
183// formula, depending upon the projectile and target.
184//
185 G4double T1 = 0.0;
186 G4double D = 0.0;
187 G4double G = 0.0;
188
189 if ((AT==1 && ZT==1) || (AP==1 && ZP==1))
190 {
191 T1 = 23.0;
192 D = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
193 }
194 else if ((AT==1 && ZT==0) || (AP==1 && ZP==0))
195 {
196 T1 = 18.0;
197 D = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
198 }
199 else if ((AT==2 && ZT==1) || (AP==2 && ZP==1))
200 {
201 T1 = 23.0;
202 D = 1.65 + 0.1/(1+std::exp((500.0-E)/200.0));
203 }
204 else if ((AT==3 && ZT==2) || (AP==3 && ZP==2))
205 {
206 T1 = 40.0;
207 D = 1.55;
208 }
209 else if (AP==4 && ZP==2)
210 {
211 if (AT==4 && ZT==2) {T1 = 40.0; G = 300.0;}
212 else if (ZT==4) {T1 = 25.0; G = 300.0;}
213 else if (ZT==7) {T1 = 40.0; G = 500.0;}
214 else if (ZT==13) {T1 = 25.0; G = 300.0;}
215 else if (ZT==26) {T1 = 40.0; G = 300.0;}
216 else {T1 = 40.0; G = 75.0;}
217 D = 2.77 - 8.0E-3*AT + 1.8E-5*AT*AT-0.8/(1.0+std::exp((250.0-E)/G));
218 }
219 else if (AT==4 && ZT==2)
220 {
221 if (AP==4 && ZP==2) {T1 = 40.0; G = 300.0;}
222 else if (ZP==4) {T1 = 25.0; G = 300.0;}
223 else if (ZP==7) {T1 = 40.0; G = 500.0;}
224 else if (ZP==13) {T1 = 25.0; G = 300.0;}
225 else if (ZP==26) {T1 = 40.0; G = 300.0;}
226 else {T1 = 40.0; G = 75.0;}
227 D = 2.77 - 8.0E-3*AP + 1.8E-5*AP*AP-0.8/(1.0+std::exp((250.0-E)/G));
228 }
229//
230//
231// C_E, S, deltaE, X1, S_L and X_m correspond directly with the original
232// formulae of Tripathi et al in his report.
233//
234 G4double C_E = D*(1.0-std::exp(-E/T1)) -
235 0.292*std::exp(-E/792.0)*std::cos(0.229*std::pow(E,0.453));
236
237 G4double S = std::pow(AP,third)*std::pow(AT,third)/(std::pow(AP,third) + std::pow(AT,third));
238
239 G4double deltaE = 0.0;
240 G4double X1 = 0.0;
241 if (AT >= AP)
242 {
243 deltaE = 1.85*S + 0.16*S/std::pow(E_cm,third) - C_E + 0.91*(AT-2.0*ZT)*ZP/AT/AP;
244 X1 = 2.83 - 3.1E-2*AT + 1.7E-4*AT*AT;
245 }
246 else
247 {
248 deltaE = 1.85*S + 0.16*S/std::pow(E_cm,third) - C_E + 0.91*(AP-2.0*ZP)*ZT/AT/AP;
249 X1 = 2.83 - 3.1E-2*AP + 1.7E-4*AP*AP;
250 }
251 G4double S_L = 1.2 + 1.6*(1.0-std::exp(-E/15.0));
252//JMQ 241110 bug fixed
253// G4double X_m = 1.0 - X1*std::exp(-E/X1*S_L);
254 G4double X_m = 1.0 - X1*std::exp(-E/(X1*S_L));
255//
256//
257// R_c is also highly dependent upon the A and Z of the projectile and
258// target.
259//
260 G4double R_c = 1.0;
261 if (AP==1 && ZP==1)
262 {
263 if (AT==2 && ZT==1) R_c = 13.5;
264 else if (AT==3 && ZT==2) R_c = 21.0;
265 else if (AT==4 && ZT==2) R_c = 27.0;
266 else if (ZT==3) R_c = 2.2;
267 }
268 else if (AT==1 && ZT==1)
269 {
270 if (AP==2 && ZP==1) R_c = 13.5;
271 else if (AP==3 && ZP==2) R_c = 21.0;
272 else if (AP==4 && ZP==2) R_c = 27.0;
273 else if (ZP==3) R_c = 2.2;
274 }
275 else if (AP==2 && ZP==1)
276 {
277 if (AT==2 && ZT==1) R_c = 13.5;
278 else if (AT==4 && ZT==2) R_c = 13.5;
279 else if (AT==12 && ZT==6) R_c = 6.0;
280 }
281 else if (AT==2 && ZT==1)
282 {
283 if (AP==2 && ZP==1) R_c = 13.5;
284 else if (AP==4 && ZP==2) R_c = 13.5;
285 else if (AP==12 && ZP==6) R_c = 6.0;
286 }
287 else if ((AP==4 && ZP==2 && (ZT==73 || ZT==79)) ||
288 (AT==4 && ZT==2 && (ZP==73 || ZP==79))) R_c = 0.6;
289//
290//
291// Find the total cross-section. Check that it's value is positive, and if
292// the energy is less that 10 MeV/nuc, find out if the cross-section is
293// increasing with decreasing energy. If so this is a sign that the function
294// is behaving badly at low energies, and the cross-section should be
295// set to zero.
296//
297 result = pi * r_0*r_0 *
298 std::pow((std::pow(AT,third) + std::pow(AP,third) + deltaE),2.0) *
299 (1.0 - R_c*B/E_cm) * X_m;
300 if (!lowEnergyCheck)
301 {
302 if (result < 0.0)
303 result = 0.0;
304 else if (E < 6.0*MeV)
305 {
306 G4double f = 0.95;
307 G4DynamicParticle slowerProjectile = *theProjectile;
308 slowerProjectile.SetKineticEnergy(f * EA * MeV);
309 // G4TripathiLightCrossSection theTripathiLightCrossSection; // MHM 20090824 Not needed
310 // theTripathiLightCrossSection.SetLowEnergyCheck(true);
311 G4bool savelowenergy=lowEnergyCheck;
312 SetLowEnergyCheck(true);
313 G4double resultp =
314 GetZandACrossSection(&slowerProjectile, ZZ, AA, 0.0);
315 SetLowEnergyCheck(savelowenergy);
316 if (resultp >result) result = 0.0;
317 }
318 }
319
320 return result;
321}
322
323
324G4double G4TripathiLightCrossSection::GetCrossSection
325 (const G4DynamicParticle* theProjectile, const G4Element* theTarget,
326 G4double theTemperature)
327{
328 G4int nIso = theTarget->GetNumberOfIsotopes();
329 G4double xsection = 0;
330
331 if (nIso) {
332 G4double sig;
333 G4IsotopeVector* isoVector = theTarget->GetIsotopeVector();
334 G4double* abundVector = theTarget->GetRelativeAbundanceVector();
335 G4int ZZ;
336 G4int AA;
337
338 for (G4int i = 0; i < nIso; i++) {
339 ZZ = (*isoVector)[i]->GetZ();
340 AA = (*isoVector)[i]->GetN();
341 sig = GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
342 xsection += sig*abundVector[i];
343 }
344
345 } else {
346 G4int ZZ = G4lrint(theTarget->GetZ());
347 G4int AA = G4lrint(theTarget->GetN());
348 xsection = GetZandACrossSection(theProjectile, ZZ, AA, theTemperature);
349 }
350
351 return xsection;
352}
353
354
355///////////////////////////////////////////////////////////////////////////////
356//
357void G4TripathiLightCrossSection::SetLowEnergyCheck (G4bool aLowEnergyCheck)
358{
359 lowEnergyCheck = aLowEnergyCheck;
360}
361
Note: See TracBrowser for help on using the repository browser.