G4ContinuumGammaTransition.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //      GEANT 4 class file 
00030 //
00031 //      CERN, Geneva, Switzerland
00032 //
00033 //      File name:     G4ContinuumGammaTransition
00034 //
00035 //      Authors:       Carlo Dallapiccola (dallapiccola@umdhep.umd.edu)
00036 //                     Maria Grazia Pia (pia@genova.infn.it)
00037 // 
00038 //      Creation date: 23 October 1998
00039 //
00040 //      Modifications: 
00041 //      
00042 //        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
00043 //              Added creation time evaluation for products of evaporation
00044 //        02 May 2003,   V. Ivanchenko change interface to G4NuclearlevelManager
00045 //        06 Oct 2010,   M. Kelsey -- follow changes to G4NuclearLevelManager
00046 //        17 Nov 2010,   V. Ivanchenko use exponential law for sampling of time
00047 //                                     and extra cleanup
00048 // ----------------------------------------------------------------------------
00049 //
00050 //  Class G4ContinuumGammaTransition.cc
00051 //
00052 
00053 #include "G4ContinuumGammaTransition.hh"
00054 #include "G4VLevelDensityParameter.hh"
00055 #include "G4ConstantLevelDensityParameter.hh"
00056 #include "G4RandGeneralTmp.hh"
00057 #include "G4PhysicalConstants.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "Randomize.hh"
00060 #include "G4Pow.hh"
00061 
00062 //
00063 // Constructor
00064 //
00065 
00066 G4ContinuumGammaTransition::G4ContinuumGammaTransition(
00067                             const G4NuclearLevelManager* levelManager,
00068                             G4int Z, G4int A,
00069                             G4double excitation,
00070                             G4int verbose):
00071   _nucleusA(A), _nucleusZ(Z), _excitation(excitation), _levelManager(levelManager) 
00072 {
00073   G4double eTolerance = 0.;
00074   G4int lastButOne = _levelManager->NumberOfLevels() - 2;
00075   if (lastButOne >= 0)
00076     {
00077       eTolerance = (_levelManager->MaxLevelEnergy() -
00078                     _levelManager->GetLevel(lastButOne)->Energy());
00079       if (eTolerance < 0.) eTolerance = 0.;
00080     }
00081   
00082 
00083   _verbose = verbose;
00084   _eGamma = 0.;
00085   _gammaCreationTime = 0.;
00086 
00087   _maxLevelE = _levelManager->MaxLevelEnergy() + eTolerance;
00088   _minLevelE = _levelManager->MinLevelEnergy();
00089 
00090   // Energy range for photon generation; upper limit is defined 5*Gamma(GDR) from GDR peak
00091   _eMin = 0.001 * MeV;
00092   // Giant Dipole Resonance energy
00093   G4double energyGDR = (40.3 / G4Pow::GetInstance()->powZ(_nucleusA,0.2) ) * MeV;
00094   // Giant Dipole Resonance width
00095   G4double widthGDR = 0.30 * energyGDR;
00096   // Extend 
00097   G4double factor = 5;
00098   _eMax = energyGDR + factor * widthGDR;
00099   if (_eMax > excitation) _eMax = _excitation;
00100 
00101 }
00102 
00103 //
00104 // Destructor
00105 //
00106 
00107 G4ContinuumGammaTransition::~G4ContinuumGammaTransition() 
00108 {}
00109 
00110 void G4ContinuumGammaTransition::SelectGamma()
00111 {
00112 
00113   _eGamma = 0.;
00114 
00115   G4int nBins = 200;
00116   G4double sampleArray[200];
00117   G4int i;
00118   for (i=0; i<nBins; i++)
00119     {
00120       G4double e = _eMin + ( (_eMax - _eMin) / nBins) * i;
00121       sampleArray[i] = E1Pdf(e);
00122 
00123       if(_verbose > 10)
00124         G4cout << "*---* G4ContinuumTransition: e = " << e 
00125                << " pdf = " << sampleArray[i] << G4endl;
00126     }
00127   G4RandGeneralTmp randGeneral(sampleArray, nBins);
00128   G4double random = randGeneral.shoot();
00129   
00130   _eGamma = _eMin + (_eMax - _eMin) * random;
00131   
00132   G4double finalExcitation = _excitation - _eGamma;
00133   
00134   if(_verbose > 10) {
00135     G4cout << "*---*---* G4ContinuumTransition: eGamma = " << _eGamma
00136            << "   finalExcitation = " << finalExcitation 
00137            << " random = " << random << G4endl;
00138   }
00139   //  if (finalExcitation < 0)
00140   if(finalExcitation < _minLevelE/2.)
00141     {
00142       _eGamma = _excitation;
00143       finalExcitation = 0.;
00144     }
00145   
00146   if (finalExcitation < _maxLevelE && finalExcitation > 0.) 
00147     {
00148       G4double levelE = _levelManager->NearestLevel(finalExcitation)->Energy();
00149       G4double diff = finalExcitation - levelE;
00150       _eGamma = _eGamma + diff;
00151     }  
00152 
00153   _gammaCreationTime = GammaTime();
00154 
00155   if(_verbose > 10) {
00156     G4cout << "*---*---* G4ContinuumTransition: _gammaCreationTime = "
00157            << _gammaCreationTime/second << G4endl;
00158   }
00159   return;  
00160 }
00161 
00162 G4double G4ContinuumGammaTransition::GetGammaEnergy()
00163 {
00164   return _eGamma;
00165 }
00166 
00167 G4double G4ContinuumGammaTransition::GetGammaCreationTime()
00168 {
00169   return _gammaCreationTime;
00170 }
00171 
00172 
00173 void G4ContinuumGammaTransition::SetEnergyFrom(G4double energy)
00174 {
00175   if (energy > 0.) _excitation = energy;
00176 }
00177 
00178 
00179 G4double G4ContinuumGammaTransition::E1Pdf(G4double e)
00180 {
00181   G4double theProb = 0.0;
00182   G4double U = std::max(0.0, _excitation - e);
00183 
00184   if(e < 0.0 || _excitation < 0.0) { return theProb; }
00185 
00186   G4ConstantLevelDensityParameter ldPar;
00187   G4double aLevelDensityParam = 
00188     ldPar.LevelDensityParameter(_nucleusA,_nucleusZ,_excitation);
00189 
00190   //G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*_excitation));
00191   //G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(_excitation - e)));
00192   G4double coeff = std::exp(2.0*(std::sqrt(aLevelDensityParam*U) 
00193                                  - std::sqrt(aLevelDensityParam*_excitation)));
00194 
00195   //if(_verbose > 20)
00196   //  G4cout << _nucleusA << " LevelDensityParameter = " <<  aLevelDensityParam
00197   //       << " Bef Aft " << levelDensBef << " " << levelDensAft << G4endl;
00198   
00199   // Now form the probability density
00200 
00201   // Define constants for the photoabsorption cross-section (the reverse
00202   // process of our de-excitation)
00203 
00204   //  G4double sigma0 = 2.5 * _nucleusA * millibarn;  
00205   G4double sigma0 = 2.5 * _nucleusA;  
00206 
00207   G4double Egdp = (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
00208   G4double GammaR = 0.30 * Egdp;
00209  
00210   const G4double normC = 1.0 / (pi * hbarc)*(pi * hbarc);
00211 
00212   G4double numerator = sigma0 * e*e * GammaR*GammaR;
00213   G4double denominator = (e*e - Egdp*Egdp)* (e*e - Egdp*Egdp) + GammaR*GammaR*e*e;
00214   //  if (denominator < 1.0e-9) denominator = 1.0e-9;
00215 
00216   G4double sigmaAbs = numerator/denominator ; 
00217 
00218   if(_verbose > 20) {
00219     G4cout << ".. " << Egdp << " .. " << GammaR 
00220            << " .. " << normC << " .. " << sigmaAbs  
00221            << " .. " << e*e << " .. " << coeff
00222            << G4endl;
00223   }
00224 
00225   //  theProb = normC * sigmaAbs * e*e * levelDensAft/levelDensBef;
00226   theProb =  sigmaAbs * e*e * coeff;
00227 
00228   return theProb;
00229 }
00230 
00231 
00232 G4double G4ContinuumGammaTransition::GammaTime()
00233 {
00234 
00235   G4double GammaR = 0.30 * (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
00236   G4double tau = hbar_Planck/GammaR;
00237   G4double creationTime = -tau*std::log(G4UniformRand());
00238   /*
00239   G4double tMin = 0;
00240   G4double tMax = 10.0 * tau;
00241   G4int nBins = 200;
00242   G4double sampleArray[200];
00243 
00244   for(G4int i = 0; i<nBins;i++)
00245   {
00246     G4double t = tMin + ((tMax-tMin)/nBins)*i;
00247     sampleArray[i] = (std::exp(-t/tau))/tau;
00248   }
00249 
00250   G4RandGeneralTmp randGeneral(sampleArray, nBins);
00251   G4double random = randGeneral.shoot();
00252   
00253   G4double creationTime = tMin + (tMax - tMin) * random;
00254   */
00255   return creationTime;
00256 }
00257 
00258 
00259 
00260 
00261 
00262 
00263 
00264 
00265 
00266 
00267 
00268 

Generated on Mon May 27 17:47:57 2013 for Geant4 by  doxygen 1.4.7