Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/pre_equilibrium/exciton_model/src/G4PreCompoundTransitions.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4PreCompoundTransitions.cc,v 1.22 2009/11/21 18:03:13 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PreCompoundTransitions.cc,v 1.27 2010/10/20 00:47:46 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    3737//
    3838// Modified: 
    39 // 16.02.2008 J. M. Quesada fixed bugs
    40 // 06.09.2008 J. M. Quesada added external choices for:
     39// 16.02.2008 J.M.Quesada fixed bugs
     40// 06.09.2008 J.M.Quesada added external choices for:
    4141//                      - "never go back"  hipothesis (useNGB=true)
    4242//                      -  CEM transition probabilities (useCEMtr=true)
    43 
    44 // 30.10.09 J.M.Quesada: CEM transition probabilities have been renormalized
     43// 30.10.2009 J.M.Quesada: CEM transition probabilities have been renormalized
    4544//                       (IAEA benchmark)
    46 //
     45// 20.08.2010 V.Ivanchenko move constructor and destructor to the source and
     46//                         optimise the code
     47//
     48
    4749#include "G4PreCompoundTransitions.hh"
    4850#include "G4HadronicException.hh"
    49 
    50 const G4PreCompoundTransitions & G4PreCompoundTransitions::
    51 operator=(const G4PreCompoundTransitions &)
     51#include "G4PreCompoundParameters.hh"
     52#include "G4Proton.hh"
     53#include "Randomize.hh"
     54#include "G4Pow.hh"
     55
     56G4PreCompoundTransitions::G4PreCompoundTransitions()
    5257{
    53   throw G4HadronicException(__FILE__, __LINE__, "G4PreCompoundTransitions::operator= meant to not be accessable");
    54   return *this;
     58  proton = G4Proton::Proton();
     59  FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
     60  r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0();
     61  aLDP = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     62  g4pow = G4Pow::GetInstance();
    5563}
    5664
    57 
    58 G4bool G4PreCompoundTransitions::operator==(const G4PreCompoundTransitions &) const
    59 {
    60   return false;
    61 }
    62 
    63 G4bool G4PreCompoundTransitions::operator!=(const G4PreCompoundTransitions &) const
    64 {
    65   return true;
    66 }
    67 
    68 
     65G4PreCompoundTransitions::~G4PreCompoundTransitions()
     66{}
     67
     68// Calculates transition probabilities with
     69// DeltaN = +2 (Trans1) -2 (Trans2) and 0 (Trans3)
    6970G4double G4PreCompoundTransitions::
    7071CalculateProbability(const G4Fragment & aFragment)
    7172{
    72   //G4cout<<"In G4PreCompoundTransitions.cc  useNGB="<<useNGB<<G4endl;
    73   //G4cout<<"In G4PreCompoundTransitions.cc  useCEMtr="<<useCEMtr<<G4endl;
    74 
    75   // Fermi energy
    76   const G4double FermiEnergy = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
     73  // Number of holes
     74  G4int H = aFragment.GetNumberOfHoles();
     75  // Number of Particles
     76  G4int P = aFragment.GetNumberOfParticles();
     77  // Number of Excitons
     78  G4int N = P+H;
     79  // Nucleus
     80  G4int A = aFragment.GetA_asInt();
     81  G4int Z = aFragment.GetZ_asInt();
     82  G4double U = aFragment.GetExcitationEnergy();
     83
     84  //G4cout << aFragment << G4endl;
    7785 
    78   // Nuclear radius
    79   const G4double r0 = G4PreCompoundParameters::GetAddress()->GetTransitionsr0();
    80    
    81   // In order to calculate the level density parameter
    82   // G4EvaporationLevelDensityParameter theLDP;
    83 
    84   // Number of holes
    85   G4double H = aFragment.GetNumberOfHoles();
    86   // Number of Particles
    87   G4double P = aFragment.GetNumberOfParticles();
    88   // Number of Excitons
    89   G4double N = P+H;
    90   // Nucleus
    91   G4double A = aFragment.GetA();
    92   G4double Z = static_cast<G4double>(aFragment.GetZ());
    93   G4double U = aFragment.GetExcitationEnergy();
    94  
    95   if(U<10*eV) return 0.0;
     86  if(U < 10*eV) { return 0.0; }
    9687 
    9788  //J. M. Quesada (Feb. 08) new physics
    98   // OPT=1 Transitions are calculated according to Gudima's paper (original in G4PreCompound from VL)
     89  // OPT=1 Transitions are calculated according to Gudima's paper
     90  //       (original in G4PreCompound from VL)
    9991  // OPT=2 Transitions are calculated according to Gupta's formulae
    10092  //
    101  
    102  
    103  
    10493  if (useCEMtr){
    10594
    106    
    10795    // Relative Energy (T_{rel})
    108     G4double RelativeEnergy = (8.0/5.0)*FermiEnergy + U/N;
     96    G4double RelativeEnergy = 1.6*FermiEnergy + U/G4double(N);
    10997   
    11098    // Sample kind of nucleon-projectile
    11199    G4bool ChargedNucleon(false);
    112100    G4double chtest = 0.5;
    113     if (P > 0) chtest = aFragment.GetNumberOfCharged()/P;
    114     if (G4UniformRand() < chtest) ChargedNucleon = true;
     101    if (P > 0) {
     102      chtest = G4double(aFragment.GetNumberOfCharged())/G4double(P);
     103    }
     104    if (G4UniformRand() < chtest) { ChargedNucleon = true; }
    115105   
    116106    // Relative Velocity:
    117107    // <V_{rel}>^2
    118108    G4double RelativeVelocitySqr(0.0);
    119     if (ChargedNucleon) RelativeVelocitySqr = 2.0*RelativeEnergy/proton_mass_c2;
    120     else RelativeVelocitySqr = 2.0*RelativeEnergy/neutron_mass_c2;
     109    if (ChargedNucleon) {
     110      RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::proton_mass_c2;
     111    } else {
     112      RelativeVelocitySqr = 2.0*RelativeEnergy/CLHEP::neutron_mass_c2;
     113    }
    121114   
    122115    // <V_{rel}>
    … …  
    124117   
    125118    // Proton-Proton Cross Section
    126     G4double ppXSection = (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)*millibarn;
     119    G4double ppXSection =
     120      (10.63/RelativeVelocitySqr - 29.92/RelativeVelocity + 42.9)
     121      * CLHEP::millibarn;
    127122    // Proton-Neutron Cross Section
    128     G4double npXSection = (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)*millibarn;
     123    G4double npXSection =
     124      (34.10/RelativeVelocitySqr - 82.20/RelativeVelocity + 82.2)
     125      * CLHEP::millibarn;
    129126   
    130127    // Averaged Cross Section: \sigma(V_{rel})
    … …  
    134131      {
    135132        //JMQ: small bug fixed
    136         //      AveragedXSection = ((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection) / (A-1.0);
    137         AveragedXSection = ((Z-1.0) * ppXSection + (A-Z) * npXSection) / (A-1.0);
     133        //AveragedXSection=((Z-1.0) * ppXSection + (A-Z-1.0) * npXSection)/(A-1.0);
     134        AveragedXSection = ((Z-1)*ppXSection + (A-Z)*npXSection)/G4double(A-1);
    138135      }
    139136    else
    140137      {
    141         AveragedXSection = ((A-Z-1.0) * ppXSection + Z * npXSection) / (A-1.0);
     138        AveragedXSection = ((A-Z-1)*ppXSection + Z*npXSection)/G4double(A-1);
     139        //AveragedXSection = ((A-Z-1)*npXSection + Z*ppXSection)/G4double(A-1);
    142140      }
    143141   
    … …  
    146144   
    147145    // This factor is introduced to take into account the Pauli principle
    148     G4double PauliFactor = 1.0 - (7.0/5.0)*FermiRelRatio;
    149     if (FermiRelRatio > 0.5) PauliFactor += (2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
    150    
     146    G4double PauliFactor = 1.0 - 1.4*FermiRelRatio;
     147    if (FermiRelRatio > 0.5) {
     148      G4double x = 2.0 - 1.0/FermiRelRatio;
     149      PauliFactor += 0.4*FermiRelRatio*x*x*std::sqrt(x);
     150      //PauliFactor +=
     151      //(2.0/5.0)*FermiRelRatio*std::pow(2.0 - (1.0/FermiRelRatio), 5.0/2.0);
     152    }
    151153    // Interaction volume
    152     //  G4double Vint = (4.0/3.0)*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
    153     G4double xx=2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity);
    154     G4double Vint = (4.0/3.0)*pi*xx*xx*xx;
     154    //  G4double Vint = (4.0/3.0)
     155    //*pi*std::pow(2.0*r0 + hbarc/(proton_mass_c2*RelativeVelocity) , 3.0);
     156    G4double xx = 2.0*r0 + hbarc/(CLHEP::proton_mass_c2*RelativeVelocity);
     157    //    G4double Vint = (4.0/3.0)*CLHEP::pi*xx*xx*xx;
     158    G4double Vint = CLHEP::pi*xx*xx*xx/0.75;
    155159   
    156160    // Transition probability for \Delta n = +2
    157161   
    158     TransitionProb1 = AveragedXSection*PauliFactor*std::sqrt(2.0*RelativeEnergy/proton_mass_c2)/Vint;
    159 
    160 //JMQ 281009  phenomenological factor in order to increase equilibrium contribution
    161 //   G4double factor=5.0;
    162 //   TransitionProb1 *= factor;
    163 //
    164     if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
    165    
    166     G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     162    TransitionProb1 = AveragedXSection*PauliFactor
     163      *std::sqrt(2.0*RelativeEnergy/CLHEP::proton_mass_c2)/Vint;
     164
     165    //JMQ 281009  phenomenological factor in order to increase
     166    //   equilibrium contribution
     167    //   G4double factor=5.0;
     168    //   TransitionProb1 *= factor;
     169    //
     170    if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }
     171   
    167172    // GE = g*E where E is Excitation Energy
    168     G4double GE = (6.0/pi2)*a*A*U;
    169    
    170     G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
    171    
    172     //G4bool NeverGoBack(false);
    173     G4bool NeverGoBack;
    174     if(useNGB)  NeverGoBack=true;
    175     else NeverGoBack=false;
    176    
     173    G4double GE = (6.0/pi2)*aLDP*A*U;
     174   
     175    //G4double Fph = ((P*P+H*H+P-H)/4.0 - H/2.0);
     176    G4double Fph = G4double(P*P+H*H+P-3*H)/4.0;
     177   
     178    G4bool NeverGoBack(false);
     179    if(useNGB) { NeverGoBack=true; }
    177180   
    178181    //JMQ/AH  bug fixed: if (U-Fph < 0.0) NeverGoBack = true;
    179     if (GE-Fph < 0.0) NeverGoBack = true;
     182    if (GE-Fph < 0.0) { NeverGoBack = true; }
    180183   
    181184    // F(p+1,h+1)
    182185    G4double Fph1 = Fph + N/2.0;
    183186   
    184     G4double ProbFactor = std::pow((GE-Fph)/(GE-Fph1),N+1.0);
    185    
     187    G4double ProbFactor = g4pow->powN((GE-Fph)/(GE-Fph1),N+1);
    186188   
    187189    if (NeverGoBack)
    188190      {
    189       TransitionProb2 = 0.0;
    190       TransitionProb3 = 0.0;
     191        TransitionProb2 = 0.0;
     192        TransitionProb3 = 0.0;
    191193      }
    192194    else
    193195      {
    194196        // Transition probability for \Delta n = -2 (at F(p,h) = 0)
    195         TransitionProb2 = TransitionProb1 * ProbFactor * (P*H*(N+1.0)*(N-2.0))/((GE-Fph)*(GE-Fph));
    196         if (TransitionProb2 < 0.0) TransitionProb2 = 0.0;
     197        TransitionProb2 =
     198          TransitionProb1 * ProbFactor * (P*H*(N+1)*(N-2))/((GE-Fph)*(GE-Fph));
     199        if (TransitionProb2 < 0.0) { TransitionProb2 = 0.0; }
    197200       
    198201        // Transition probability for \Delta n = 0 (at F(p,h) = 0)
    199         TransitionProb3 = TransitionProb1* ((N+1.0)/N) * ProbFactor  * (P*(P-1.0) + 4.0*P*H + H*(H-1.0))/(GE-Fph);
    200         if (TransitionProb3 < 0.0) TransitionProb3 = 0.0;
    201       }
    202    
    203     //  G4cout<<"U = "<<U<<G4endl;
    204     //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
    205     //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
    206     return TransitionProb1 + TransitionProb2 + TransitionProb3;
    207   }
    208  
    209   else  {
    210     //JMQ: Transition probabilities from Gupta's work
    211    
    212     G4double a = G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     202        TransitionProb3 = TransitionProb1*(N+1)* ProbFactor 
     203          * (P*(P-1) + 4.0*P*H + H*(H-1))/(N*(GE-Fph));
     204        if (TransitionProb3 < 0.0) { TransitionProb3 = 0.0; }
     205      }
     206  } else {
     207    //JMQ: Transition probabilities from Gupta's work   
    213208    // GE = g*E where E is Excitation Energy
    214     G4double GE = (6.0/pi2)*a*A*U;
     209    G4double GE = (6.0/pi2)*aLDP*A*U;
    215210 
    216211    G4double Kmfp=2.;
    217212       
    218     TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
    219     if (TransitionProb1 < 0.0) TransitionProb1 = 0.0;
    220    
    221     if (useNGB){
    222       TransitionProb2=0.;
    223       TransitionProb3=0.;
    224     }
    225     else{       
    226       if (N<=1) TransitionProb2=0. ;   
    227       else  TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);     
     213    //TransitionProb1=1./Kmfp*3./8.*1./c_light*1.0e-9*(1.4e+21*U-2./(N+1)*6.0e+18*U*U);
     214    TransitionProb1 = 3.0e-9*(1.4e+21*U - 1.2e+19*U*U/G4double(N+1))
     215      /(8*Kmfp*CLHEP::c_light);
     216    if (TransitionProb1 < 0.0) { TransitionProb1 = 0.0; }
     217
     218    TransitionProb2=0.;
     219    TransitionProb3=0.;
     220   
     221    if (!useNGB && N > 1) {
     222      // TransitionProb2=1./Kmfp*3./8.*1./c_light*1.0e-9*(N-1.)*(N-2.)*P*H/(GE*GE)*(1.4e+21*U - 2./(N-1)*6.0e+18*U*U);     
     223      TransitionProb2 =
     224        3.0e-9*(N-2)*P*H*(1.4e+21*U*(N-1) - 1.2e+19*U*U)/(8*Kmfp*c_light*GE*GE);     
    228225      if (TransitionProb2 < 0.0) TransitionProb2 = 0.0;
    229       TransitionProb3=0.;
    230     }
    231    
    232       //  G4cout<<"U = "<<U<<G4endl;
    233     //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
    234     //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2<<"  l0 ="<< TransitionProb3<<G4endl;
    235     return TransitionProb1 + TransitionProb2 + TransitionProb3;
     226    }
    236227  }
     228  //  G4cout<<"U = "<<U<<G4endl;
     229  //  G4cout<<"N="<<N<<"  P="<<P<<"  H="<<H<<G4endl;
     230  //  G4cout<<"l+ ="<<TransitionProb1<<"  l- ="<< TransitionProb2
     231  //   <<"  l0 ="<< TransitionProb3<<G4endl;
     232  return TransitionProb1 + TransitionProb2 + TransitionProb3;
    237233}
    238234
    239 
    240 G4Fragment G4PreCompoundTransitions::PerformTransition(const G4Fragment & aFragment)
     235void G4PreCompoundTransitions::PerformTransition(G4Fragment & result)
    241236{
    242   G4Fragment result(aFragment);
    243   G4double ChosenTransition = G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
     237  G4double ChosenTransition =
     238    G4UniformRand()*(TransitionProb1 + TransitionProb2 + TransitionProb3);
    244239  G4int deltaN = 0;
    245   G4int Nexcitons = result.GetNumberOfExcitons();
     240  //  G4int Nexcitons = result.GetNumberOfExcitons();
     241  G4int Npart     = result.GetNumberOfParticles();
     242  G4int Ncharged  = result.GetNumberOfCharged();
     243  G4int Nholes    = result.GetNumberOfHoles();
    246244  if (ChosenTransition <= TransitionProb1)
    247245    {
    … …  
    255253    }
    256254
    257   // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and in proportion
    258   // to the number charges w.r.t. number of particles,  PROVIDED that there are charged particles
    259   if(deltaN < 0 && G4UniformRand() <=
    260      static_cast<G4double>(result.GetNumberOfCharged())/static_cast<G4double>(result.GetNumberOfParticles())
    261      && (result.GetNumberOfCharged() >= 1)) {
    262     result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2); // deltaN is negative!
    263   }
     255  // AH/JMQ: Randomly decrease the number of charges if deltaN is -2 and 
     256  // in proportion to the number charges w.r.t. number of particles, 
     257  // PROVIDED that there are charged particles
     258  deltaN /= 2;
     259
     260  //G4cout << "deltaN= " << deltaN << G4endl;
    264261
    265262  // JMQ the following lines have to be before SetNumberOfCharged, otherwise the check on
    266263  // number of charged vs. number of particles fails
    267   result.SetNumberOfParticles(result.GetNumberOfParticles()+deltaN/2);
    268   result.SetNumberOfHoles(result.GetNumberOfHoles()+deltaN/2);
    269 
    270   // With weight Z/A, number of charged particles is increased with +1
    271   if ( ( deltaN > 0 ) &&
    272       (G4UniformRand() <= static_cast<G4double>(result.GetZ()-result.GetNumberOfCharged())/
    273        std::max(static_cast<G4double>(result.GetA()-Nexcitons),1.)))
    274     {
    275       result.SetNumberOfCharged(result.GetNumberOfCharged()+deltaN/2);
    276     }
     264  result.SetNumberOfParticles(Npart+deltaN);
     265  result.SetNumberOfHoles(Nholes+deltaN);
     266
     267  if(deltaN < 0) {
     268    if( Ncharged >= 1 && G4int(Npart*G4UniformRand()) <= Ncharged)
     269      {
     270        result.SetNumberOfCharged(Ncharged+deltaN); // deltaN is negative!
     271      }
     272
     273  } else if ( deltaN > 0 ) {
     274    // With weight Z/A, number of charged particles is increased with +1
     275    G4int A = result.GetA_asInt();
     276    G4int Z = result.GetZ_asInt();
     277    if( G4int(std::max(1, A - Npart)*G4UniformRand()) <= Z)
     278      {
     279        result.SetNumberOfCharged(Ncharged+deltaN);
     280      }
     281  }
    277282 
    278283  // Number of charged can not be greater that number of particles
    279   if ( result.GetNumberOfParticles() < result.GetNumberOfCharged() )
     284  if ( Npart < Ncharged )
    280285    {
    281       result.SetNumberOfCharged(result.GetNumberOfParticles());
    282     }
    283  
    284   return result;
     286      result.SetNumberOfCharged(Npart);
     287    }
     288  //G4cout << "### After transition" << G4endl;
     289  //G4cout << result << G4endl;
    285290}
    286291
Note: See TracChangeset for help on using the changeset viewer.