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

Last change on this file since 1347 was 1347, checked in by garnier, 14 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.