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/G4PreCompoundEmission.cc

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4PreCompoundEmission.cc,v 1.28 2010/02/25 10:27:36 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4PreCompoundEmission.cc,v 1.32 2010/09/01 15:11:10 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-ref-09 $
    2928//
    3029// -------------------------------------------------------------------
    … …  
    4140// 19.01.2010 V.Ivanchenko simplified computation of parameter an, sample cosTheta
    4241//                         instead of theta; protect all calls to sqrt
     42// 20.08.2010 V.Ivanchenko added G4Pow and G4PreCompoundParameters pointers
     43//                         use int Z and A and cleanup
    4344//
    4445
    4546#include "G4PreCompoundEmission.hh"
    4647#include "G4PreCompoundParameters.hh"
    47 
    4848#include "G4PreCompoundEmissionFactory.hh"
    4949#include "G4HETCEmissionFactory.hh"
    50 
    51 const G4PreCompoundEmission &
    52 G4PreCompoundEmission::operator=(const G4PreCompoundEmission &)
    53 {
    54   throw G4HadronicException(__FILE__, __LINE__,
    55                             "G4PreCompoundEmission::operator= meant to not be accessable");
    56   return *this;
    57 }
    58 
    59 
    60 G4bool G4PreCompoundEmission::operator==(const G4PreCompoundEmission &) const
    61 {
    62   return false;
    63 }
    64 
    65 G4bool G4PreCompoundEmission::operator!=(const G4PreCompoundEmission &) const
    66 {
    67   return true;
    68 }
     50#include "G4HadronicException.hh"
     51#include "G4Pow.hh"
     52#include "Randomize.hh"
    6953
    7054G4PreCompoundEmission::G4PreCompoundEmission()
    7155{
    7256  theFragmentsFactory = new G4PreCompoundEmissionFactory();
    73   theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
     57  theFragmentsVector =
     58    new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
     59  g4pow = G4Pow::GetInstance();
     60  theParameters = G4PreCompoundParameters::GetAddress();
    7461}
    7562
    7663G4PreCompoundEmission::~G4PreCompoundEmission()
    7764{
    78   if (theFragmentsFactory) delete theFragmentsFactory;
    79   if (theFragmentsVector) delete theFragmentsVector;
     65  if (theFragmentsFactory) { delete theFragmentsFactory; }
     66  if (theFragmentsVector)  { delete theFragmentsVector; }
    8067}
    8168
    8269void G4PreCompoundEmission::SetDefaultModel()
    8370{
    84   if (theFragmentsFactory) delete theFragmentsFactory;
     71  if (theFragmentsFactory) { delete theFragmentsFactory; }
    8572  theFragmentsFactory = new G4PreCompoundEmissionFactory();
    8673  if (theFragmentsVector)
    … …  
    9077  else
    9178    {
    92       theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
    93     }
    94   theFragmentsVector->ResetStage();
     79      theFragmentsVector =
     80        new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
     81    }
    9582  return;
    9683}
    … …  
    10693  else
    10794    {
    108       theFragmentsVector = new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
    109     }
    110   theFragmentsVector->ResetStage();
     95      theFragmentsVector =
     96        new G4PreCompoundFragmentVector(theFragmentsFactory->GetFragmentVector());
     97    }
    11198  return;
    11299}
    … …  
    114101G4ReactionProduct * G4PreCompoundEmission::PerformEmission(G4Fragment & aFragment)
    115102{
    116 #ifdef debug
    117   G4Fragment InitialState(aFragment);
    118 #endif
    119103  // Choose a Fragment for emission
    120   G4VPreCompoundFragment * theFragment = theFragmentsVector->ChooseFragment();
    121   if (theFragment == 0)
    122     {
    123       G4cerr <<  "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n"
     104  G4VPreCompoundFragment * thePreFragment = theFragmentsVector->ChooseFragment();
     105  if (thePreFragment == 0)
     106    {
     107      G4cout <<  "G4PreCompoundEmission::PerformEmission : I couldn't choose a fragment\n"
    124108             << "while trying to de-excite\n"
    125              << aFragment << '\n';
     109             << aFragment << G4endl;
    126110      throw G4HadronicException(__FILE__, __LINE__, "");
    127111    }
    128112  // Kinetic Energy of emitted fragment
    129   G4double KineticEnergyOfEmittedFragment = theFragment->GetKineticEnergy(aFragment);
     113  G4double kinEnergyOfEmittedFragment = thePreFragment->GetKineticEnergy(aFragment);
     114  if(kinEnergyOfEmittedFragment < 0.0) { kinEnergyOfEmittedFragment = 0.0; }
    130115 
    131116  // Calculate the fragment momentum (three vector)
    132   G4ThreeVector momentum = AngularDistribution(theFragment,aFragment,KineticEnergyOfEmittedFragment);
     117  AngularDistribution(thePreFragment,aFragment,kinEnergyOfEmittedFragment);
    133118 
    134119  // Mass of emittef fragment
    135   G4double EmittedMass = theFragment->GetNuclearMass();
     120  G4double EmittedMass = thePreFragment->GetNuclearMass();
    136121 
    137122  // Now we can calculate the four momentum
    138123  // both options are valid and give the same result but 2nd one is faster
    139   // G4LorentzVector EmittedMomentum(momentum,std::sqrt(momentum.mag2()+EmittedMass*EmittedMass));
    140   G4LorentzVector EmittedMomentum(momentum,EmittedMass+KineticEnergyOfEmittedFragment);
     124  G4LorentzVector Emitted4Momentum(theFinalMomentum,
     125                                   EmittedMass + kinEnergyOfEmittedFragment);
    141126   
    142127  // Perform Lorentz boost
    143   EmittedMomentum.boost(aFragment.GetMomentum().boostVector()); 
     128  G4LorentzVector Rest4Momentum = aFragment.GetMomentum();
     129  Emitted4Momentum.boost(Rest4Momentum.boostVector()); 
    144130
    145131  // Set emitted fragment momentum
    146   theFragment->SetMomentum(EmittedMomentum);   
    147 
     132  thePreFragment->SetMomentum(Emitted4Momentum);       
    148133
    149134  // NOW THE RESIDUAL NUCLEUS
    150135  // ------------------------
    151    
    152   // Now the residual nucleus.
    153   // The energy conservation says that
    154   G4double ResidualEcm =
    155     //    aFragment.GetGroundStateMass() + aFragment.GetExcitationEnergy() // initial energy in cm
    156     aFragment.GetMomentum().m()
    157     - (EmittedMass+KineticEnergyOfEmittedFragment);
    158 
    159   // Then the four momentum for residual is
    160   G4LorentzVector RestMomentum(-momentum,ResidualEcm);
    161   // This could save a Lorentz boost
    162   // G4LorentzVector RestMomentum2(aFragment.GetMomentum()-EmittedMomentum);
    163 
    164   // Just for test
    165   // Excitation energy
    166   //  G4double anU = ResidualEcm - theFragment->GetRestNuclearMass();
    167   // This is equivalent
    168   //  G4double anU = theFragment->GetMaximalKineticEnergy() - KineticEnergyOfEmittedFragment +
    169   //    theFragment->GetCoulombBarrier();
    170    
    171   // check that Excitation energy is >= 0
    172   G4double anU = RestMomentum.m()-theFragment->GetRestNuclearMass();
    173   if (anU < 0.0) {
    174     throw G4HadronicException(__FILE__, __LINE__,
    175                               "G4PreCompoundModel::DeExcite: Excitation energy less than 0!");
    176   }
    177      
     136
     137  Rest4Momentum -= Emitted4Momentum;
     138   
    178139  // Update nucleus parameters:
    179140  // --------------------------
    … …  
    181142  // Number of excitons
    182143  aFragment.SetNumberOfParticles(aFragment.GetNumberOfParticles()-
    183                                  static_cast<G4int>(theFragment->GetA()));
     144                                 thePreFragment->GetA());
    184145  // Number of charges
    185146  aFragment.SetNumberOfCharged(aFragment.GetNumberOfCharged()-
    186                                static_cast<G4int>(theFragment->GetZ()));
    187    
    188   // Atomic number
    189   aFragment.SetA(theFragment->GetRestA());
    190    
    191   // Charge
    192   aFragment.SetZ(theFragment->GetRestZ());
    193 
    194    
    195   // Perform Lorentz boosts
    196   RestMomentum.boost(aFragment.GetMomentum().boostVector());
    197 
    198   // Update nucleus momentum
    199   aFragment.SetMomentum(RestMomentum);
     147                               thePreFragment->GetZ());
     148   
     149  // Z and A
     150  aFragment.SetZandA_asInt(thePreFragment->GetRestZ(),
     151                           thePreFragment->GetRestA());
     152   
     153  // Update nucleus momentum
     154  // A check on consistence of Z, A, and mass will be performed
     155  aFragment.SetMomentum(Rest4Momentum);
    200156       
    201157  // Create a G4ReactionProduct
    202   G4ReactionProduct * MyRP = theFragment->GetReactionProduct();
    203 #ifdef PRECOMPOUND_TEST
    204   MyRP->SetCreatorModel("G4PreCompoundModel");
    205 #endif
    206 #ifdef debug
    207   CheckConservation(InitialState,aFragment,MyRP);
    208 #endif
     158  G4ReactionProduct * MyRP = thePreFragment->GetReactionProduct();
     159
     160  //G4cout << "G4PreCompoundEmission::Fragment emitted" << G4endl;
     161  //G4cout << thePreFragment << G4endl;
     162
    209163  return MyRP;
    210164}
    211165
    212 G4ThreeVector
    213 G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment * theFragment,
     166void
     167G4PreCompoundEmission::AngularDistribution(G4VPreCompoundFragment* thePreFragment,
    214168                                           const G4Fragment& aFragment,
    215                                            const G4double kinEnergyOfEmittedFrag) const
    216 {
    217   G4double p = aFragment.GetNumberOfParticles();
    218   G4double h = aFragment.GetNumberOfHoles();
     169                                           G4double ekin)
     170{
     171  G4int p = aFragment.GetNumberOfParticles();
     172  G4int h = aFragment.GetNumberOfHoles();
    219173  G4double U = aFragment.GetExcitationEnergy();
    220174
    221   G4double ekin = std::max(0.0, kinEnergyOfEmittedFrag);
    222        
    223175  // Emission particle separation energy
    224   G4double Bemission = theFragment->GetBindingEnergy();
     176  G4double Bemission = thePreFragment->GetBindingEnergy();
    225177
    226178  // Fermi energy
    227   G4double Ef = G4PreCompoundParameters::GetAddress()->GetFermiEnergy();
     179  G4double Ef = theParameters->GetFermiEnergy();
    228180       
    229181  //
    … …  
    231183  //  G4double g = (6.0/pi2)*aFragment.GetA()*
    232184
    233   G4double g = (6.0/pi2)*aFragment.GetA()
    234     *G4PreCompoundParameters::GetAddress()->GetLevelDensity();
     185  G4double g = (6.0/pi2)*aFragment.GetA_asInt()*theParameters->GetLevelDensity();
    235186       
    236187  // Average exciton energy relative to bottom of nuclear well
    237   G4double Eav = 2.0*p*(p+1.0)/((p+h)*g);
     188  G4double Eav = 2*p*(p+1)/((p+h)*g);
    238189       
    239190  // Excitation energy relative to the Fermi Level
    … …  
    241192  //  G4double Uf = U - KineticEnergyOfEmittedFragment - Bemission;
    242193
    243   G4double w_num = rho(p+1,h,g,Uf,Ef);
    244   G4double w_den = rho(p,h,g,Uf,Ef);
     194  G4double w_num = rho(p+1, h, g, Uf, Ef);
     195  G4double w_den = rho(p,   h, g, Uf, Ef);
    245196  if (w_num > 0.0 && w_den > 0.0)
    246197    {
    … …  
    253204    }
    254205 
    255 
    256206  // VI + JMQ 19/01/2010 update computation of the parameter an
    257207  //
    … …  
    262212    G4double zeta = std::max(1.0,9.3/std::sqrt(ekin/MeV));
    263213 
    264     an = 3.0*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
     214    // This should be the projectile energy. If I would know which is
     215    // the projectile (proton, neutron) I could remove the binding energy.
     216    // But, what happens if INC precedes precompound? This approximation
     217    // seems to work well enough
     218    G4double ProjEnergy = aFragment.GetExcitationEnergy();
     219
     220    an = 3*std::sqrt((ProjEnergy+Ef)*Eeff)/(zeta*Eav);
    265221
    266222    G4int ne = aFragment.GetNumberOfExcitons() - 1;
    … …  
    283239  } 
    284240
    285   G4double phi = twopi*G4UniformRand();
     241  G4double phi = CLHEP::twopi*G4UniformRand();
    286242 
    287243  // Calculate the momentum magnitude of emitted fragment       
    288   G4double pmag = std::sqrt(ekin*(ekin + 2.0*theFragment->GetNuclearMass()));
     244  G4double pmag = std::sqrt(ekin*(ekin + 2.0*thePreFragment->GetNuclearMass()));
    289245 
    290246  G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
    291247
    292   G4ThreeVector momentum(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
     248  theFinalMomentum.set(pmag*std::cos(phi)*sint,pmag*std::sin(phi)*sint,pmag*cost);
     249
    293250  // theta is the angle wrt the incident direction
    294   momentum.rotateUz(theIncidentDirection);
    295 
    296   return momentum;
    297 }
    298 
    299 G4double G4PreCompoundEmission::rho(const G4double p, const G4double h, const G4double g,
    300                                     const G4double E, const G4double Ef) const
     251  G4ThreeVector theIncidentDirection = aFragment.GetMomentum().vect().unit();
     252  theFinalMomentum.rotateUz(theIncidentDirection);
     253}
     254
     255G4double G4PreCompoundEmission::rho(G4int p, G4int h, G4double g,
     256                                    G4double E, G4double Ef) const
    301257{       
    302258  // 25.02.2010 V.Ivanchenko added more protections
    … …  
    306262  if ( E - Aph < 0.0) { return 0.0; }
    307263 
    308   G4double logConst =  (p+h)*std::log(g) - logfactorial(p+h-1) - logfactorial(p) - logfactorial(h);
     264  G4double logConst =  (p+h)*std::log(g)
     265    - g4pow->logfactorial(p+h-1) - g4pow->logfactorial(p) - g4pow->logfactorial(h);
    309266
    310267  // initialise values using j=0
    … …  
    312269  G4double t1=1;
    313270  G4double t2=1;
    314   G4double logt3=(p+h-1) * std::log(E-Aph) + logConst;
     271  G4double logt3 = (p+h-1) * std::log(E-Aph) + logConst;
    315272  const G4double logmax = 200.;
    316273  if(logt3 > logmax) { logt3 = logmax; }
    … …  
    333290  return tot;
    334291}
    335 
    336 G4double G4PreCompoundEmission::factorial(G4double a) const
    337 {
    338   // Values of factorial function from 0 to 60
    339   const G4int factablesize = 61;
    340   static const G4double fact[factablesize] =
    341     {
    342       1.0, // 0!
    343       1.0, // 1!
    344       2.0, // 2!
    345       6.0, // 3!
    346       24.0, // 4!
    347       120.0, // 5!
    348       720.0, // 6!
    349       5040.0, // 7!
    350       40320.0, // 8!
    351       362880.0, // 9!
    352       3628800.0, // 10!
    353       39916800.0, // 11!
    354       479001600.0, // 12!
    355       6227020800.0, // 13!
    356       87178291200.0, // 14!
    357       1307674368000.0, // 15!
    358       20922789888000.0, // 16!
    359       355687428096000.0, // 17!
    360       6402373705728000.0, // 18!
    361       121645100408832000.0, // 19!
    362       2432902008176640000.0, // 20!
    363       51090942171709440000.0, // 21!
    364       1124000727777607680000.0, // 22!
    365       25852016738884976640000.0, // 23!
    366       620448401733239439360000.0, // 24!
    367       15511210043330985984000000.0, // 25!
    368       403291461126605635584000000.0, // 26!
    369       10888869450418352160768000000.0, // 27!
    370       304888344611713860501504000000.0, // 28!
    371       8841761993739701954543616000000.0, // 29!
    372       265252859812191058636308480000000.0, // 30!
    373       8222838654177922817725562880000000.0, // 31!
    374       263130836933693530167218012160000000.0, // 32!
    375       8683317618811886495518194401280000000.0, // 33!
    376       295232799039604140847618609643520000000.0, // 34!
    377       10333147966386144929666651337523200000000.0, // 35!
    378       371993326789901217467999448150835200000000.0, // 36!
    379       13763753091226345046315979581580902400000000.0, // 37!
    380       523022617466601111760007224100074291200000000.0, // 38!
    381       20397882081197443358640281739902897356800000000.0, // 39!
    382       815915283247897734345611269596115894272000000000.0, // 40!
    383       33452526613163807108170062053440751665152000000000.0, // 41!
    384       1405006117752879898543142606244511569936384000000000.0, // 42!
    385       60415263063373835637355132068513997507264512000000000.0, // 43!
    386       2658271574788448768043625811014615890319638528000000000.0, // 44!
    387       119622220865480194561963161495657715064383733760000000000.0, // 45!
    388       5502622159812088949850305428800254892961651752960000000000.0, // 46!
    389       258623241511168180642964355153611979969197632389120000000000.0, // 47!
    390       12413915592536072670862289047373375038521486354677760000000000.0, // 48!
    391       608281864034267560872252163321295376887552831379210240000000000.0, // 49!
    392       30414093201713378043612608166064768844377641568960512000000000000.0, // 50!
    393       1551118753287382280224243016469303211063259720016986112000000000000.0, // 51!
    394       80658175170943878571660636856403766975289505440883277824000000000000.0, // 52!
    395       4274883284060025564298013753389399649690343788366813724672000000000000.0, // 53!
    396       230843697339241380472092742683027581083278564571807941132288000000000000.0, // 54!
    397       12696403353658275925965100847566516959580321051449436762275840000000000000.0, // 55!
    398       710998587804863451854045647463724949736497978881168458687447040000000000000.0, // 56!
    399       40526919504877216755680601905432322134980384796226602145184481280000000000000.0, // 57!
    400       2350561331282878571829474910515074683828862318181142924420699914240000000000000.0, // 58!
    401       138683118545689835737939019720389406345902876772687432540821294940160000000000000.0, // 59!
    402       8320987112741390144276341183223364380754172606361245952449277696409600000000000000.0  // 60!
    403     };
    404   //    fact[0] = 1;
    405   //    for (G4int n = 1; n < 21; n++) {
    406   //      fact[n] = fact[n-1]*static_cast<G4double>(n);
    407   //    }
    408   G4double result(0.0);
    409   G4int ia = static_cast<G4int>(a);
    410   if (ia < factablesize)
    411     {
    412       result = fact[ia];
    413     }
    414   else
    415     {
    416       result = fact[factablesize-1];
    417       for (G4int n = factablesize; n < ia+1; ++n)
    418         {
    419           result *= static_cast<G4double>(n);
    420         }
    421     }
    422    
    423     return result;
    424 }
    425 G4double G4PreCompoundEmission::logfactorial(G4double a) const
    426 {
    427   // Values of logs of factorial function from 0 to 60
    428 
    429   G4double result(0.0);
    430   const G4int factablesize = 61;
    431   const G4double halfLn2pi = 0.918938533;      // 0.5 log(2* pi)
    432   static G4double logfact[factablesize];
    433   static bool needinit=true;
    434  
    435   if (needinit)
    436   {
    437       needinit=false;
    438       for ( G4int n=0; n < factablesize; ++n)
    439       {
    440          logfact[n]=std::log(factorial(n));
    441       }
    442   }
    443 
    444   G4int ia = static_cast<G4int>(a);
    445   if (ia < factablesize)
    446   {
    447       result = logfact[ia];
    448   } else {
    449       result = (ia+0.5)*std::log(G4double(ia)) - ia + halfLn2pi;
    450   }
    451    
    452   return result;
    453 }
    454 
    455 #ifdef debug
    456 void G4PreCompoundEmission::CheckConservation(const G4Fragment & theInitialState,
    457                                               const G4Fragment & theResidual,
    458                                               G4ReactionProduct * theEmitted) const
    459 {
    460   G4double ProductsEnergy = theEmitted->GetTotalEnergy() + theResidual.GetMomentum().e();
    461   G4ThreeVector ProductsMomentum(theEmitted->GetMomentum()+theResidual.GetMomentum().vect());
    462   G4int ProductsA = theEmitted->GetDefinition()->GetBaryonNumber() + theResidual.GetA();
    463   G4int ProductsZ = theEmitted->GetDefinition()->GetPDGCharge() + theResidual.GetZ();
    464 
    465   if (ProductsA != theInitialState.GetA()) {
    466     G4cout << "!!!!!!!!!! Baryonic Number Conservation Violation !!!!!!!!!!" << G4endl;
    467     G4cout << "G4PreCompoundEmission.cc: Barionic Number Conservation"
    468            << G4endl;
    469     G4cout << "Initial A = " << theInitialState.GetA()
    470            << "   Fragments A = " << ProductsA << "   Diference --> "
    471            << theInitialState.GetA() - ProductsA << G4endl;
    472   }
    473   if (ProductsZ != theInitialState.GetZ()) {
    474     G4cout << "!!!!!!!!!! Charge Conservation Violation !!!!!!!!!!" << G4endl;
    475     G4cout << "G4PreCompoundEmission.cc: Charge Conservation test"
    476            << G4endl;
    477     G4cout << "Initial Z = " << theInitialState.GetZ()
    478            << "   Fragments Z = " << ProductsZ << "   Diference --> "
    479            << theInitialState.GetZ() - ProductsZ << G4endl;
    480   }
    481   if (std::abs(ProductsEnergy-theInitialState.GetMomentum().e()) > 10.0*eV) {
    482     G4cout << "!!!!!!!!!! Energy Conservation Violation !!!!!!!!!!" << G4endl;
    483     G4cout << "G4PreCompoundEmission.cc: Energy Conservation test"
    484            << G4endl;
    485     G4cout << "Initial E = " << theInitialState.GetMomentum().e()/MeV << " MeV"
    486            << "   Fragments E = " << ProductsEnergy/MeV  << " MeV   Diference --> "
    487            << (theInitialState.GetMomentum().e() - ProductsEnergy)/MeV << " MeV" << G4endl;
    488   }
    489   if (std::abs(ProductsMomentum.x()-theInitialState.GetMomentum().x()) > 10.0*eV ||
    490       std::abs(ProductsMomentum.y()-theInitialState.GetMomentum().y()) > 10.0*eV ||
    491       std::abs(ProductsMomentum.z()-theInitialState.GetMomentum().z()) > 10.0*eV) {
    492     G4cout << "!!!!!!!!!! Momentum Conservation Violation !!!!!!!!!!" << G4endl;
    493     G4cout << "G4PreCompoundEmission.cc: Momentum Conservation test"
    494            << G4endl;
    495     G4cout << "Initial P = " << theInitialState.GetMomentum().vect() << " MeV"
    496            << "   Fragments P = " << ProductsMomentum  << " MeV   Diference --> "
    497            << theInitialState.GetMomentum().vect() - ProductsMomentum << " MeV" << G4endl;
    498   }
    499   return;
    500 }
    501 
    502 #endif
Note: See TracChangeset for help on using the changeset viewer.