Changeset 1347 for trunk/source/processes/hadronic/models/de_excitation/photon_evaporation/src/G4E1Probability.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/G4E1Probability.cc (modified) (6 diffs)

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

@@ -24 +24 @@
 // ********************************************************************
 //
-//
-// $Id: G4E1Probability.cc,v 1.9 2010/05/25 10:47:24 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-ref-09 $
+// $Id: G4E1Probability.cc,v 1.11 2010/11/23 18:05:07 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-ref-00 $
 //
 //---------------------------------------------------------------------
@@ -37 +36 @@
 // 18.05.2010 V.Ivanchenko trying to speedup the most slow method
 //            by usage of G4Pow, integer A and introduction of const members
-//
+// 17.11.2010 V.Ivanchenko perform general cleanup and simplification
+//            of integration method; low-limit of integration is defined
+//            by gamma energy or is zero (was always zero before)
 //
 
 #include "G4E1Probability.hh"
-//#include "G4ConstantLevelDensityParameter.hh"
 #include "Randomize.hh"
 #include "G4Pow.hh"
@@ -63 +63 @@
 
 G4double G4E1Probability::EmissionProbDensity(const G4Fragment& frag, 
-                                              const G4double gammaE)
+                                              G4double gammaE)
 {
 
@@ -77 +77 @@
   G4int Afrag = frag.GetA_asInt();
   G4double Uexcite = frag.GetExcitationEnergy();
+  G4double U = std::max(0.0,Uexcite-gammaE);
 
-  if( (Uexcite-gammaE) < 0.0 || gammaE < 0) { return theProb; }
+  if(gammaE < 0.0) { return theProb; }
 
   // Need a level density parameter.
@@ -89 +90 @@
   //  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))); 
-  }
+  G4double levelDens = 
+    std::exp(2*(std::sqrt(aLevelDensityParam*U)-std::sqrt(aLevelDensityParam*Uexcite)));
   // Now form the probability density
 
@@ -140 +138 @@
 
 G4double G4E1Probability::EmissionProbability(const G4Fragment& frag, 
-                                                 const G4double gammaE)
+                                              G4double gammaE)
 {
-
   // From nuclear fragment properties and the excitation energy, calculate
   // the probability for photon evaporation down to last ground level.
   // fragment = nuclear fragment BEFORE de-excitation
 
-  G4double theProb = 0.0;
+  G4double upperLim = gammaE;
+  G4double lowerLim = 0.0; 
 
-  G4double Uafter = 0.0;
-  const G4double Uexcite = frag.GetExcitationEnergy();
-
-  G4double normC = 3.0;
-
-  const G4double upperLim = Uexcite;
-  const G4double lowerLim = Uafter;
-  const G4int numIters = 100;
-
-  // Fall-back is a uniform random number
-
-  //G4double uniformNum = G4UniformRand();
-  //theProb = uniformNum;
+  //G4cout << "G4E1Probability::EmissionProbability:  Emin= " << lowerLim
+  //     << " Emax= " << upperLim << G4endl;
+  if( upperLim - lowerLim <= CLHEP::keV ) { return 0.0; } 
 
   // Need to integrate EmissionProbDensity from lowerLim to upperLim 
-  // and multiply by normC
+  // and multiply by factor 3 (?!)
 
-  G4double integ = normC *
-           EmissionIntegration(frag,gammaE,lowerLim,upperLim,numIters);
-  if(integ > 0.0) theProb = integ/(upperLim-lowerLim);
+  G4double integ = 3.0 * EmissionIntegration(frag,lowerLim,upperLim);
 
-  return theProb;
+  return integ;
 
 }
 
 G4double G4E1Probability::EmissionIntegration(const G4Fragment& frag, 
-                             const G4double ,
-                             const G4double lowLim, const G4double upLim,
-                             const G4int numIters)
+                                              G4double lowLim, G4double upLim)
 
 {
+  // Simple integration
+  // VI replace by direct integration over 100 point
 
-  // Simple Gaussian quadrature integration
+  const G4int numIters = 100;
+  G4double Step = (upLim-lowLim)/G4double(numIters);
 
-  G4double x;
-  G4double root3 = 1.0/std::sqrt(3.0);
+  G4double res = 0.0;
+  G4double x = lowLim - 0.5*Step;
 
-  G4double Step = (upLim-lowLim)/(2.0*numIters);
-  G4double Delta = Step*root3;
-
-  G4double mean = 0.0;
-
-  G4double theInt = 0.0;
-
-  for(G4int i = 0; i < numIters; i++) {
-
-    x = (2*i + 1)*Step;
-    G4double E1ProbDensityA = EmissionProbDensity(frag,x+Delta);
-    G4double E1ProbDensityB = EmissionProbDensity(frag,x-Delta);
-
-    mean += E1ProbDensityA + E1ProbDensityB;
-
+  for(G4int i = 0; i < numIters; ++i) {
+    x += Step;
+    res += EmissionProbDensity(frag, x);
   }
 
-  if(mean*Step > 0.0) theInt = mean*Step;
+  if(res > 0.0) { res /= G4double(numIters); }
+  else { res = 0.0; }
 
-  return theInt;
+  return res;
 
 }