Geant4  10.01
G4ContinuumGammaTransition.cc
Go to the documentation of this file.
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 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4ContinuumGammaTransition.cc 87443 2014-12-04 12:26:31Z gunter $
27 //
28 // -------------------------------------------------------------------
29 // GEANT 4 class file
30 //
31 // CERN, Geneva, Switzerland
32 //
33 // File name: G4ContinuumGammaTransition
34 //
35 // Authors: Carlo Dallapiccola (dallapiccola@umdhep.umd.edu)
36 // Maria Grazia Pia (pia@genova.infn.it)
37 //
38 // Creation date: 23 October 1998
39 //
40 // Modifications:
41 //
42 // 15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
43 // Added creation time evaluation for products of evaporation
44 // 02 May 2003, V. Ivanchenko change interface to G4NuclearlevelManager
45 // 06 Oct 2010, M. Kelsey -- follow changes to G4NuclearLevelManager
46 // 17 Nov 2010, V. Ivanchenko use exponential law for sampling of time
47 // and extra cleanup
48 // ----------------------------------------------------------------------------
49 //
50 // Class G4ContinuumGammaTransition.cc
51 //
52 
54 #include "G4NuclearLevelManager.hh"
55 #include "G4PhysicalConstants.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "Randomize.hh"
58 #include "G4Pow.hh"
59 #include "G4Log.hh"
60 
61 //static const G4double normC = CLHEP::millibarn/
62 // (CLHEP::pi*CLHEP::hbarc*CLHEP::pi*CLHEP::hbarc);
63 static const G4double factor = 5;
64 static const G4double tolerance = 2*CLHEP::keV;
65 
67  const G4NuclearLevelManager* manager,
68  G4int Z, G4int A, G4double exc, G4int verb)
69 {
70  nBins = 100;
71  verbose = verb;
72  eMin = keV;
74  sampleArray.resize(nBins+1,0.0);
76  Update(manager, Z, A, exc);
77 }
78 
80 {}
81 
83  G4int Z, G4int A, G4double exc)
84 {
85  levelManager = manager;
86  nucleusZ = Z;
87  nucleusA = A;
88  excitation = exc;
89  eGamma = 0.;
90  gammaCreationTime = 0.;
91 
93  if(levelManager) {
95  }
96  // Energy range for photon generation;
97  // upper limit is defined 5*Gamma(GDR) from GDR peak
98  // Giant Dipole Resonance energy
99  G4double energyGDR = (40.3 / g4pow->powZ(nucleusA,0.2) ) * MeV;
100  energyGDR2 = energyGDR*energyGDR;
101  // Giant Dipole Resonance width
102  widthGDR = 0.30 * energyGDR;
104  // Extend
105  eMax = energyGDR + factor * widthGDR;
106  if (eMax > excitation) { eMax = excitation; }
107 }
108 
110 {
111  eGamma = 0.;
112  sampleArray[0] = 0.0;
113  G4int i;
114  G4double del = (eMax - eMin) / G4double(nBins);
115  G4double sum = 0;
116  G4double w1 = E1Pdf(eMin);
117  G4double w2;
118  //G4cout << eMin << " " << eMax << " " << del << G4endl;
119  for (i=1; i<=nBins; i++) {
120  G4double e = eMin + del * i;
121  w2 = E1Pdf(e);
122  sum += 0.5*(w1 + w2);
123  w1 = w2;
124  sampleArray[i] = sum;
125  if(verbose > 1) {
126  G4cout << "*---* G4ContinuumTransition: e = " << e
127  << " pdf = " << sampleArray[i] << G4endl;
128  }
129  }
130  sum *= G4UniformRand();
131  eGamma = eMax;
132  for (i=1; i<=nBins; i++) {
133  if(sum <= sampleArray[i]) {
134  eGamma = eMin + del * i;
135  G4double w = sampleArray[i] - sampleArray[i-1];
136  //G4cout << eGamma << " " << w << G4endl;
137  if(w != 0.0) {
138  eGamma -= (sampleArray[i] - sum)*del/w;
139  }
140  break;
141  }
142  }
143 
144  G4double finalExcitation = excitation - eGamma;
145 
146  if(verbose > 1) {
147  G4cout << "*---*---* G4ContinuumTransition: eGamma = " << eGamma
148  << " finalExcitation = " << finalExcitation << G4endl;
149  }
150 
151  if(finalExcitation <= minLevelE) {
152  eGamma = excitation;
153 
154  } else {
155  finalExcitation = levelManager->NearestLevel(finalExcitation)->Energy();
156  eGamma = excitation - finalExcitation;
157  }
158 
160 
161  if(verbose > 1) {
162  G4cout << "*---*---* G4ContinuumTransition: gammaCreationTime = "
164  }
165 }
166 
168 {
169  return eGamma;
170 }
171 
173 {
174  return gammaCreationTime;
175 }
176 
178 {
179  excitation = energy;
180 }
181 
183 {
184  G4double theProb = 0.0;
185  G4double U = excitation - e;
186 
187  if(U < 0.0) { return theProb; }
188 
189  G4double aLevelDensityParam =
191 
192  //G4double levelDensBef = G4Exp(2.0*std::sqrt(aLevelDensityParam*excitation));
193  //G4double levelDensAft = G4Exp(2.0*std::sqrt(aLevelDensityParam*(excitation - e)));
194  G4double coeff = G4Exp(2.0*(std::sqrt(aLevelDensityParam*U)
195  - std::sqrt(aLevelDensityParam*excitation)));
196 
197  //if(verbose > 20)
198  // G4cout << nucleusA << " LevelDensityParameter = " << aLevelDensityParam
199  // << " Bef Aft " << levelDensBef << " " << levelDensAft << G4endl;
200 
201  // Now form the probability density
202  // Define constants for the photoabsorption cross-section (the reverse
203  // process of our de-excitation)
204 
205  G4double sigma0 = 2.5 * nucleusA;
206 
207  G4double e2 = e*e;
208  G4double numerator = sigma0 * e2 * widthGDR2;
209  G4double denominator = (e2 - energyGDR2)* (e2 - energyGDR2) + widthGDR2*e2;
210  // if (denominator < 1.0e-9) denominator = 1.0e-9;
211 
212  G4double sigmaAbs = numerator/denominator;
213 
214  if(verbose > 2) {
215  G4cout << "E_GDR(MeV)= " << std::sqrt(energyGDR2) << " W_GDR(MeV)= " << widthGDR
216  << " sigAbs= " << sigmaAbs
217  << " E(MeV)= " << e << " coeff= " << coeff
218  << G4endl;
219  }
220 
221  theProb = sigmaAbs * e2 * coeff;
222 
223  return theProb;
224 }
225 
227 {
228  G4double tau = hbar_Planck/widthGDR;
229  G4double creationTime = -tau*G4Log(G4UniformRand());
230 
231  return creationTime;
232 }
233 
234 
235 
236 
237 
238 
239 
240 
241 
242 
243 
244 
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static const double MeV
Definition: G4SIunits.hh:193
G4double Energy() const
G4double MinLevelEnergy() const
static const G4double e2
static const G4double tolerance
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
static const double second
Definition: G4SIunits.hh:138
virtual void SetEnergyFrom(G4double energy)
static const G4double A[nN]
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static const G4double factor
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
void Update(const G4NuclearLevelManager *manager, G4int Z, G4int A, G4double exc)
G4ConstantLevelDensityParameter ldPar
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:195
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:256
G4double LevelDensityParameter(const G4int A, const G4int, const G4double) const
double G4double
Definition: G4Types.hh:76
const G4NuclearLevelManager * levelManager
G4ContinuumGammaTransition(const G4NuclearLevelManager *manager, G4int Z, G4int A, G4double exc, G4int ver)
#define DBL_MAX
Definition: templates.hh:83
const G4NuclearLevel * NearestLevel(G4double energy, G4double eDiffMax=1.e+8) const