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

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

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