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

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