Ignore:
Timestamp:
Dec 22, 2010, 3:52:27 PM (14 years ago)
Author:
garnier
Message:

geant4 tag 9.4

Location:
trunk/source/processes/hadronic/models/de_excitation/gem_evaporation
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4AlphaGEMChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4AlphaGEMChannel.hh,v 1.4 2009/09/15 12:54:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4AlphaGEMChannel.hh,v 1.5 2010/10/29 17:35:04 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    5757private:
    5858    const G4AlphaGEMChannel & operator=(const G4AlphaGEMChannel & right); 
    59    
    60     G4AlphaGEMChannel(const G4AlphaGEMChannel & right);
    61    
    62 public:
     59    G4AlphaGEMChannel(const G4AlphaGEMChannel & right);   
    6360    G4bool operator==(const G4AlphaGEMChannel & right) const;
    6461    G4bool operator!=(const G4AlphaGEMChannel & right) const;
    6562   
    66 private:
    6763// JMQ 190709
    6864//       G4AlphaCoulombBarrier theCoulombBarrier;

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMChannel.hh

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4GEMChannel.hh,v 1.5 2009/09/15 12:54:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4GEMChannel.hh,v 1.7 2010/11/18 16:21:17 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4848//#define debug
    4949
     50class G4Pow;
     51
    5052class G4GEMChannel : public G4VEvaporationChannel
    5153{
    5254public:
    53     // Available constructors
    54     G4GEMChannel(const G4int theA, const G4int theZ,
    55                  G4GEMProbability * aEmissionStrategy,
    56                  G4VCoulombBarrier * aCoulombBarrier);
     55
     56  G4GEMChannel(const G4int theA, const G4int theZ, const G4String & aName,
     57               G4GEMProbability * aEmissionStrategy,
     58               G4VCoulombBarrier * aCoulombBarrier);
     59
     60  // destructor
     61  virtual ~G4GEMChannel();
    5762   
    58     G4GEMChannel(const G4int theA, const G4int theZ, const G4String & aName,
    59                  G4GEMProbability * aEmissionStrategy,
    60                  G4VCoulombBarrier * aCoulombBarrier);
     63  void Initialize(const G4Fragment & fragment);
     64
     65  G4FragmentVector * BreakUp(const G4Fragment & theNucleus);
     66
     67  inline void SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
     68  {
     69    if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
     70    theLevelDensityPtr = aLevelDensity;
     71    MyOwnLevelDensity = false;
     72  }
     73
     74  inline G4double GetEmissionProbability(void) const
     75  { return EmissionProbability; }
     76 
     77  inline G4double GetMaximalKineticEnergy(void) const
     78  { return MaximalKineticEnergy; }
     79 
     80private:
    6181   
    62     G4GEMChannel(const G4int theA, const G4int theZ, const G4String * aName,
    63                  G4GEMProbability * aEmissionStrategy,
    64                  G4VCoulombBarrier * aCoulombBarrier);
    65    
    66 public:
    67     // destructor
    68     ~G4GEMChannel();
    69  
    70     void SetEmissionStrategy(G4GEMProbability * aEmissionStrategy)
    71         {
    72             theEvaporationProbabilityPtr = aEmissionStrategy;
    73         }
    74  
    75     void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
    76         {
    77             theCoulombBarrierPtr = aCoulombBarrier;
    78         }
     82  // Calculate Binding Energy for separate fragment from nucleus
     83  G4double CalcBindingEnergy(G4int anA, G4int aZ);
     84
     85  // Calculate maximal kinetic energy that can be carried by fragment (in MeV)
     86  G4double CalcMaximalKineticEnergy(G4double U);
     87
     88  // Samples fragment kinetic energy.
     89  G4double CalcKineticEnergy(const G4Fragment & fragment);
     90
     91  // This has to be removed and put in Random Generator
     92  G4ThreeVector IsotropicVector(G4double Magnitude  = 1.0);
     93
     94  G4GEMChannel(const G4GEMChannel & right); 
     95  const G4GEMChannel & operator=(const G4GEMChannel & right);
     96  G4bool operator==(const G4GEMChannel & right) const;
     97  G4bool operator!=(const G4GEMChannel & right) const;
    7998 
    8099protected:
    81     // default constructor
    82     G4GEMChannel() {};
    83  
     100  G4GEMChannel();
     101
     102  // Data Members ************
     103
    84104private:
    85     // copy constructor
    86     G4GEMChannel(const G4GEMChannel & right);
    87  
    88 private:
    89     const G4GEMChannel & operator=(const G4GEMChannel & right);
    90  
    91 public:
    92     G4bool operator==(const G4GEMChannel & right) const;
    93     G4bool operator!=(const G4GEMChannel & right) const;
    94105
    95 public:
    96     void Initialize(const G4Fragment & fragment);
     106  // This data member define the channel.
     107  // They are intializated at object creation (constructor) time.
     108   
     109  // Atomic Number
     110  G4int A;
     111   
     112  // Charge
     113  G4int Z;
    97114
    98     G4FragmentVector * BreakUp(const G4Fragment & theNucleus);
     115  G4double EvaporatedMass;
     116  G4double ResidualMass;
    99117
    100     inline void SetLevelDensityParameter(G4VLevelDensityParameter * aLevelDensity)
    101         {
    102             if (MyOwnLevelDensity) delete theLevelDensityPtr;
    103             theLevelDensityPtr = aLevelDensity;
    104             MyOwnLevelDensity = false;
    105         }
    106  
    107 public:
     118  G4Pow* fG4pow;
     119   
     120  // For evaporation probability calcualtion
     121  G4GEMProbability * theEvaporationProbabilityPtr;
     122   
     123  // For Level Density calculation
     124  G4bool MyOwnLevelDensity;
     125  G4VLevelDensityParameter * theLevelDensityPtr;
     126   
     127  // For Coulomb Barrier calculation
     128  G4VCoulombBarrier * theCoulombBarrierPtr;
     129  G4double CoulombBarrier;
     130   
     131  //---------------------------------------------------
     132   
     133  // These values depend on the nucleus that is being evaporated.
     134  // They are calculated through the Initialize method which takes as parameters
     135  // the atomic number, charge and excitation energy of nucleus.
     136   
     137  // Residual Atomic Number
     138  G4int ResidualA;
     139
     140  // Residual Charge
     141  G4int ResidualZ;
     142   
     143  // Emission Probability
     144  G4double EmissionProbability;
     145
     146  // Maximal Kinetic Energy that can be carried by fragment
     147  G4double MaximalKineticEnergy;
    108148
    109149
    110     inline G4double GetEmissionProbability(void) const
    111         {
    112             return EmissionProbability;
    113         }
    114  
    115  
    116     inline G4double GetMaximalKineticEnergy(void) const
    117         {
    118             return MaximalKineticEnergy;
    119         }
    120  
    121     // ----------------------
    122    
    123 private:
    124    
    125     // Calculate Binding Energy for separate fragment from nucleus
    126     G4double CalcBindingEnergy(const G4int anA, const G4int aZ);
    127 
    128     // Calculate maximal kinetic energy that can be carried by fragment (in MeV)
    129     G4double CalcMaximalKineticEnergy(const G4double U);
    130 
    131     // Samples fragment kinetic energy.
    132     G4double CalcKineticEnergy(const G4Fragment & fragment);
    133 
    134     // This has to be removed and put in Random Generator
    135     G4ThreeVector IsotropicVector(const G4double Magnitude  = 1.0);
    136 
    137         // Data Members
    138         // ************
    139 private:
    140 
    141     // This data member define the channel.
    142   // They are intializated at object creation (constructor) time.
    143    
    144     // Atomic Number
    145     G4int A;
    146    
    147     // Charge
    148     G4int Z;
    149    
    150 
    151     // For evaporation probability calcualtion
    152     G4GEMProbability * theEvaporationProbabilityPtr;
    153    
    154     // For Level Density calculation
    155     G4bool MyOwnLevelDensity;
    156     G4VLevelDensityParameter * theLevelDensityPtr;
    157    
    158     // For Coulomb Barrier calculation
    159     G4VCoulombBarrier * theCoulombBarrierPtr;
    160     G4double CoulombBarrier;
    161    
    162     //---------------------------------------------------
    163    
    164     // These values depend on the nucleus that is being evaporated.
    165     // They are calculated through the Initialize method which takes as parameters
    166     // the atomic number, charge and excitation energy of nucleus.
    167    
    168     // Residual Atomic Number
    169     G4int AResidual;
    170 
    171     // Residual Charge
    172     G4int ZResidual;
    173    
    174     // Emission Probability
    175     G4double EmissionProbability;
    176 
    177 
    178     // Maximal Kinetic Energy that can be carried by fragment
    179     G4double MaximalKineticEnergy;
    180150};
    181151

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/include/G4GEMProbability.hh

    r1340 r1347  
    2424// ********************************************************************
    2525//
    26 // $Id: G4GEMProbability.hh,v 1.4 2010/05/19 10:21:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     26// $Id: G4GEMProbability.hh,v 1.6 2010/11/05 14:42:52 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2828//
    2929//---------------------------------------------------------------------
    … …  
    4343#define G4GEMProbability_h 1
    4444
    45 
    4645#include "G4VEmissionProbability.hh"
    4746#include "G4VLevelDensityParameter.hh"
    … …  
    6059
    6160  // Only available constructor
    62   G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin);
     61  G4GEMProbability(G4int anA, G4int aZ, G4double aSpin);
    6362   
    6463  virtual ~G4GEMProbability();
    … …  
    114113  }
    115114
    116   inline void SetCoulombBarrierStrategy(G4VCoulombBarrier * aCoulombBarrier)
    117   {
    118     theCoulombBarrierPtr = aCoulombBarrier;
    119   }
    120 
    121115private:
    122116
    … …  
    130124public:
    131125
    132   G4double EmissionProbability(const G4Fragment & fragment, const G4double anEnergy);
     126  G4double EmissionProbability(const G4Fragment & fragment, G4double anEnergy);
    133127 
    134128private:
    135129
    136   G4double CalcProbability(const G4Fragment & fragment, const G4double MaximalKineticEnergy,
    137                            const G4double V);
     130  G4double CalcProbability(const G4Fragment & fragment, G4double MaximalKineticEnergy,
     131                           G4double V);
    138132
    139   virtual G4double CCoeficient(const G4double ) const;
     133  virtual G4double CCoeficient(G4double ) const;
    140134
    141   inline G4double I0(const G4double t);
    142   inline G4double I1(const G4double t, const G4double tx);
    143   inline G4double I2(const G4double s, const G4double sx);
    144   G4double I3(const G4double s, const G4double sx);
     135  inline G4double I0(G4double t);
     136  inline G4double I1(G4double t, G4double tx);
     137  inline G4double I2(G4double s, G4double sx);
     138  G4double I3(G4double s, G4double sx);
    145139   
    146140  // Data Members
    … …  
    174168};
    175169
    176 inline G4double G4GEMProbability::I0(const G4double t)
     170inline G4double G4GEMProbability::I0(G4double t)
    177171{
    178172  return std::exp(t) - 1.0;
    179173}
    180174
    181 inline G4double G4GEMProbability::I1(const G4double t, const G4double tx)
     175inline G4double G4GEMProbability::I1(G4double t, G4double tx)
    182176{
    183177  return (t - tx + 1.0)*std::exp(tx) - t - 1.0;
    … …  
    185179
    186180
    187 inline G4double G4GEMProbability::I2(const G4double s, const G4double sx)
     181inline G4double G4GEMProbability::I2(G4double s, G4double sx)
    188182{
    189183  G4double S = 1.0/std::sqrt(s);

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4AlphaGEMChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4AlphaGEMChannel.cc,v 1.5 2009/09/15 12:54:16 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4AlphaGEMChannel.cc,v 1.6 2010/10/29 17:35:04 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3535
    3636
    37 const G4AlphaGEMChannel & G4AlphaGEMChannel::operator=(const G4AlphaGEMChannel & )
    38 {
    39     throw G4HadronicException(__FILE__, __LINE__, "G4AlphaGEMChannel::operator= meant to not be accessable");
    40     return *this;
    41 }
    42 
    43 G4AlphaGEMChannel::G4AlphaGEMChannel(const G4AlphaGEMChannel & ): G4GEMChannel()
    44 {
    45     throw G4HadronicException(__FILE__, __LINE__, "G4AlphaGEMChannel::CopyConstructor meant to not be accessable");
    46 }
    47 
    48 G4bool G4AlphaGEMChannel::operator==(const G4AlphaGEMChannel & right) const
    49 {
    50     return (this == (G4AlphaGEMChannel *) &right);
    51     //  return false;
    52 }
    53 
    54 G4bool G4AlphaGEMChannel::operator!=(const G4AlphaGEMChannel & right) const
    55 {
    56     return (this != (G4AlphaGEMChannel *) &right);
    57     //  return true;
    58 }
    59 

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMChannel.cc

    r1340 r1347  
    2525//
    2626//
    27 // $Id: G4GEMChannel.cc,v 1.10 2009/10/08 07:55:39 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     27// $Id: G4GEMChannel.cc,v 1.12 2010/11/18 16:21:17 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    4040#include "G4GEMChannel.hh"
    4141#include "G4PairingCorrection.hh"
    42 
    43 
    44 G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ,
    45                            G4GEMProbability * aEmissionStrategy,
    46                            G4VCoulombBarrier * aCoulombBarrier):
    47     A(theA),
    48     Z(theZ),
    49     theEvaporationProbabilityPtr(aEmissionStrategy),
    50     theCoulombBarrierPtr(aCoulombBarrier),
    51     EmissionProbability(0.0),
    52     MaximalKineticEnergy(-1000.0)
    53 {
    54     theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    55     MyOwnLevelDensity = true;
    56 }
    57 
    58 G4GEMChannel::G4GEMChannel(const G4int theA, const G4int theZ, const G4String & aName,
     42#include "G4Pow.hh"
     43
     44G4GEMChannel::G4GEMChannel(G4int theA, G4int theZ, const G4String & aName,
    5945                           G4GEMProbability * aEmissionStrategy,
    6046                           G4VCoulombBarrier * aCoulombBarrier) :
    … …  
    6955  theLevelDensityPtr = new G4EvaporationLevelDensityParameter;
    7056  MyOwnLevelDensity = true;
     57  EvaporatedMass = G4NucleiProperties::GetNuclearMass(A, Z);
     58  ResidualMass = CoulombBarrier = 0.0;
     59  fG4pow = G4Pow::GetInstance();
     60}
     61
     62G4GEMChannel::G4GEMChannel() :
     63  G4VEvaporationChannel("dummy"),
     64  A(0),
     65  Z(0),
     66  theEvaporationProbabilityPtr(0),
     67  theCoulombBarrierPtr(0),
     68  EmissionProbability(0.0),
     69  MaximalKineticEnergy(-1000.0)
     70{
     71  theLevelDensityPtr = 0;
     72  MyOwnLevelDensity = true;
     73  EvaporatedMass = 0.0;
     74  ResidualMass = CoulombBarrier = 0.0;
     75  fG4pow = G4Pow::GetInstance();
    7176}
    7277
    7378G4GEMChannel::~G4GEMChannel()
    7479{
    75     if (MyOwnLevelDensity) delete theLevelDensityPtr;
    76 }
    77 
    78 G4GEMChannel::G4GEMChannel(const G4GEMChannel & ) : G4VEvaporationChannel()
    79 {
    80     throw G4HadronicException(__FILE__, __LINE__, "G4GEMChannel::copy_costructor meant to not be accessable");
    81 }
    82 
    83 const G4GEMChannel & G4GEMChannel::operator=(const G4GEMChannel & )
    84 {
    85     throw G4HadronicException(__FILE__, __LINE__, "G4GEMChannel::operator= meant to not be accessable");
    86     return *this;
    87 }
    88 
    89 G4bool G4GEMChannel::operator==(const G4GEMChannel & right) const
    90 {
    91     return (this == (G4GEMChannel *) &right);
    92     //  return false;
    93 }
    94 
    95 G4bool G4GEMChannel::operator!=(const G4GEMChannel & right) const
    96 {
    97     return (this != (G4GEMChannel *) &right);
    98     //  return true;
     80  if (MyOwnLevelDensity) { delete theLevelDensityPtr; }
    9981}
    10082
    10183void G4GEMChannel::Initialize(const G4Fragment & fragment)
    10284{
    103   G4int anA = static_cast<G4int>(fragment.GetA());
    104   G4int aZ = static_cast<G4int>(fragment.GetZ());
    105   AResidual = anA - A;
    106   ZResidual = aZ - Z;
     85  G4int anA = fragment.GetA_asInt();
     86  G4int aZ  = fragment.GetZ_asInt();
     87  ResidualA = anA - A;
     88  ResidualZ = aZ - Z;
    10789  //G4cout << "G4GEMChannel::Initialize: Z= " << aZ << " A= " << anA
    108   //       << " Zres= " << ZResidual << " Ares= " << AResidual << G4endl;
     90  //       << " Zres= " << ResidualZ << " Ares= " << ResidualA << G4endl;
    10991
    11092  // We only take into account channels which are physically allowed
    111   if (AResidual <= 0 || ZResidual <= 0 || AResidual < ZResidual ||
    112       (AResidual == ZResidual && AResidual > 1))
     93  if (ResidualA <= 0 || ResidualZ <= 0 || ResidualA < ResidualZ ||
     94      (ResidualA == ResidualZ && ResidualA > 1))
    11395    {
    11496      CoulombBarrier = 0.0;
    … …  
    118100  else
    119101    {
    120 
    121102      // Effective excitation energy
    122103      // JMQ 071009: pairing in ExEnergy should be the one of parent compound nucleus
    … …  
    125106      // param for protons must be the same)   
    126107      //    G4double ExEnergy = fragment.GetExcitationEnergy() -
    127       //    G4PairingCorrection::GetInstance()->GetPairingCorrection(AResidual,ZResidual);
     108      //    G4PairingCorrection::GetInstance()->GetPairingCorrection(ResidualA,ResidualZ);
    128109      G4double ExEnergy = fragment.GetExcitationEnergy() -
    129110        G4PairingCorrection::GetInstance()->GetPairingCorrection(anA,aZ);
    … …  
    138119      } else {
    139120
     121        ResidualMass = G4NucleiProperties::GetNuclearMass(ResidualA, ResidualZ);
     122
    140123        // Coulomb Barrier calculation
    141         CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(AResidual,ZResidual,ExEnergy);
     124        CoulombBarrier = theCoulombBarrierPtr->GetCoulombBarrier(ResidualA,ResidualZ,ExEnergy);
    142125        //G4cout << "CBarrier(MeV)= " << CoulombBarrier/MeV << G4endl;
    143126
    144127        //Maximal kinetic energy (JMQ : at the Coulomb barrier)
    145128        MaximalKineticEnergy =
    146           CalcMaximalKineticEnergy(G4ParticleTable::GetParticleTable()->
    147                                    GetIonTable()->GetNucleusMass(aZ,anA)+ExEnergy);
     129          CalcMaximalKineticEnergy(fragment.GetGroundStateMass()+ExEnergy);
    148130        //G4cout << "MaxE(MeV)= " << MaximalKineticEnergy/MeV << G4endl;
    149 
    150131               
    151132        // Emission probability
    … …  
    161142          }
    162143      }
    163     }
    164    
     144    }   
    165145    //G4cout << "Prob= " << EmissionProbability << G4endl;
    166146    return;
    … …  
    169149G4FragmentVector * G4GEMChannel::BreakUp(const G4Fragment & theNucleus)
    170150{
    171     G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
    172     G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z,A);
    173     G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
    174 
    175    
    176     G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
    177                                                 (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
    178    
    179     momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
    180 
    181     G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
    182     EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
    183     G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
    184 #ifdef PRECOMPOUND_TEST
    185     EvaporatedFragment->SetCreatorModel(G4String("G4Evaporation"));
     151  G4double EvaporatedKineticEnergy = CalcKineticEnergy(theNucleus);
     152  G4double EvaporatedEnergy = EvaporatedKineticEnergy + EvaporatedMass;
     153 
     154  G4ThreeVector momentum(IsotropicVector(std::sqrt(EvaporatedKineticEnergy*
     155                                                   (EvaporatedKineticEnergy+2.0*EvaporatedMass))));
     156   
     157  momentum.rotateUz(theNucleus.GetMomentum().vect().unit());
     158
     159  G4LorentzVector EvaporatedMomentum(momentum,EvaporatedEnergy);
     160  EvaporatedMomentum.boost(theNucleus.GetMomentum().boostVector());
     161  G4Fragment * EvaporatedFragment = new G4Fragment(A,Z,EvaporatedMomentum);
     162  // ** And now the residual nucleus **
     163  G4double theExEnergy = theNucleus.GetExcitationEnergy();
     164  G4double theMass = theNucleus.GetGroundStateMass();
     165  G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
     166       
     167  G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
     168  ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
     169       
     170  G4Fragment * ResidualFragment = new G4Fragment( ResidualA, ResidualZ, ResidualMomentum );
     171   
     172  G4FragmentVector * theResult = new G4FragmentVector;
     173   
     174#ifdef debug
     175  G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
     176  G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
     177  if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 10.0*eV) {
     178    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
     179    G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :"
     180           <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
     181           << " MeV" << G4endl;
     182  }
     183  if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 10.0*eV ||
     184      std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 10.0*eV ||
     185      std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 10.0*eV ) {
     186    G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
     187    G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :"
     188           <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
     189           << " MeV" << G4endl;   
     190       
     191  }
    186192#endif
    187     // ** And now the residual nucleus **
    188     G4double theExEnergy = theNucleus.GetExcitationEnergy();
    189     G4double theMass = G4ParticleTable::GetParticleTable()->
    190         GetIonTable()->GetNucleusMass(static_cast<G4int>(theNucleus.GetZ()),
    191                                       static_cast<G4int>(theNucleus.GetA()));
    192     G4double ResidualEnergy = theMass + (theExEnergy - EvaporatedKineticEnergy) - EvaporatedMass;
    193        
    194     G4LorentzVector ResidualMomentum(-momentum,ResidualEnergy);
    195     ResidualMomentum.boost(theNucleus.GetMomentum().boostVector());
    196        
    197     G4Fragment * ResidualFragment = new G4Fragment( AResidual, ZResidual, ResidualMomentum );
    198 
    199 #ifdef PRECOMPOUND_TEST
    200     ResidualFragment->SetCreatorModel(G4String("ResidualNucleus"));
    201 #endif
    202    
    203    
    204     G4FragmentVector * theResult = new G4FragmentVector;
    205    
    206 #ifdef debug
    207     G4double Efinal = ResidualMomentum.e() + EvaporatedMomentum.e();
    208     G4ThreeVector Pfinal = ResidualMomentum.vect() + EvaporatedMomentum.vect();
    209     if (std::abs(Efinal-theNucleus.GetMomentum().e()) > 10.0*eV) {
    210         G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: ENERGY @@@@@@@@@@@@@@@@" << G4endl;
    211         G4cout << "Initial : " << theNucleus.GetMomentum().e()/MeV << " MeV    Final :"
    212                <<Efinal/MeV << " MeV    Delta: " <<  (Efinal-theNucleus.GetMomentum().e())/MeV
    213                << " MeV" << G4endl;
    214     }
    215     if (std::abs(Pfinal.x()-theNucleus.GetMomentum().x()) > 10.0*eV ||
    216         std::abs(Pfinal.y()-theNucleus.GetMomentum().y()) > 10.0*eV ||
    217         std::abs(Pfinal.z()-theNucleus.GetMomentum().z()) > 10.0*eV ) {
    218         G4cout << "@@@@@@@@@@@@@@@@@@@@@ G4GEMChanel: MOMENTUM @@@@@@@@@@@@@@@@" << G4endl;
    219         G4cout << "Initial : " << theNucleus.GetMomentum().vect() << " MeV    Final :"
    220                <<Pfinal/MeV << " MeV    Delta: " <<  Pfinal-theNucleus.GetMomentum().vect()
    221                << " MeV" << G4endl;   
    222        
    223     }
    224 #endif
    225     theResult->push_back(EvaporatedFragment);
    226     theResult->push_back(ResidualFragment);
    227     return theResult;
     193  theResult->push_back(EvaporatedFragment);
     194  theResult->push_back(ResidualFragment);
     195  return theResult;
    228196}
    229 
    230 
    231197
    232198G4double G4GEMChannel::CalcMaximalKineticEnergy(const G4double NucleusTotalE)
    … …  
    234200//JMQ this is not the assimptotic kinetic energy but the one at the Coulomb barrier
    235201{
    236   G4double ResidualMass = G4ParticleTable::GetParticleTable()->
    237     GetIonTable()->GetNucleusMass( ZResidual, AResidual );
    238   G4double EvaporatedMass = G4ParticleTable::GetParticleTable()->
    239     GetIonTable()->GetNucleusMass( Z, A );
    240        
    241202  G4double T = (NucleusTotalE*NucleusTotalE + EvaporatedMass*EvaporatedMass - ResidualMass*ResidualMass)/
    242     (2.0*NucleusTotalE) -
    243     EvaporatedMass - CoulombBarrier;
     203    (2.0*NucleusTotalE) - EvaporatedMass - CoulombBarrier;
    244204       
    245205  return T;
    246206}
    247207
    248 
    249 
    250 
    251208G4double G4GEMChannel::CalcKineticEnergy(const G4Fragment & fragment)
    252209// Samples fragment kinetic energy.
    253210{
    254211  G4double U = fragment.GetExcitationEnergy();
    255   G4double NuclearMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetNucleusMass(Z,A);
    256212 
    257213  G4double Alpha = theEvaporationProbabilityPtr->CalcAlphaParam(fragment);
    … …  
    260216  G4double Normalization = theEvaporationProbabilityPtr->GetNormalization();
    261217
    262 //                       ***RESIDUAL***
     218  //                       ***RESIDUAL***
    263219  //JMQ (September 2009) the following quantities  refer to the RESIDUAL:
    264   G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(AResidual,ZResidual);
    265   G4double Ux = (2.5 + 150.0/AResidual)*MeV;
     220  G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection(ResidualA,ResidualZ);
     221  G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
    266222  G4double Ex = Ux + delta0;
    267223  G4double InitialLevelDensity;
    … …  
    271227  //JMQ (September 2009) the following quantities   refer to the PARENT:
    272228 
    273   G4double deltaCN = G4PairingCorrection::GetInstance()->GetPairingCorrection(static_cast<G4int>(fragment.GetA()),
    274                                                                               static_cast<G4int>(fragment.GetZ()));
    275   G4double aCN = theLevelDensityPtr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
    276                                                            static_cast<G4int>(fragment.GetZ()),U-deltaCN);   
     229  G4double deltaCN =
     230    G4PairingCorrection::GetInstance()->GetPairingCorrection(fragment.GetA_asInt(),
     231                                                             fragment.GetZ_asInt());
     232  G4double aCN = theLevelDensityPtr->LevelDensityParameter(fragment.GetA_asInt(),
     233                                                           fragment.GetZ_asInt(),U-deltaCN);   
    277234  G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
    278235  G4double ExCN = UxCN + deltaCN;
    279236  G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
    280   G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
    281237  //                       ***end PARENT***
    282238 
    … …  
    284240  if ( U < ExCN )
    285241    {
     242      G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0
     243                                  - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
    286244      InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
    287245    }
    288246  else
    289247    {
    290       InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
     248      G4double x  = U-deltaCN;
     249      G4double x1 = std::sqrt(aCN*x);
     250      InitialLevelDensity = (pi/12.0)*std::exp(2*x1)/(x*std::sqrt(x1));
     251      //InitialLevelDensity =
     252      //(pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
    291253    }
    292254 
    293255  const G4double Spin = theEvaporationProbabilityPtr->GetSpin();
    294 //JMQ  BIG BUG fixed: hbarc instead of hbar_Planck !!!!
    295 //     it was fixed in total probability (for this channel) but remained still here!!
    296 //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
    297     G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
    298 //
    299 //JMQ  fix on Rb and  geometrical cross sections according to Furihata's paper
    300 //                      (JAERI-Data/Code 2001-105, p6)
     256  //JMQ  BIG BUG fixed: hbarc instead of hbar_Planck !!!!
     257  //     it was fixed in total probability (for this channel) but remained still here!!
     258  //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
     259  G4double g = (2.0*Spin+1.0)*EvaporatedMass/(pi2* hbarc*hbarc);
     260  //
     261  //JMQ  fix on Rb and  geometrical cross sections according to Furihata's paper
     262  //                      (JAERI-Data/Code 2001-105, p6)
    301263  G4double Rb = 0.0;
    302     if (A > 4)
    303     {
    304       G4double Ad = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
    305       G4double Aj = std::pow(G4double(A),1.0/3.0);
     264  if (A > 4)
     265    {
     266      G4double Ad = fG4pow->Z13(ResidualA);
     267      G4double Aj = fG4pow->Z13(A);
     268      //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
    306269      Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
    307270      Rb *= fermi;
    308271    }
    309     else if (A>1)
    310       {
    311         //      G4double R1 = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
    312         //      G4double R2 = std::pow(G4double(A),1.0/3.0);
    313       G4double Ad = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
    314       G4double Aj = std::pow(G4double(A),1.0/3.0);
    315       //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
     272  else if (A>1)
     273    {
     274      G4double Ad = fG4pow->Z13(ResidualA);
     275      G4double Aj = fG4pow->Z13(A);
    316276      Rb=1.5*(Aj+Ad)*fermi;
    317277    }
    318     else
    319     {
    320       G4double Ad = std::pow(G4double(fragment.GetA()-A),1.0/3.0);
    321       //        RN = 1.5*fermi;
     278  else
     279    {
     280      G4double Ad = fG4pow->Z13(ResidualA);
    322281      Rb = 1.5*Ad*fermi;
    323282    }
    324     //   G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.);
    325     G4double GeometricalXS = pi*Rb*Rb;
    326    
    327 
    328     G4double ConstantFactor = g*GeometricalXS*Alpha/InitialLevelDensity;
    329     ConstantFactor *= pi/12.0;
    330     //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
    331     //      (obtained by energy-momentum conservation).
    332     //      In general smaller than U-theSeparationEnergy
    333     G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
    334     G4double KineticEnergy;
    335     G4double Probability;
    336 
    337     do
    338       {
    339         KineticEnergy =  CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
    340         Probability = ConstantFactor*(KineticEnergy + Beta);
    341         G4double a = theLevelDensityPtr->LevelDensityParameter(AResidual,ZResidual,theEnergy-KineticEnergy-delta0);
    342         G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
    343         //JMQ fix in units
    344         G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
     283  //   G4double GeometricalXS = pi*RN*RN*std::pow(G4double(fragment.GetA()-A),2./3.);
     284  G4double GeometricalXS = pi*Rb*Rb;
     285   
     286  G4double ConstantFactor = g*GeometricalXS*Alpha/InitialLevelDensity;
     287  ConstantFactor *= pi/12.0;
     288  //JMQ : this is the assimptotic maximal kinetic energy of the ejectile
     289  //      (obtained by energy-momentum conservation).
     290  //      In general smaller than U-theSeparationEnergy
     291  G4double theEnergy = MaximalKineticEnergy + CoulombBarrier;
     292  G4double KineticEnergy;
     293  G4double Probability;
     294
     295  do
     296    {
     297      KineticEnergy =  CoulombBarrier + G4UniformRand()*MaximalKineticEnergy;
     298      Probability = ConstantFactor*(KineticEnergy + Beta);
     299      G4double a =
     300        theLevelDensityPtr->LevelDensityParameter(ResidualA,ResidualZ,theEnergy-KineticEnergy-delta0);
     301      G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
     302      //JMQ fix in units
    345303       
    346         if ( theEnergy-KineticEnergy < Ex)
    347           {
    348             Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
    349           }
    350         else
    351           {
    352             Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
    353               std::pow(a*std::pow(theEnergy-KineticEnergy-delta0,5.0),1.0/4.0);
    354           }
    355         Probability /= Normalization;
    356       }
    357     while (G4UniformRand() > Probability);
    358    
    359     return KineticEnergy;
     304      if ( theEnergy-KineticEnergy < Ex)
     305        {
     306          G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0
     307                                - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
     308          Probability *= std::exp((theEnergy-KineticEnergy-E0)/T)/T;
     309        }
     310      else
     311        {
     312          Probability *= std::exp(2*std::sqrt(a*(theEnergy-KineticEnergy-delta0)))/
     313            std::pow(a*fG4pow->powN(theEnergy-KineticEnergy-delta0,5), 0.25);
     314        }
     315    }
     316  while (Normalization*G4UniformRand() > Probability);
     317   
     318  return KineticEnergy;
    360319}
    361320
    … …  
    365324    // By default Magnitude = 1.0
    366325{
    367     G4double CosTheta = 1.0 - 2.0*G4UniformRand();
    368     G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
    369     G4double Phi = twopi*G4UniformRand();
    370     G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
    371                          Magnitude*std::sin(Phi)*SinTheta,
    372                          Magnitude*CosTheta);
    373     return Vector;
    374 }
    375 
    376 
    377 
     326  G4double CosTheta = 1.0 - 2.0*G4UniformRand();
     327  G4double SinTheta = std::sqrt(1.0 - CosTheta*CosTheta);
     328  G4double Phi = twopi*G4UniformRand();
     329  G4ThreeVector Vector(Magnitude*std::cos(Phi)*SinTheta,
     330                       Magnitude*std::sin(Phi)*SinTheta,
     331                       Magnitude*CosTheta);
     332  return Vector;
     333}
     334
     335
     336

  • trunk/source/processes/hadronic/models/de_excitation/gem_evaporation/src/G4GEMProbability.cc

    r1340 r1347  
    2424// ********************************************************************
    2525//
    26 // $Id: G4GEMProbability.cc,v 1.13 2010/05/19 10:21:16 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03-ref-09 $
     26// $Id: G4GEMProbability.cc,v 1.15 2010/11/05 14:43:27 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-ref-00 $
    2828//
    2929//---------------------------------------------------------------------
    … …  
    6565}
    6666
    67 G4GEMProbability:: G4GEMProbability(const G4int anA, const G4int aZ, const G4double aSpin) :
     67G4GEMProbability:: G4GEMProbability(G4int anA, G4int aZ, G4double aSpin) :
    6868  theA(anA), theZ(aZ), Spin(aSpin), theCoulombBarrierPtr(0),
    6969  ExcitationEnergies(0), ExcitationSpins(0), ExcitationLifetimes(0), Normalization(1.0)
    … …  
    9595
    9696G4double G4GEMProbability::EmissionProbability(const G4Fragment & fragment,
    97                                                const G4double MaximalKineticEnergy)
     97                                               G4double MaximalKineticEnergy)
    9898{
    9999  G4double EmissionProbability = 0.0;
    … …  
    130130
    131131G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
    132                                            const G4double MaximalKineticEnergy,
    133                                            const G4double V)
     132                                           G4double MaximalKineticEnergy,
     133                                           G4double V)
    134134
    135135// Calculate integrated probability (width) for evaporation channel
    … …  
    146146  G4double Alpha = CalcAlphaParam(fragment);
    147147  G4double Beta = CalcBetaParam(fragment);
    148  
    149  
     148   
    150149  //                       ***RESIDUAL***
    151150  //JMQ (September 2009) the following quantities refer to the RESIDUAL:
    … …  
    159158  G4double T  = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
    160159  //JMQ fixed bug in units
    161   G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
     160  G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0
     161        - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
    162162  //                      ***end RESIDUAL ***
    163163 
    … …  
    165165  //JMQ (September 2009) the following quantities refer to the PARENT:
    166166     
    167   G4double deltaCN=fPairCorr->GetPairingCorrection(A, Z);                                     
    168   G4double aCN = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
    169   G4double UxCN = (2.5 + 150.0/G4double(A))*MeV;
    170   G4double ExCN = UxCN + deltaCN;
    171   G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
    172   //JMQ fixed bug in units
    173   G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
    174 //                       ***end PARENT***
     167  G4double deltaCN = fPairCorr->GetPairingCorrection(A, Z);                                   
     168  G4double aCN     = theEvapLDPptr->LevelDensityParameter(A, Z, U-deltaCN);
     169  G4double UxCN    = (2.5 + 150.0/G4double(A))*MeV;
     170  G4double ExCN    = UxCN + deltaCN;
     171  G4double TCN     = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
     172  //                     ***end PARENT***
    175173
    176174  G4double Width;
    … …  
    231229  //end of JMQ fix on Rb by 190709
    232230 
    233 
    234 
    235 //JMQ 160909 fix: initial level density must be calculated according to the
    236 // conditions at the initial compound nucleus
    237 // (it has been removed from previous "if" for the residual)
    238 
    239    if ( U < ExCN )
    240       {
    241         InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
    242       }
    243     else
    244       {
    245         //        InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
    246         //VI speedup
    247         G4double x  = U-deltaCN;
    248         G4double x2 = x*x;
    249         G4double x3 = x2*x;
    250         InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*x))/std::sqrt(std::sqrt(aCN*x2*x3));
    251       }
    252  //
    253 
     231  //JMQ 160909 fix: initial level density must be calculated according to the
     232  // conditions at the initial compound nucleus
     233  // (it has been removed from previous "if" for the residual)
     234
     235  if ( U < ExCN )
     236    {
     237      //JMQ fixed bug in units
     238      //VI moved the computation here
     239      G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0
     240                                  - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
     241      InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
     242    }
     243  else
     244    {
     245      //InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
     246      //VI speedup
     247      G4double x  = U-deltaCN;
     248      G4double x1 = std::sqrt(aCN*x);
     249      InitialLevelDensity = (pi/12.0)*std::exp(2*x1)/(x*std::sqrt(x1));
     250    }
    254251
    255252  //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
    … …  
    257254  Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
    258255   
    259 
    260256  return Width;
    261257}
    262258
    263 /*
    264 
    265 G4double G4GEMProbability::CalcProbability(const G4Fragment & fragment,
    266                                            const G4double MaximalKineticEnergy,
    267                                            const G4double V)
    268 // Calculate integrated probability (width) for evaporation channel
    269 {
    270   G4double ResidualA = static_cast<G4double>(fragment.GetA() - theA);
    271   G4double ResidualZ = static_cast<G4double>(fragment.GetZ() - theZ);
    272   G4double U = fragment.GetExcitationEnergy();
    273  
    274   G4double NuclearMass = G4ParticleTable::GetParticleTable()->
    275     GetIonTable()->GetNucleusMass(theZ,theA);
    276  
    277  
    278   G4double Alpha = CalcAlphaParam(fragment);
    279   G4double Beta = CalcBetaParam(fragment);
    280  
    281  
    282   //                       ***RESIDUAL***
    283   //JMQ (September 2009) the following quantities refer to the RESIDUAL:
    284  
    285   G4double delta0 = G4PairingCorrection::GetInstance()->GetPairingCorrection( static_cast<G4int>( ResidualA ) , static_cast<G4int>( ResidualZ ) ); 
    286  
    287   G4double a = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(ResidualA),
    288                                                     static_cast<G4int>(ResidualZ),MaximalKineticEnergy+V-delta0);
    289   G4double Ux = (2.5 + 150.0/ResidualA)*MeV;
    290   G4double Ex = Ux + delta0;
    291   G4double T = 1.0/(std::sqrt(a/Ux) - 1.5/Ux);
    292   //JMQ fixed bug in units
    293   G4double E0 = Ex - T*(std::log(T/MeV) - std::log(a*MeV)/4.0 - 1.25*std::log(Ux/MeV) + 2.0*std::sqrt(a*Ux));
    294   //                      ***end RESIDUAL ***
    295  
    296   //                       ***PARENT***
    297   //JMQ (September 2009) the following quantities refer to the PARENT:
    298      
    299   G4double deltaCN=G4PairingCorrection::GetInstance()->
    300     GetPairingCorrection(static_cast<G4int>(fragment.GetA()),static_cast<G4int>(fragment.GetZ()));                                     
    301   G4double aCN = theEvapLDPptr->LevelDensityParameter(static_cast<G4int>(fragment.GetA()),
    302                                                       static_cast<G4int>(fragment.GetZ()),U-deltaCN);
    303   G4double UxCN = (2.5 + 150.0/fragment.GetA())*MeV;
    304   G4double ExCN = UxCN + deltaCN;
    305   G4double TCN = 1.0/(std::sqrt(aCN/UxCN) - 1.5/UxCN);
    306   //JMQ fixed bug in units
    307   G4double E0CN = ExCN - TCN*(std::log(TCN/MeV) - std::log(aCN*MeV)/4.0 - 1.25*std::log(UxCN/MeV) + 2.0*std::sqrt(aCN*UxCN));
    308 //                       ***end PARENT***
    309 
    310   G4double Width;
    311   G4double InitialLevelDensity;
    312   G4double t = MaximalKineticEnergy/T;
    313   if ( MaximalKineticEnergy < Ex ) {
    314     //JMQ 190709 bug in I1 fixed (T was  missing)
    315     Width = (I1(t,t)*T + (Beta+V)*I0(t))/std::exp(E0/T);
    316 //JMQ 160909 fix:  InitialLevelDensity has been taken away (different conditions for initial CN..)
    317   } else {
    318     G4double tx = Ex/T;
    319     G4double s = 2.0*std::sqrt(a*(MaximalKineticEnergy-delta0));
    320     G4double sx = 2.0*std::sqrt(a*(Ex-delta0));
    321     Width = I1(t,tx)*T/std::exp(E0/T) + I3(s,sx)*std::exp(s)/(std::sqrt(2.0)*a);
    322     // For charged particles (Beta+V) = 0 beacuse Beta = -V
    323     if (theZ == 0) {
    324       Width += (Beta+V)*(I0(tx)/std::exp(E0/T) + 2.0*std::sqrt(2.0)*I2(s,sx)*std::exp(s));
    325     }
    326   }
    327  
    328   //JMQ 14/07/2009 BIG BUG : NuclearMass is in MeV => hbarc instead of hbar_planck must be used
    329   //    G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbar_Planck*hbar_Planck);
    330   G4double g = (2.0*Spin+1.0)*NuclearMass/(pi2* hbarc*hbarc);
    331  
    332   //JMQ 190709 fix on Rb and  geometrical cross sections according to Furihata's paper
    333   //                      (JAERI-Data/Code 2001-105, p6)
    334   //    G4double RN = 0.0;
    335   G4double Rb = 0.0;
    336   if (theA > 4)
    337     {
    338       //        G4double R1 = std::pow(ResidualA,1.0/3.0);
    339       //        G4double R2 = std::pow(G4double(theA),1.0/3.0);
    340       G4double Ad = std::pow(ResidualA,1.0/3.0);
    341       G4double Aj = std::pow(G4double(theA),1.0/3.0);
    342       //        RN = 1.12*(R1 + R2) - 0.86*((R1+R2)/(R1*R2));
    343       Rb = 1.12*(Aj + Ad) - 0.86*((Aj+Ad)/(Aj*Ad))+2.85;
    344       Rb *= fermi;
    345     }
    346   else if (theA>1)
    347     {
    348       G4double Ad = std::pow(ResidualA,1.0/3.0);
    349       G4double Aj = std::pow(G4double(theA),1.0/3.0);
    350       Rb=1.5*(Aj+Ad)*fermi;
    351     }
    352   else
    353     {
    354       G4double Ad = std::pow(ResidualA,1.0/3.0);
    355       Rb = 1.5*Ad*fermi;
    356     }
    357   //   G4double GeometricalXS = pi*RN*RN*std::pow(ResidualA,2./3.);
    358   G4double GeometricalXS = pi*Rb*Rb;
    359   //end of JMQ fix on Rb by 190709
    360  
    361 
    362 
    363 //JMQ 160909 fix: initial level density must be calculated according to the
    364 // conditions at the initial compound nucleus
    365 // (it has been removed from previous "if" for the residual)
    366 
    367    if ( U < ExCN )
    368       {
    369         InitialLevelDensity = (pi/12.0)*std::exp((U-E0CN)/TCN)/TCN;
    370       }
    371     else
    372       {
    373         InitialLevelDensity = (pi/12.0)*std::exp(2*std::sqrt(aCN*(U-deltaCN)))/std::pow(aCN*std::pow(U-deltaCN,5.0),1.0/4.0);
    374       }
    375  //
    376 
    377 
    378   //JMQ 190709 BUG : pi instead of sqrt(pi) must be here according to Furihata's report:
    379   //    Width *= std::sqrt(pi)*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
    380   Width *= pi*g*GeometricalXS*Alpha/(12.0*InitialLevelDensity);
    381    
    382 
    383   return Width;
    384 }
    385 */
    386 
    387 G4double G4GEMProbability::I3(const G4double s, const G4double sx)
     259G4double G4GEMProbability::I3(G4double s, G4double sx)
    388260{
    389261  G4double s2 = s*s;
Note: See TracChangeset for help on using the changeset viewer.