Changeset 1347 for trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc

Timestamp:

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

Author:

garnier

Message:

geant4 tag 9.4

File:

• trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc (modified) (15 diffs)

trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationChannel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
+// $Id: G4EvaporationChannel.cc,v 1.19 2010/11/24 15:30:49 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-04-ref-00 $
 //
 //J.M. Quesada (August2008). Based on:
@@ -30 +32 @@
 // by V. Lara (Oct 1998)
 //
-// Modif (03 September 2008) by J. M. Quesada for external choice of inverse 
-// cross section option
-// JMQ (06 September 2008) Also external choices have been added for 
-// superimposed Coulomb barrier (if useSICB is set true, by default is false) 
-
+// Modified:
+// 03-09-2008 J.M. Quesada for external choice of inverse cross section option
+// 06-09-2008 J.M. Quesada Also external choices have been added for superimposed 
+//                 Coulomb barrier (if useSICB is set true, by default is false) 
+// 17-11-2010 V.Ivanchenko in constructor replace G4VEmissionProbability by 
+//            G4EvaporationProbability and do not new and delete probability
+//            object at each call; use G4Pow
 
 #include "G4EvaporationChannel.hh"
 #include "G4PairingCorrection.hh"
-
-
-
-G4EvaporationChannel::G4EvaporationChannel(const G4int anA, const G4int aZ, const G4String & aName,
-                                           G4VEmissionProbability * aEmissionStrategy,
-                                           G4VCoulombBarrier * aCoulombBarrier):
+#include "G4NucleiProperties.hh"
+#include "G4Pow.hh"
+#include "G4EvaporationLevelDensityParameter.hh"
+#include "Randomize.hh"
+#include "G4Alpha.hh"
+
+G4EvaporationChannel::G4EvaporationChannel(G4int anA, G4int aZ, 
+                                           const G4String & aName,
+                                           G4EvaporationProbability* aEmissionStrategy,
+                                           G4VCoulombBarrier* aCoulombBarrier):
     G4VEvaporationChannel(aName),
     theA(anA),
@@ -52 +60 @@
     MaximalKineticEnergy(-1000.0)
 { 
-    theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
-//    MyOwnLevelDensity = true;  
+  ResidualA = 0;
+  ResidualZ = 0;
+  ResidualMass = CoulombBarrier=0.0;
+  EvaporatedMass = G4NucleiProperties::GetNuclearMass(theA, theZ);
+  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
+}
+
+G4EvaporationChannel::G4EvaporationChannel():
+    G4VEvaporationChannel(""),
+    theA(0),
+    theZ(0),
+    theEvaporationProbabilityPtr(0),
+    theCoulombBarrierPtr(0),
+    EmissionProbability(0.0),
+    MaximalKineticEnergy(-1000.0)
+{ 
+  ResidualA = 0;
+  ResidualZ = 0;
+  EvaporatedMass = ResidualMass = CoulombBarrier = 0.0;
+  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
 }
 
@@ -60 +86 @@
   delete theLevelDensityPtr;
 }
-
-G4EvaporationChannel::G4EvaporationChannel(const G4EvaporationChannel & ) : G4VEvaporationChannel()
-{
-    throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::copy_costructor meant to not be accessable");
-}
-
-const G4EvaporationChannel & G4EvaporationChannel::operator=(const G4EvaporationChannel & )
-{
-    throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationChannel::operator= meant to not be accessable");
-    return *this;
-}
-
-G4bool G4EvaporationChannel::operator==(const G4EvaporationChannel & right) const 
-{
-    return (this == (G4EvaporationChannel *) &right);
-    //  return false;
-}
-
-G4bool G4EvaporationChannel::operator!=(const G4EvaporationChannel & right) const 
-{
-    return (this != (G4EvaporationChannel *) &right);
-    //  return true;
-}
-
-//void G4EvaporationChannel::SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
-//  {
-//    if (MyOwnLevelDensity) delete theLevelDensityPtr;
-//    theLevelDensityPtr = aLevelDensity;
-//    MyOwnLevelDensity = false;
-//  }
 
 void G4EvaporationChannel::Initialize(const G4Fragment & fragment)
@@ -98 +94 @@
   theEvaporationProbabilityPtr->UseSICB(useSICB);
   
-  
-  G4int FragmentA = static_cast<G4int>(fragment.GetA());
-  G4int FragmentZ = static_cast<G4int>(fragment.GetZ());
+  G4int FragmentA = fragment.GetA_asInt();
+  G4int FragmentZ = fragment.GetZ_asInt();
   ResidualA = FragmentA - theA;
   ResidualZ = FragmentZ - theZ;
+  //G4cout << "G4EvaporationChannel::Initialize Z= " << theZ << " A= " << theA 
+  //     << " FragZ= " << FragmentZ << " FragA= " << FragmentA << G4endl;
   
   //Effective excitation energy
@@ -108 +105 @@
     G4PairingCorrection::GetInstance()->GetPairingCorrection(FragmentA,FragmentZ);
   
-  
   // Only channels which are physically allowed are taken into account 
   if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
       (ResidualA == ResidualZ && ResidualA > 1) || ExEnergy <= 0.0) {
-    CoulombBarrier=0.0;
+    CoulombBarrier = ResidualMass = 0.0;
     MaximalKineticEnergy = -1000.0*MeV;
     EmissionProbability = 0.0;
   } else {
+    ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
+    G4double FragmentMass = fragment.GetGroundStateMass();
     CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
     // Maximal Kinetic Energy
-    MaximalKineticEnergy = CalcMaximalKineticEnergy
-      (G4ParticleTable::GetParticleTable()->
-       GetIonTable()->GetNucleusMass(FragmentZ,FragmentA)+ExEnergy);
+    MaximalKineticEnergy = CalcMaximalKineticEnergy(FragmentMass + ExEnergy);
+    //MaximalKineticEnergy = ExEnergy + fragment.GetGroundStateMass() 
+    //  - EvaporatedMass - ResidualMass;
     
     // Emission probability
@@ -131 +129 @@
       limit= CoulombBarrier;
     else limit=0.;
+    limit = std::max(limit, FragmentMass - ResidualMass - EvaporatedMass);
   
     // The threshold for charged particle emission must be  set to 0 if Coulomb 
     //cutoff  is included in the cross sections
-    //if (MaximalKineticEnergy <= 0.0) EmissionProbability = 0.0;  
     if (MaximalKineticEnergy <= limit) EmissionProbability = 0.0;
     else { 
@@ -142 +140 @@
     }
   }
+  //G4cout << "G4EvaporationChannel:: probability= " << EmissionProbability << G4endl; 
   
   return;
@@ -148 +147 @@
 G4FragmentVector * G4EvaporationChannel::BreakUp(const G4Fragment & theNucleus)
 {
-  G4double EvaporatedKineticEnergy=GetKineticEnergy(theNucleus);
-  
-  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
-    GetIonMass(theZ,theA);
-  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
-  
+  /*
+  G4double Ecm = GetKineticEnergy(theNucleus) + ResidualMass + EvaporatedMass;
+  
+  G4double EvaporatedEnergy = 
+    ((Ecm-ResidualMass)*(Ecm+ResidualMass) + EvaporatedMass*EvaporatedMass)/(2*Ecm);
+  */
+  G4double EvaporatedEnergy = GetKineticEnergy(theNucleus) + EvaporatedMass;
+
   G4ThreeVector momentum(IsotropicVector
-                         (std::sqrt(EvaporatedKineticEnergy*
-                                    (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
-  
-  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
+                         (std::sqrt((EvaporatedEnergy - EvaporatedMass)*
+                                    (EvaporatedEnergy + EvaporatedMass))));
   
   G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
-  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
+  G4LorentzVector ResidualMomentum = theNucleus.GetMomentum();
+  EvaporatedMomentum.boost(ResidualMomentum.boostVector());
   
   G4Fragment * EvaporatedFragment = new G4Fragment(theA,theZ,EvaporatedMomentum);
-#ifdef PRECOMPOUND_TEST
-  EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
-#endif
-  // ** And now the residual nucleus ** 
-  G4double theExEnergy = theNucleus.GetExcitationEnergy();
-  G4double theMass = G4ParticleTable::GetParticleTable()->GetIonTable()->
-    GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
-                   static_cast<G4int>(theNucleus.GetA()));
-  G4double ResidualEnergy = theMass + 
-    (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
-  
-  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
-  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
-  
-  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum );
-  
-#ifdef PRECOMPOUND_TEST
-  ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
-#endif
+  ResidualMomentum -= EvaporatedMomentum;
+
+  G4Fragment * ResidualFragment = new G4Fragment(ResidualA, ResidualZ, ResidualMomentum);
+ 
   G4FragmentVector * theResult = new G4FragmentVector;
   
@@ -211 +196 @@
 /////////////////////////////////////////
 // Calculates the maximal kinetic energy that can be carried by fragment.
-G4double G4EvaporationChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
-{
-  G4double ResidualMass = G4ParticleTable::GetParticleTable()->
-    GetIonTable()->GetNucleusMass( ResidualZ, ResidualA );
-  G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
-    GetIonTable()->GetNucleusMass( theZ, theA );
-
+G4double G4EvaporationChannel::CalcMaximalKineticEnergy(G4double NucleusTotalE)
+{
   // This is the "true" assimptotic kinetic energy (from energy conservation)   
-  G4double Tmax = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass -                
-                   ResidualMass*ResidualMass)/(2.0*NucleusTotalE) - EvaporatedMass;
+  G4double Tmax = 
+    ((NucleusTotalE-ResidualMass)*(NucleusTotalE+ResidualMass) + EvaporatedMass*EvaporatedMass)
+    /(2.0*NucleusTotalE) - EvaporatedMass;
   
   //JMQ (13-09-08) bug fixed: in the original version the Tmax is calculated
@@ -245 +226 @@
     
     
-    if (MaximalKineticEnergy < 0.0)   
+    if (MaximalKineticEnergy < 0.0) {
       throw G4HadronicException(__FILE__, __LINE__, 
                                 "G4EvaporationChannel::CalcKineticEnergy: maximal kinetic at the Coulomb barrier is less than 0");
-    
+    }
     G4double Rb = 4.0*theLevelDensityPtr->
       LevelDensityParameter(ResidualA+theA,ResidualZ+theZ,MaximalKineticEnergy)*
@@ -263 +244 @@
       G4double Q2 = 1.0;
       if (theZ == 0) { // for emitted neutron
-        G4double Beta = (2.12/std::pow(G4double(ResidualA),2./3.) - 0.05)*MeV/
-          (0.76 + 2.2/std::pow(G4double(ResidualA),1./3.));
+        G4double Beta = (2.12/G4Pow::GetInstance()->Z23(ResidualA) - 0.05)*MeV/
+          (0.76 + 2.2/G4Pow::GetInstance()->Z13(ResidualA));
         Q1 = 1.0 + Beta/(MaximalKineticEnergy);
         Q2 = Q1*std::sqrt(Q1);
@@ -277 +258 @@
   } else if (OPTxs==1 || OPTxs==2 || OPTxs==3 || OPTxs==4) {
     
-    G4double V;
-    if(useSICB) V= CoulombBarrier;
-    else V=0.;
-    //If Coulomb barrier is just included  in the cross sections
-    //  G4double V=0.;
+    // Coulomb barrier is just included  in the cross sections
+    G4double V = 0;
+    if(useSICB) { V= CoulombBarrier; }
+
+    V = std::max(V, aFragment.GetGroundStateMass()-EvaporatedMass-ResidualMass);
 
     G4double Tmax=MaximalKineticEnergy;
     G4double T(0.0);
     G4double NormalizedProbability(1.0);
-    
+
+    // VI: This is very ineffective - create new objects at each call to the method    
+    /*
     // A pointer is created in order to access the distribution function.
-    G4EvaporationProbability * G4EPtemp;
+    G4EvaporationProbability * G4EPtemp = 0;
     
     if (theA==1 && theZ==0) G4EPtemp=new G4NeutronEvaporationProbability();
@@ -305 +288 @@
       G4EPtemp->SetOPTxs(OPTxs);
       G4EPtemp->UseSICB(useSICB);
-
+    */
+
+    // use local pointer and not create a new one
     do
       {  
         T=V+G4UniformRand()*(Tmax-V);
-        NormalizedProbability=G4EPtemp->ProbabilityDistributionFunction(aFragment,T)/
-          (this->GetEmissionProbability());
+        NormalizedProbability = 
+          theEvaporationProbabilityPtr->ProbabilityDistributionFunction(aFragment,T)/
+          GetEmissionProbability();
         
       }
     while (G4UniformRand() > NormalizedProbability);
-   delete G4EPtemp;
-   return T;
+    //   delete G4EPtemp;
+    return T;
   } else{
     std::ostringstream errOs;
@@ -323 +309 @@
 }
 
-G4ThreeVector G4EvaporationChannel::IsotropicVector(const G4double Magnitude)
+G4ThreeVector G4EvaporationChannel::IsotropicVector(G4double Magnitude)
     // Samples a isotropic random vectorwith a magnitud given by Magnitude.
     // By default Magnitude = 1.0
 {
-    G4double CosTheta = 1.0 - 2.0*G4UniformRand();
-    G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
-    G4double Phi = twopi*G4UniformRand();
-    G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
-                         Magnitude*std::sin(Phi)*SinTheta,
-                         Magnitude*CosTheta);
-    return Vector;
-            }
+  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
+  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
+  G4double Phi = twopi*G4UniformRand();
+  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
+                       Magnitude*std::sin(Phi)*SinTheta,
+                       Magnitude*CosTheta);
+  return Vector;
+}