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

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

update CVS release candidate geant4.9.3.01
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
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  if(E_cm <= DBL_MIN) return result; 
156
157
158//
159//
160// Determine nuclear radii.  Note that the r_p and r_T are defined differently
161// from Wilson et al.
162//
163  G4WilsonRadius theWilsonNuclearRadius;
164  G4double r_rms_p = theWilsonRadius->GetWilsonRMSRadius(AP);
165  G4double r_rms_t = theWilsonRadius->GetWilsonRMSRadius(AT);
166
167  G4double r_p = 1.29*r_rms_p;
168  G4double r_t = 1.29*r_rms_t;
169
170  G4double Radius = (r_p + r_t)/fermi + 1.2*(std::pow(AT, third) + std::pow(AP, third))/
171    std::pow(E_cm, third);
172
173  G4double B = 1.44 * ZP * ZT / Radius;
174  if(E_cm <= B) return result; 
175
176//
177// Now determine other parameters associated with the parametric
178// formula, depending upon the projectile and target.
179//
180  G4double T1 = 0.0;
181  G4double D  = 0.0;
182  G4double G  = 0.0;
183
184  if ((AT==1 && ZT==1) || (AP==1 && ZP==1))
185  {
186    T1 = 23.0;
187    D  = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
188  }
189  else if ((AT==1 && ZT==0) || (AP==1 && ZP==0))
190  {
191    T1 = 18.0;
192    D  = 1.85 + 0.16/(1+std::exp((500.0-E)/200.0));
193  }
194  else if ((AT==2 && ZT==1) || (AP==2 && ZP==1))
195  {
196    T1 = 23.0;
197    D  = 1.65 + 0.1/(1+std::exp((500.0-E)/200.0));
198  }
199  else if ((AT==3 && ZT==2) || (AP==3 && ZP==2))
200  {
201    T1 = 40.0;
202    D  = 1.55;
203  }
204  else if (AP==4 && ZP==2)
205  {
206    if      (AT==4 && ZT==2) {T1 = 40.0; G = 300.0;}
207    else if (ZT==4)          {T1 = 25.0; G = 300.0;}
208    else if (ZT==7)          {T1 = 40.0; G = 500.0;}
209    else if (ZT==13)         {T1 = 25.0; G = 300.0;}
210    else if (ZT==26)         {T1 = 40.0; G = 300.0;}
211    else                     {T1 = 40.0; G = 75.0;}
212    D = 2.77 - 8.0E-3*AT + 1.8E-5*AT*AT-0.8/(1.0+std::exp((250.0-E)/G));
213  }
214  else if (AT==4 && ZT==2)
215  {
216    if      (AP==4 && ZP==2) {T1 = 40.0; G = 300.0;}
217    else if (ZP==4)          {T1 = 25.0; G = 300.0;}
218    else if (ZP==7)          {T1 = 40.0; G = 500.0;}
219    else if (ZP==13)         {T1 = 25.0; G = 300.0;}
220    else if (ZP==26)         {T1 = 40.0; G = 300.0;}
221    else                     {T1 = 40.0; G = 75.0;}
222    D = 2.77 - 8.0E-3*AP + 1.8E-5*AP*AP-0.8/(1.0+std::exp((250.0-E)/G));
223  }
224//
225//
226// C_E, S, deltaE, X1, S_L and X_m correspond directly with the original
227// formulae of Tripathi et al in his report.
228//
229  G4double C_E = D*(1.0-std::exp(-E/T1)) -
230                 0.292*std::exp(-E/792.0)*std::cos(0.229*std::pow(E,0.453));
231
232  G4double S = std::pow(AP,third)*std::pow(AT,third)/(std::pow(AP,third) + std::pow(AT,third));
233
234  G4double deltaE = 0.0;
235  G4double X1     = 0.0;
236  if (AT >= AP)
237  {
238    deltaE = 1.85*S + 0.16*S/std::pow(E_cm,third) - C_E + 0.91*(AT-2.0*ZT)*ZP/AT/AP;
239    X1     = 2.83 - 3.1E-2*AT + 1.7E-4*AT*AT;
240  }
241  else
242  {
243    deltaE = 1.85*S + 0.16*S/std::pow(E_cm,third) - C_E + 0.91*(AP-2.0*ZP)*ZT/AT/AP;
244    X1     = 2.83 - 3.1E-2*AP + 1.7E-4*AP*AP;
245  }
246  G4double S_L = 1.2 + 1.6*(1.0-std::exp(-E/15.0));
247  G4double X_m = 1.0 - X1*std::exp(-E/X1*S_L);
248//
249//
250// R_c is also highly dependent upon the A and Z of the projectile and
251// target.
252//
253  G4double R_c = 1.0;
254  if (AP==1 && ZP==1)
255  {
256    if      (AT==2 && ZT==1) R_c = 13.5;
257    else if (AT==3 && ZT==2) R_c = 21.0;
258    else if (AT==4 && ZT==2) R_c = 27.0;
259    else if (ZT==3)          R_c = 2.2;
260  }
261  else if (AT==1 && ZT==1)
262  {
263    if      (AP==2 && ZP==1) R_c = 13.5;
264    else if (AP==3 && ZP==2) R_c = 21.0;
265    else if (AP==4 && ZP==2) R_c = 27.0;
266    else if (ZP==3)          R_c = 2.2;
267  }
268  else if (AP==2 && ZP==1)
269  {
270    if       (AT==2 && ZT==1) R_c = 13.5;
271    else if (AT==4 && ZT==2)  R_c = 13.5;
272    else if (AT==12 && ZT==6) R_c = 6.0;
273  }
274  else if (AT==2 && ZT==1)
275  {
276    if       (AP==2 && ZP==1) R_c = 13.5;
277    else if (AP==4 && ZP==2)  R_c = 13.5;
278    else if (AP==12 && ZP==6) R_c = 6.0;
279  }
280  else if ((AP==4 && ZP==2 && (ZT==73 || ZT==79)) ||
281           (AT==4 && ZT==2 && (ZP==73 || ZP==79))) R_c = 0.6;
282//
283//
284// Find the total cross-section.  Check that it's value is positive, and if
285// the energy is less that 10 MeV/nuc, find out if the cross-section is
286// increasing with decreasing energy.  If so this is a sign that the function
287// is behaving badly at low energies, and the cross-section should be
288// set to zero.
289//
290  result = pi * r_0*r_0 *
291           std::pow((std::pow(AT,third) + std::pow(AP,third) + deltaE),2.0) *
292           (1.0 - R_c*B/E_cm) * X_m;
293  if(result < 0.0) result = 0.0;
294  /*
295  if (!lowEnergyCheck)
296  {
297    if (result < 0.0)  result = 0.0;
298    else if (E < 6.0*MeV)
299    {
300      G4double f  = 0.95;
301      G4DynamicParticle slowerProjectile = *theProjectile;
302      slowerProjectile.SetKineticEnergy(f * EA * MeV);
303      G4TripathiLightCrossSection theTripathiLightCrossSection;
304      theTripathiLightCrossSection.SetLowEnergyCheck(true);
305      G4double resultp =
306        theTripathiLightCrossSection.GetIsoZACrossSection
307        (&slowerProjectile, ZZ, AA, 0.0);
308      if (resultp >result) result = 0.0;
309    }
310  }
311  */
312  return result;
313}
314
315
316G4double G4TripathiLightCrossSection::GetCrossSection
317  (const G4DynamicParticle* theProjectile, const G4Element* theTarget,
318  G4double theTemperature)
319{
320  G4int nIso = theTarget->GetNumberOfIsotopes();
321  G4double xsection = 0;
322     
323  if (nIso) {
324    G4double sig;
325    G4IsotopeVector* isoVector = theTarget->GetIsotopeVector();
326    G4double* abundVector = theTarget->GetRelativeAbundanceVector();
327    G4double ZZ;
328    G4double AA;
329     
330    for (G4int i = 0; i < nIso; i++) {
331      ZZ = G4double( (*isoVector)[i]->GetZ() );
332      AA = G4double( (*isoVector)[i]->GetN() );
333      sig = GetIsoZACrossSection(theProjectile, ZZ, AA, theTemperature);
334      xsection += sig*abundVector[i];
335    }
336   
337  } else {
338    xsection =
339      GetIsoZACrossSection(theProjectile, theTarget->GetZ(), theTarget->GetN(),
340                           theTemperature);
341  }
342   
343  return xsection;
344}
345
346
347///////////////////////////////////////////////////////////////////////////////
348//
349void G4TripathiLightCrossSection::SetLowEnergyCheck (G4bool aLowEnergyCheck)
350{
351  lowEnergyCheck = aLowEnergyCheck;
352}
353///////////////////////////////////////////////////////////////////////////////
354//
Note: See TracBrowser for help on using the repository browser.