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

Last change on this file since 1315 was 1228, checked in by garnier, 16 years ago

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