Geant4_10
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 74869 2013-10-23 09:26:17Z gcosmo $
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 
56 #include "G4PhysicalConstants.hh"
57 #include "G4SystemOfUnits.hh"
58 #include "Randomize.hh"
59 #include "G4Pow.hh"
60 #include "G4Log.hh"
61 
62 //
63 // Constructor
64 //
65 
67  const G4NuclearLevelManager* levelManager,
68  G4int Z, G4int A,
69  G4double excitation,
70  G4int verbose):
71  _nucleusA(A), _nucleusZ(Z), _excitation(excitation), _levelManager(levelManager)
72 {
73  G4double eTolerance = 0.;
74  G4int lastButOne = _levelManager->NumberOfLevels() - 2;
75  if (lastButOne >= 0)
76  {
77  eTolerance = (_levelManager->MaxLevelEnergy() -
78  _levelManager->GetLevel(lastButOne)->Energy());
79  if (eTolerance < 0.) eTolerance = 0.;
80  }
81 
82 
83  _verbose = verbose;
84  _eGamma = 0.;
85  _gammaCreationTime = 0.;
86 
87  _maxLevelE = _levelManager->MaxLevelEnergy() + eTolerance;
88  _minLevelE = _levelManager->MinLevelEnergy();
89 
90  // Energy range for photon generation; upper limit is defined 5*Gamma(GDR) from GDR peak
91  _eMin = 0.001 * MeV;
92  // Giant Dipole Resonance energy
93  G4double energyGDR = (40.3 / G4Pow::GetInstance()->powZ(_nucleusA,0.2) ) * MeV;
94  // Giant Dipole Resonance width
95  G4double widthGDR = 0.30 * energyGDR;
96  // Extend
97  G4double factor = 5;
98  _eMax = energyGDR + factor * widthGDR;
99  if (_eMax > excitation) _eMax = _excitation;
100 
101 }
102 
103 //
104 // Destructor
105 //
106 
108 {}
109 
111 {
112 
113  _eGamma = 0.;
114 
115  G4int nBins = 100;
116  G4double sampleArray[101];
117  sampleArray[0] = 0.0;
118  G4int i;
119  G4double del = (_eMax - _eMin) / G4double(nBins);
120  G4double sum = 0;
121  G4double w1 = E1Pdf(_eMin);
122  G4double w2;
123  //G4cout << _eMin << " " << _eMax << " " << del << G4endl;
124  for (i=1; i<=nBins; i++) {
125  G4double e = _eMin + del * i;
126  w2 = E1Pdf(e);
127  sum += 0.5*(w1 + w2);
128  w1 = w2;
129  sampleArray[i] = sum;
130  if(_verbose > 10) {
131  G4cout << "*---* G4ContinuumTransition: e = " << e
132  << " pdf = " << sampleArray[i] << G4endl;
133  }
134  }
135  sum *= G4UniformRand();
136  _eGamma = _eMax;
137  for (i=1; i<=nBins; i++) {
138  if(sum <= sampleArray[i]) {
139  _eGamma = _eMin + del * i;
140  G4double w = sampleArray[i] - sampleArray[i-1];
141  //G4cout << _eGamma << " " << w << G4endl;
142  if(w != 0.0) {
143  _eGamma -= (sampleArray[i] - sum)*del/w;
144  }
145  break;
146  }
147  }
148 
149  G4double finalExcitation = _excitation - _eGamma;
150 
151  if(_verbose > 10) {
152  G4cout << "*---*---* G4ContinuumTransition: eGamma = " << _eGamma
153  << " finalExcitation = " << finalExcitation << G4endl;
154  }
155  // if (finalExcitation < 0)
156  if(finalExcitation < _minLevelE/2.)
157  {
158  _eGamma = _excitation;
159  finalExcitation = 0.;
160  }
161 
162  if (finalExcitation < _maxLevelE && finalExcitation > 0.)
163  {
164  G4double levelE = _levelManager->NearestLevel(finalExcitation)->Energy();
165  G4double diff = finalExcitation - levelE;
166  _eGamma = _eGamma + diff;
167  }
168 
169  _gammaCreationTime = GammaTime();
170 
171  if(_verbose > 10) {
172  G4cout << "*---*---* G4ContinuumTransition: _gammaCreationTime = "
173  << _gammaCreationTime/second << G4endl;
174  }
175  return;
176 }
177 
179 {
180  return _eGamma;
181 }
182 
184 {
185  return _gammaCreationTime;
186 }
187 
188 
190 {
191  if (energy > 0.) _excitation = energy;
192 }
193 
194 
195 G4double G4ContinuumGammaTransition::E1Pdf(G4double e)
196 {
197  G4double theProb = 0.0;
198  G4double U = std::max(0.0, _excitation - e);
199 
200  if(e < 0.0 || _excitation < 0.0) { return theProb; }
201 
203  G4double aLevelDensityParam =
204  ldPar.LevelDensityParameter(_nucleusA,_nucleusZ,_excitation);
205 
206  //G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*_excitation));
207  //G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(_excitation - e)));
208  G4double coeff = std::exp(2.0*(std::sqrt(aLevelDensityParam*U)
209  - std::sqrt(aLevelDensityParam*_excitation)));
210 
211  //if(_verbose > 20)
212  // G4cout << _nucleusA << " LevelDensityParameter = " << aLevelDensityParam
213  // << " Bef Aft " << levelDensBef << " " << levelDensAft << G4endl;
214 
215  // Now form the probability density
216 
217  // Define constants for the photoabsorption cross-section (the reverse
218  // process of our de-excitation)
219 
220  // G4double sigma0 = 2.5 * _nucleusA * millibarn;
221  G4double sigma0 = 2.5 * _nucleusA;
222 
223  G4double Egdp = (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
224  G4double GammaR = 0.30 * Egdp;
225 
226  static const G4double normC = 1.0 / (pi * hbarc)*(pi * hbarc);
227 
228  G4double numerator = sigma0 * e*e * GammaR*GammaR;
229  G4double denominator = (e*e - Egdp*Egdp)* (e*e - Egdp*Egdp) + GammaR*GammaR*e*e;
230  // if (denominator < 1.0e-9) denominator = 1.0e-9;
231 
232  G4double sigmaAbs = numerator/denominator ;
233 
234  if(_verbose > 20) {
235  G4cout << ".. " << Egdp << " .. " << GammaR
236  << " .. " << normC << " .. " << sigmaAbs
237  << " .. " << e*e << " .. " << coeff
238  << G4endl;
239  }
240 
241  // theProb = normC * sigmaAbs * e*e * levelDensAft/levelDensBef;
242  theProb = sigmaAbs * e*e * coeff;
243 
244  return theProb;
245 }
246 
247 
248 G4double G4ContinuumGammaTransition::GammaTime()
249 {
250 
251  G4double GammaR = 0.30 * (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
252  G4double tau = hbar_Planck/GammaR;
253  G4double creationTime = -tau*G4Log(G4UniformRand());
254 
255  return creationTime;
256 }
257 
258 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
static G4Pow * GetInstance()
Definition: G4Pow.cc:53
G4double Energy() const
G4double MinLevelEnergy() const
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
Float_t Z
Definition: plot.C:39
const G4NuclearLevel * GetLevel(G4int i) const
virtual void SetEnergyFrom(G4double energy)
G4double G4Log(G4double x)
Definition: G4Log.hh:227
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
G4ContinuumGammaTransition(const G4NuclearLevelManager *levelManager, G4int Z, G4int A, G4double excitation, G4int verbose)
G4double powZ(G4int Z, G4double y) const
Definition: G4Pow.hh:258
G4double LevelDensityParameter(const G4int A, const G4int, const G4double) const
double G4double
Definition: G4Types.hh:76
const G4NuclearLevel * NearestLevel(G4double energy, G4double eDiffMax=9999.*CLHEP::GeV) const
G4double MaxLevelEnergy() const