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/src
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • 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.