Changeset 1228 for trunk/source/processes/hadronic/models/de_excitation/ablation

Timestamp:

Jan 8, 2010, 11:56:51 AM (15 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

Location:

trunk/source/processes/hadronic/models/de_excitation/ablation

Files:

• include/G4WilsonAblationModel.hh (modified) (4 diffs)
• src/G4WilsonAblationModel.cc (modified) (19 diffs)

trunk/source/processes/hadronic/models/de_excitation/ablation/include/G4WilsonAblationModel.hh

@@ -40 +40 @@
 // MODULE:              G4WilsonAblationModel.hh
 //
-// Version:             B.1
-// Date:                15/04/04
+// Version:             1.0
+// Date:                08/12/2009
 // Author:              P R Truscott
 // Organisation:        QinetiQ Ltd, UK
@@ -58 +58 @@
 // Beta release
 //
+// 08 December 2009, P R Truscott, QinetiQ Ltd, UK
+// Ver 1.0
+// Introduced vector of evaporation channels and evaporation factory.  Also
+// copied directly over the SumProbabilities class in G4Evaporation.hh at
+// version 9.2.r9, just in cases there's any subtle differences.  See .cc
+// file comments to see impact of the rest of the changes.
+//
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ////////////////////////////////////////////////////////////////////////////////
@@ -67 +74 @@
 #include "G4ParticleDefinition.hh"
 #include "globals.hh"
+#include "G4VEvaporationFactory.hh"
+#include "G4EvaporationFactory.hh"
+
 
 #include <vector>
@@ -99 +109 @@
     VectorOfFragmentTypes  evapType;
 
-    class SumProbabilities : 
-      public std::binary_function<G4double,G4double,G4double>
-    {
-      public:
-        SumProbabilities() : total(0.0) {}
-        G4double operator() (G4double& /* probSoFar */, G4VEvaporationChannel*& frag)
-        {
-          total += frag->GetEmissionProbability();
-          return total;
-        }
+    std::vector<G4VEvaporationChannel*> * theChannels;
+    G4VEvaporationFactory * theChannelFactory;
 
-        G4double GetTotal() { return total; }
-      public:
-      G4double total;
+  class SumProbabilities : public std::binary_function<G4double,G4double,G4double>
+  {
+  public:
+    SumProbabilities() : total(0.0) {}
+    G4double operator() (G4double& /* probSoFar */, G4VEvaporationChannel*& frag)
+    { 
+      G4double temp_prob = frag->GetEmissionProbability();
+      if(temp_prob >= 0.0) total += temp_prob;
+      //      total += frag->GetEmissionProbability();
+      return total;
+    }
+    
+    G4double GetTotal() { return total; }
+  public:
+    G4double total;
+    
   };
-                                                                                
-
 
 };

trunk/source/processes/hadronic/models/de_excitation/ablation/src/G4WilsonAblationModel.cc

@@ -38 +38 @@
 // MODULE:              G4WilsonAblationModel.cc
 //
-// Version:             B.1
-// Date:                15/04/04
+// Version:             1.0
+// Date:                08/12/2009
 // Author:              P R Truscott
 // Organisation:        QinetiQ Ltd, UK
@@ -56 +56 @@
 // Beta release
 //
+// 08 December 2009, P R Truscott, QinetiQ Ltd, UK
+// Ver 1.0
+// Updated as a result of changes in the G4Evaporation classes.  These changes
+// affect mostly SelectSecondariesByEvaporation, and now you have variables
+// associated with the evaporation model which can be changed:
+//    OPTxs to select the inverse cross-section
+//    OPTxs = 0      => Dostrovski's parameterization
+//    OPTxs = 1 or 2 => Chatterjee's paramaterization
+//    OPTxs = 3 or 4 => Kalbach's parameterization
+//    useSICB        => use superimposed Coulomb Barrier for inverse cross
+//                      sections
+// Other problem found with G4Fragment definition using Lorentz vector and
+// **G4ParticleDefinition**.  This does not allow A and Z to be defined for the
+// fragment for some reason.  Now the fragment is defined more explicitly:
+//    G4Fragment *fragment = new G4Fragment(A, Z, lorentzVector);
+// to avoid this quirk.  Bug found in SelectSecondariesByDefault: *type is now
+// equated to evapType[i] whereas previously it was equated to fragType[i].
+// 
 // %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ////////////////////////////////////////////////////////////////////////////////
@@ -123 +141 @@
 //
   verboseLevel = 0;
+  theChannelFactory = new G4EvaporationFactory();
+  theChannels = theChannelFactory->GetChannel();
+//
+//
+// Set defaults for evaporation classes.  These can be overridden by user
+// "set" methods.
+//
+  OPTxs   = 3;
+  useSICB = false;
 }
 ////////////////////////////////////////////////////////////////////////////////
 //
 G4WilsonAblationModel::~G4WilsonAblationModel()
-{;}
+{
+  if (theChannels != 0) theChannels = 0;
+  if (theChannelFactory != 0) delete theChannelFactory;
+}
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -207 +237 @@
   G4int DAabl = (G4int) (ex / B);
   if (DAabl > A) DAabl = A;
-  if (verboseLevel >= 2)
-    G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;
+// The following lines are no longer accurate given we now treat the final fragment
+//  if (verboseLevel >= 2)
+//    G4cout <<"Number of nucleons ejected = " <<DAabl <<G4endl;
 
 //
@@ -219 +250 @@
   if (AF > 0)
   {
-    G4double AFd = static_cast<G4double>(AF);
+    G4double AFd = (G4double) AF;
     G4double R = 11.8 / std::pow(AFd, 0.45);
     G4int minZ = Z - DAabl;
@@ -253 +284 @@
   }
   G4int DZabl = Z - ZF;
-  if (verboseLevel >= 2)
-    G4cout <<"Final fragment      A=" <<AF
-           <<", Z=" <<ZF
-           <<G4endl;
 //
 //
@@ -285 +312 @@
       DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
       DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
-      if (verboseLevel >= 2)
-        G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
-               <<", number of particles emitted = " <<n
-               <<G4endl;
-    }
-  }
-//
-//
-// Determine the properties of the final nuclear fragment.
+    }
+  }
+//
+//
+// Determine the properties of the final nuclear fragment.  Note that if
+// the final fragment is predicted to have a nucleon number of zero, then
+// really it's the particle last in the vector evapType which becomes the
+// final fragment.  Therefore delete this from the vector if this is the
+// case.
 //
   G4double massFinalFrag = 0.0;
-  if (AF > 0.0)
+  if (AF > 0)
     massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()->
       GetIonMass(ZF,AF);
+  else
+  {
+    G4ParticleDefinition *type = evapType[evapType.size()-1];
+    AF                         = type->GetBaryonNumber();
+    ZF                         = (G4int) (type->GetPDGCharge() + 1.0E-10);
+    evapType.erase(evapType.end()-1);
+  }
   totalEpost   += massFinalFrag;
 //
+//
+// Provide verbose output on the nuclear fragment if requested.
+//
+  if (verboseLevel >= 2)
+  {
+    G4cout <<"Final fragment      A=" <<AF
+           <<", Z=" <<ZF
+           <<G4endl;
+    for (G4int ift=0; ift<nFragTypes; ift++)
+    {
+      G4ParticleDefinition *type = fragType[ift];
+      G4int n                    = std::count(evapType.begin(),evapType.end(),type);
+      if (n > 0) 
+        G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
+               <<", number of particles emitted = " <<n <<G4endl;
+    }
+  }
 //
 // Add the total energy from the fragment.  Note that the fragment is assumed
@@ -326 +377 @@
       SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
   }
+
   if (AF > 0)
   {
@@ -352 +404 @@
     G4int ie = 0;
     G4FragmentVector::iterator iter;
-    for (iter = fragmentVector->begin(); iter != fragmentVector->end(); ++iter)
+    for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++)
     {
       if (ie == nEvap)
       {
-        G4cout <<*iter  <<G4endl;
+//        G4cout <<*iter  <<G4endl;
         G4cout <<"---------------------------------" <<G4endl;
         G4cout <<"Particles from default emission :" <<G4endl;
@@ -375 +427 @@
   (G4Fragment *intermediateNucleus)
 {
+  G4Fragment theResidualNucleus = *intermediateNucleus;
   G4bool evaporate = true;
   while (evaporate && evapType.size() != 0)
@@ -386 +439 @@
     std::vector <G4VEvaporationChannel*>  theChannels;
     theChannels.clear();
+    std::vector <G4VEvaporationChannel*>::iterator i;
     VectorOfFragmentTypes::iterator iter;
     std::vector <VectorOfFragmentTypes::iterator> iters;
@@ -393 +447 @@
     {
       theChannels.push_back(new G4AlphaEvaporationChannel);
+      i = theChannels.end() - 1;
+      (*i)->SetOPTxs(OPTxs);
+      (*i)->UseSICB(useSICB);
+//      (*i)->Initialize(theResidualNucleus);
       iters.push_back(iter);
     }
@@ -399 +457 @@
     {
       theChannels.push_back(new G4He3EvaporationChannel);
+      i = theChannels.end() - 1;
+      (*i)->SetOPTxs(OPTxs);
+      (*i)->UseSICB(useSICB);
+//      (*i)->Initialize(theResidualNucleus);
       iters.push_back(iter);
     }
@@ -405 +467 @@
     {
       theChannels.push_back(new G4TritonEvaporationChannel);
+      i = theChannels.end() - 1;
+      (*i)->SetOPTxs(OPTxs);
+      (*i)->UseSICB(useSICB);
+//      (*i)->Initialize(theResidualNucleus);
       iters.push_back(iter);
     }
@@ -411 +477 @@
     {
       theChannels.push_back(new G4DeuteronEvaporationChannel);
+      i = theChannels.end() - 1;
+      (*i)->SetOPTxs(OPTxs);
+      (*i)->UseSICB(useSICB);
+//      (*i)->Initialize(theResidualNucleus);
       iters.push_back(iter);
     }
@@ -417 +487 @@
     {
       theChannels.push_back(new G4ProtonEvaporationChannel);
+      i = theChannels.end() - 1;
+      (*i)->SetOPTxs(OPTxs);
+      (*i)->UseSICB(useSICB);
+//      (*i)->Initialize(theResidualNucleus);
       iters.push_back(iter);
     }
@@ -423 +497 @@
     {
       theChannels.push_back(new G4NeutronEvaporationChannel);
+      i = theChannels.end() - 1;
+      (*i)->SetOPTxs(OPTxs);
+      (*i)->UseSICB(useSICB);
+//      (*i)->Initialize(theResidualNucleus);
       iters.push_back(iter);
     }
@@ -479 +557 @@
   for (unsigned i=0; i<evapType.size(); i++)
   {
-    G4ParticleDefinition *type = fragType[i];
+    G4ParticleDefinition *type = evapType[i];
     G4double mass              = type->GetPDGMass();
     G4double e                 = mass + 10.0*eV;
@@ -489 +567 @@
     G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
     lorentzVector.boost(-boost);
+// Possibility that the following line is not correctly carrying over A and Z
+// from particle definition.  Force values.  PRT 03/12/2009.
+//    G4Fragment *fragment          = 
+//      new G4Fragment(lorentzVector, type);
+    G4int A = type->GetBaryonNumber();
+    G4int Z = (G4int) (type->GetPDGCharge() + 1.0E-10);
     G4Fragment *fragment          = 
-      new G4Fragment(lorentzVector, type);
+      new G4Fragment(A, Z, lorentzVector);
+
     fragmentVector->push_back(fragment);
   }