Changeset 1315 for trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.cc

Timestamp:

Jun 18, 2010, 11:42:07 AM (14 years ago)

Author:

garnier

Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:

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

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

@@ -25 +25 @@
 //
 //
-//  Class G4E1Probability.cc
+// $Id: G4E1Probability.cc,v 1.9 2010/05/25 10:47:24 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
+//
+//---------------------------------------------------------------------
+//
+// Geant4 class G4E1Probability
+//
+// by V. Lara (May 2003)
+//
+// Modifications:
+// 18.05.2010 V.Ivanchenko trying to speedup the most slow method
+//            by usage of G4Pow, integer A and introduction of const members
+//
 //
 
 #include "G4E1Probability.hh"
-#include "G4ConstantLevelDensityParameter.hh"
+//#include "G4ConstantLevelDensityParameter.hh"
 #include "Randomize.hh"
+#include "G4Pow.hh"
 
 // Constructors and operators
 //
 
-G4E1Probability::G4E1Probability(const G4E1Probability& ) : G4VEmissionProbability()
-{
-
-  throw G4HadronicException(__FILE__, __LINE__, "G4E1Probability::copy_constructor meant to not be accessible");
-
-}
-
-const G4E1Probability& G4E1Probability::
-operator=(const G4E1Probability& ) 
-{
-
-  throw G4HadronicException(__FILE__, __LINE__, "G4E1Probability::operator= meant to not be accessible");
-  return *this;
-}
-
-G4bool G4E1Probability::operator==(const G4E1Probability& ) const
-{
-
-  return false;
-
-}
-
-G4bool G4E1Probability::operator!=(const G4E1Probability& ) const
-{
-
-  return true;
-
-}
+G4E1Probability::G4E1Probability():G4VEmissionProbability()
+{
+  G4double x = CLHEP::pi*CLHEP::hbarc;
+  normC = 1.0 / (x*x);
+  theLevelDensityParameter = 0.125/MeV;
+  fG4pow = G4Pow::GetInstance(); 
+}
+
+G4E1Probability::~G4E1Probability()
+{}
 
 // Calculate the emission probability
@@ -68 +63 @@
 
 G4double G4E1Probability::EmissionProbDensity(const G4Fragment& frag, 
-                                                 const G4double gammaE)
+                                              const G4double gammaE)
 {
 
@@ -76 +71 @@
   // the probability density for photon evaporation from U to U - gammaE
   // (U = nucleus excitation energy, gammaE = total evaporated photon
-  // energy).
-  // fragment = nuclear fragment BEFORE de-excitation
+  // energy). Fragment = nuclear fragment BEFORE de-excitation
 
   G4double theProb = 0.0;
 
-  const G4double Afrag = frag.GetA();
-  const G4double Zfrag = frag.GetZ();
-  const G4double Uexcite = frag.GetExcitationEnergy();
-
-  if( (Uexcite-gammaE) < 0.0 || gammaE < 0 || Uexcite <= 0) return theProb;
+  G4int Afrag = frag.GetA_asInt();
+  G4double Uexcite = frag.GetExcitationEnergy();
+
+  if( (Uexcite-gammaE) < 0.0 || gammaE < 0) { return theProb; }
 
   // Need a level density parameter.
@@ -91 +84 @@
   // nuclei).
 
-  static G4ConstantLevelDensityParameter a;
-
-  G4double aLevelDensityParam = a.LevelDensityParameter(static_cast<G4int>(Afrag),
-                                                        static_cast<G4int>(Zfrag),
-                                                        Uexcite);
-
-  G4double levelDensBef = std::exp(2.0*std::sqrt(aLevelDensityParam*Uexcite));
-  G4double levelDensAft = std::exp(2.0*std::sqrt(aLevelDensityParam*(Uexcite-gammaE)));
-
+  G4double aLevelDensityParam = Afrag*theLevelDensityParameter;
+
+  //  G4double levelDensBef = std::exp(2*std::sqrt(aLevelDensityParam*Uexcite));
+  //  G4double levelDensAft = std::exp(2*std::sqrt(aLevelDensityParam*(Uexcite-gammaE)));
+  // VI reduce number of calls to exp 
+  G4double levelDens = 1.0;
+  if( aLevelDensityParam > 0.0 ) {
+    levelDens = std::exp(2*(std::sqrt(aLevelDensityParam*(Uexcite-gammaE)) 
+      - std::sqrt(aLevelDensityParam*Uexcite))); 
+  }
   // Now form the probability density
 
@@ -107 +101 @@
   G4double sigma0 = 2.5 * Afrag * millibarn;  // millibarns
 
-  G4double Egdp = (40.3 / std::pow(Afrag,0.2) )*MeV;
+  G4double Egdp   = (40.3 / fG4pow->powZ(Afrag,0.2) )*MeV;
   G4double GammaR = 0.30 * Egdp;
  
-  static G4double normC = 1.0 / ((pi * hbarc)*(pi * hbarc));
-
   // CD
   //cout<<"  PROB TESTS "<<G4endl;
@@ -124 +116 @@
   //cout<<" normC = "<<normC<<G4endl;
 
-  G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR;
-  G4double denominator = (gammaE*gammaE - Egdp*Egdp)*
-           (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE;
-
-  G4double sigmaAbs = numerator/denominator; 
-
-  theProb = normC * sigmaAbs * gammaE*gammaE *
-            levelDensAft/levelDensBef;
+  // VI implementation 18.05.2010
+  G4double gammaE2 = gammaE*gammaE;
+  G4double gammaR2 = gammaE2*GammaR*GammaR;
+  G4double egdp2   = gammaE2 - Egdp*Egdp;
+  G4double sigmaAbs = sigma0*gammaR2/(egdp2*egdp2 + gammaR2); 
+  theProb = normC * sigmaAbs * gammaE2 * levelDens;
+
+  // old implementation
+  //  G4double numerator = sigma0 * gammaE*gammaE * GammaR*GammaR;
+  // G4double denominator = (gammaE*gammaE - Egdp*Egdp)*
+  //         (gammaE*gammaE - Egdp*Egdp) + GammaR*GammaR*gammaE*gammaE;
+
+  //G4double sigmaAbs = numerator/denominator; 
+  //theProb = normC * sigmaAbs * gammaE2 * levelDensAft/levelDensBef;
 
   // CD
@@ -211 +209 @@
 }
 
-G4E1Probability::~G4E1Probability() {}
-
-
+