Changeset 1347 for trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaTransition.cc

Timestamp:

Dec 22, 2010, 3:52:27 PM (13 years ago)

Author:

garnier

Message:

geant4 tag 9.4

File:

• trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaTransition.cc (modified) (14 diffs)

trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4ContinuumGammaTransition.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4ContinuumGammaTransition.cc,v 1.4 2010/10/07 07:50:13 mkelsey Exp $
+// $Id: G4ContinuumGammaTransition.cc,v 1.5 2010/11/17 19:17:17 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-ref-00 $
+//
 // -------------------------------------------------------------------
 //      GEANT 4 class file 
@@ -41 +43 @@
 //        15 April 1999, Alessandro Brunengo (Alessandro.Brunengo@ge.infn.it)
 //              Added creation time evaluation for products of evaporation
-//        02 May 2003,   Vladimir Ivanchenko change interface to G4NuclearlevelManager
-//        06 Oct 2010, M. Kelsey -- follow changes to G4NuclearLevelManager
-// -------------------------------------------------------------------
+//        02 May 2003,   V. Ivanchenko change interface to G4NuclearlevelManager
+//        06 Oct 2010,   M. Kelsey -- follow changes to G4NuclearLevelManager
+//        17 Nov 2010,   V. Ivanchenko use exponential law for sampling of time
+//                                     and extra cleanup
+// ----------------------------------------------------------------------------
 //
 //  Class G4ContinuumGammaTransition.cc
@@ -52 +56 @@
 #include "G4ConstantLevelDensityParameter.hh"
 #include "G4RandGeneralTmp.hh"
+#include "Randomize.hh"
+#include "G4Pow.hh"
+
 //
 // Constructor
@@ -83 +90 @@
   _eMin = 0.001 * MeV;
   // Giant Dipole Resonance energy
-  G4double energyGDR = (40.3 / std::pow(G4double(_nucleusA),0.2) ) * MeV;
+  G4double energyGDR = (40.3 / G4Pow::GetInstance()->powZ(_nucleusA,0.2) ) * MeV;
   // Giant Dipole Resonance width
   G4double widthGDR = 0.30 * energyGDR;
@@ -97 +104 @@
 //
 
-G4ContinuumGammaTransition::~G4ContinuumGammaTransition() {}
-
-//
-// Override GammaEnergy function from G4VGammaTransition
-//
+G4ContinuumGammaTransition::~G4ContinuumGammaTransition() 
+{}
 
 void G4ContinuumGammaTransition::SelectGamma()
@@ -127 +131 @@
   G4double finalExcitation = _excitation - _eGamma;
   
-  if(_verbose > 10)
+  if(_verbose > 10) {
     G4cout << "*---*---* G4ContinuumTransition: eGamma = " << _eGamma
            << "   finalExcitation = " << finalExcitation 
            << " random = " << random << G4endl;
-
-//  if (finalExcitation < 0)
+  }
+  //  if (finalExcitation < 0)
   if(finalExcitation < _minLevelE/2.)
     {
@@ -148 +152 @@
   _gammaCreationTime = GammaTime();
 
-  if(_verbose > 10)
+  if(_verbose > 10) {
     G4cout << "*---*---* G4ContinuumTransition: _gammaCreationTime = "
            << _gammaCreationTime/second << G4endl;
-
+  }
   return;  
 }
@@ -166 +170 @@
 
 
-void G4ContinuumGammaTransition::SetEnergyFrom(const G4double energy)
-{
-
+void G4ContinuumGammaTransition::SetEnergyFrom(G4double energy)
+{
   if (energy > 0.) _excitation = energy;
-  return;  
-
 }
 
@@ -178 +179 @@
 {
   G4double theProb = 0.0;
-
-  if( (_excitation - e) < 0.0 || e < 0 || _excitation < 0) return theProb;
+  G4double U = std::max(0.0, _excitation - e);
+
+  if(e < 0.0 || _excitation < 0.0) { return theProb; }
 
   G4ConstantLevelDensityParameter ldPar;
-  G4double aLevelDensityParam = ldPar.LevelDensityParameter(_nucleusA,_nucleusZ,_excitation);
-
-  G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*_excitation));
-  G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(_excitation - e)));
-
-  if(_verbose > 20)
-    G4cout << _nucleusA << " LevelDensityParameter = " <<  aLevelDensityParam
-           << " Bef Aft " << levelDensBef << " " << levelDensAft << G4endl;
+  G4double aLevelDensityParam = 
+    ldPar.LevelDensityParameter(_nucleusA,_nucleusZ,_excitation);
+
+  //G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*_excitation));
+  //G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(_excitation - e)));
+  G4double coeff = std::exp(2.0*(std::sqrt(aLevelDensityParam*U) 
+                                 - std::sqrt(aLevelDensityParam*_excitation)));
+
+  //if(_verbose > 20)
+  //  G4cout << _nucleusA << " LevelDensityParameter = " <<  aLevelDensityParam
+  //       << " Bef Aft " << levelDensBef << " " << levelDensAft << G4endl;
   
   // Now form the probability density
@@ -199 +204 @@
   G4double sigma0 = 2.5 * _nucleusA;  
 
-  G4double Egdp = (40.3 / std::pow(G4double(_nucleusA),0.2) )*MeV;
+  G4double Egdp = (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
   G4double GammaR = 0.30 * Egdp;
  
-  G4double normC = 1.0 / (pi * hbarc)*(pi * hbarc);
+  const G4double normC = 1.0 / (pi * hbarc)*(pi * hbarc);
 
   G4double numerator = sigma0 * e*e * GammaR*GammaR;
@@ -210 +215 @@
   G4double sigmaAbs = numerator/denominator ; 
 
-  if(_verbose > 20)
+  if(_verbose > 20) {
     G4cout << ".. " << Egdp << " .. " << GammaR 
            << " .. " << normC << " .. " << sigmaAbs  
-           << " .. " << e*e << " .. " << levelDensAft/levelDensBef
+           << " .. " << e*e << " .. " << coeff
            << G4endl;
+  }
 
   //  theProb = normC * sigmaAbs * e*e * levelDensAft/levelDensBef;
-  theProb =  sigmaAbs * e*e * levelDensAft/levelDensBef;
+  theProb =  sigmaAbs * e*e * coeff;
 
   return theProb;
@@ -226 +232 @@
 {
 
-  G4double GammaR = 0.30 * (40.3 / std::pow(G4double(_nucleusA),0.2) )*MeV;
+  G4double GammaR = 0.30 * (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
   G4double tau = hbar_Planck/GammaR;
-
+  G4double creationTime = -tau*std::log(G4UniformRand());
+  /*
   G4double tMin = 0;
   G4double tMax = 10.0 * tau;
@@ -234 +241 @@
   G4double sampleArray[200];
 
-  for(G4int i = 0;i<nBins;i++)
+  for(G4int i = 0; i<nBins;i++)
   {
     G4double t = tMin + ((tMax-tMin)/nBins)*i;
@@ -244 +251 @@
   
   G4double creationTime = tMin + (tMax - tMin) * random;
-
+  */
   return creationTime;
 }