Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

Location:
trunk/source/processes/hadronic/models/de_excitation/evaporation/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4AlphaEvaporationProbability.cc

    r1055 r1315  
    225225  ecut = (ecut-a) / (p+p);
    226226  ecut2 = ecut;
    227   if (cut < 0.) ecut2 = ecut - 2.;
     227//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     228//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     229//  if (cut < 0.) ecut2 = ecut - 2.;
     230  if (cut < 0.) ecut2 = ecut;
    228231  elab = K * FragmentA / ResidualA;
    229232  sig = 0.;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4DeuteronEvaporationProbability.cc

    r1055 r1315  
    206206       
    207207  //JMQ 13/02/09 increase of reduced radius to lower the barrier
    208   // ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
    209   ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     208  //  ec = 1.44 * theZ * ResidualZ / (1.5*ResidualAthrd+ra);
     209   ec = 1.44 * theZ * ResidualZ / (1.7*ResidualAthrd+ra);
     210
    210211  ecsq = ec * ec;
    211212  p = p0 + p1/ec + p2/ecsq;
    … …  
    225226  ecut = (ecut-a) / (p+p);
    226227  ecut2 = ecut;
    227   if (cut < 0.) ecut2 = ecut - 2.;
     228//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     229//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     230//  if (cut < 0.) ecut2 = ecut - 2.;
     231  if (cut < 0.) ecut2 = ecut;
    228232  elab = K * FragmentA / ResidualA;
    229233  sig = 0.;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4Evaporation.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4Evaporation.cc,v 1.15 2009/09/16 15:32:25 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4Evaporation.cc,v 1.25 2010/06/09 11:56:47 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3838// superimposed Coulomb barrier (if useSICBis set true, by default is false)
    3939//
    40 // V.Ivanchenko added G4EvaporationDefaultGEMFactory option, 27/07/09
     40// V.Ivanchenko (27 July 2009) added G4EvaporationDefaultGEMFactory option
     41// V.Ivanchenko (10 May 2010) rewrited BreakItUp method: do not make new/delete
     42//                            photon channel first, fission second,
     43//                            added G4UnstableFragmentBreakUp to decay
     44//                            unphysical fragments (like 2n or 2p),
     45//                            use Z and A integer
    4146
    4247#include "G4Evaporation.hh"
    … …  
    4550#include "G4EvaporationDefaultGEMFactory.hh"
    4651#include "G4HadronicException.hh"
    47 #include <numeric>
     52#include "G4NistManager.hh"
    4853
    4954G4Evaporation::G4Evaporation()
    5055{
    51   theChannelFactory = new G4EvaporationFactory();
    52   //theChannelFactory = new G4EvaporationDefaultGEMFactory();
    53   theChannels = theChannelFactory->GetChannel();
     56  //theChannelFactory = new G4EvaporationFactory();
     57  theChannelFactory = new G4EvaporationDefaultGEMFactory();
     58  InitialiseEvaporation();
     59}
     60
     61G4Evaporation::G4Evaporation(std::vector<G4VEvaporationChannel*> * aChannelsVector)
     62  : theChannels(aChannelsVector), theChannelFactory(0), nChannels(0)
     63{
     64  InitialiseEvaporation();
    5465}
    5566
    5667G4Evaporation::~G4Evaporation()
    5768{
    58   if (theChannels != 0) theChannels = 0;
    59   if (theChannelFactory != 0) delete theChannelFactory;
     69  if (theChannels != 0) { theChannels = 0; }
     70  if (theChannelFactory != 0) { delete theChannelFactory; }
     71}
     72
     73void G4Evaporation::InitialiseEvaporation()
     74{
     75  nist = G4NistManager::Instance();
     76  minExcitation = CLHEP::keV;
     77  if(theChannelFactory) { theChannels = theChannelFactory->GetChannel(); }
     78  nChannels = theChannels->size();
     79  probabilities.resize(nChannels, 0.0);
     80  Initialise();
     81}
     82
     83void G4Evaporation::Initialise()
     84{
     85  // loop over evaporation channels
     86  std::vector<G4VEvaporationChannel*>::iterator i;
     87  for (i=theChannels->begin(); i != theChannels->end(); ++i)
     88    {
     89      // for inverse cross section choice
     90      (*i)->SetOPTxs(OPTxs);
     91      // for superimposed Coulomb Barrier for inverse cross sections
     92      (*i)->UseSICB(useSICB);
     93    }
    6094}
    6195
    … …  
    6498  if (theChannelFactory != 0) delete theChannelFactory;
    6599  theChannelFactory = new G4EvaporationFactory();
    66   theChannels = theChannelFactory->GetChannel();
     100  InitialiseEvaporation();
    67101}
    68102
    … …  
    71105  if (theChannelFactory != 0) delete theChannelFactory;
    72106  theChannelFactory = new G4EvaporationGEMFactory();
    73   theChannels = theChannelFactory->GetChannel();
     107  InitialiseEvaporation();
    74108}
    75109
    … …  
    78112  if (theChannelFactory != 0) delete theChannelFactory;
    79113  theChannelFactory = new G4EvaporationDefaultGEMFactory();
    80   theChannels = theChannelFactory->GetChannel();
    81 }
    82 
     114  InitialiseEvaporation();
     115}
    83116
    84117G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
    85118{
    86     G4FragmentVector * theResult = new G4FragmentVector;
     119  G4FragmentVector * theResult = new G4FragmentVector;
     120  G4FragmentVector * theTempResult;
     121
     122  // The residual nucleus (after evaporation of each fragment)
     123  G4Fragment* theResidualNucleus = new G4Fragment(theNucleus);
     124
     125  G4double totprob, prob, oldprob = 0.0;
     126  G4int maxchannel, i;
     127
     128  G4int Amax = theResidualNucleus->GetA_asInt();
     129
     130  // Starts loop over evaporated particles, loop is limited by number
     131  // of nucleons
     132  for(G4int ia=0; ia<Amax; ++ia) {
     133 
     134    // g,n,p - evaporation is finished
     135    G4int A = theResidualNucleus->GetA_asInt();
     136    if(1 >= A) {
     137      theResult->push_back(theResidualNucleus);
     138      return theResult;
     139    }
     140
     141    // check if it is stable, then finish evaporation
     142    G4int Z = theResidualNucleus->GetZ_asInt();
     143    G4double abun = nist->GetIsotopeAbundance(Z, A);
     144    /*
     145    G4cout << "### G4Evaporation::BreakItUp step " << ia << " Z= " << Z
     146           << " A= " << A << " Eex(MeV)= "
     147           << theResidualNucleus->GetExcitationEnergy()
     148           << " aban= " << abun << G4endl;
     149    */
     150    if(theResidualNucleus->GetExcitationEnergy() <= minExcitation &&
     151       (abun > 0.0))
     152      {
     153        theResult->push_back(theResidualNucleus);
     154        return theResult;
     155      }
     156
     157    totprob = 0.0;
     158    maxchannel = nChannels;
     159
     160    //G4cout << "### Evaporation loop #" << ia
     161    //     << "  Fragment: " << theResidualNucleus << G4endl;
     162
     163    // loop over evaporation channels
     164    for(i=0; i<nChannels; ++i) {
     165      (*theChannels)[i]->Initialize(*theResidualNucleus);
     166      prob = (*theChannels)[i]->GetEmissionProbability();
     167      //G4cout << "  Channel# " << i << "  prob= " << prob << G4endl;
     168
     169      //if(0 == i && 0.0 == abun) { prob = 0.0; }
     170
     171      totprob += prob;
     172      probabilities[i] = totprob;
     173      // if two recent probabilities are near zero stop computations
     174      if(i>=8) {
     175        if(prob <= totprob*1.e-8 && oldprob <= totprob*1.e-8) {
     176          maxchannel = i+1;
     177          break;
     178        }
     179      }
     180      oldprob = prob;
     181    }
     182
     183    // photon evaporation in the case of no other channels available
     184    // do evaporation chain and reset total probability
     185    if(0.0 < totprob && probabilities[0] == totprob) {
     186      theTempResult = (*theChannels)[0]->BreakUpFragment(theResidualNucleus);
     187      if(theTempResult) {
     188        size_t nsec = theTempResult->size();
     189        for(size_t j=0; j<nsec; ++j) {
     190          theResult->push_back((*theTempResult)[j]);
     191        }
     192        delete theTempResult;
     193      }
     194      totprob = 0.0;
     195    }
     196
     197    // stable fragnent - evaporation is finished
     198    if(0.0 == totprob) {
     199
     200      // if fragment is exotic, then try to decay it
     201      if(0.0 == abun && Z < 20) {
     202        //G4cout << "$$$ Decay exotic fragment" << G4endl;
     203        theTempResult = unstableBreakUp.BreakUpFragment(theResidualNucleus);
     204        if(theTempResult) {
     205          size_t nsec = theTempResult->size();
     206          for(size_t j=0; j<nsec; ++j) {
     207            theResult->push_back((*theTempResult)[j]);
     208          }
     209          delete theTempResult;
     210        }
     211      }
     212
     213      // save residual fragment
     214      theResult->push_back(theResidualNucleus);
     215      return theResult;
     216    }
     217
     218
     219    // select channel
     220    totprob *= G4UniformRand();
     221    // loop over evaporation channels
     222    for(i=0; i<maxchannel; ++i) { if(probabilities[i] >= totprob) { break; } }
     223
     224    // this should not happen
     225    if(i >= nChannels) { i = nChannels - 1; }
     226
     227
     228    // single photon evaporation, primary pointer is kept
     229    if(0 == i) {
     230      G4Fragment* gamma = (*theChannels)[0]->EmittedFragment(theResidualNucleus);
     231      if(gamma) { theResult->push_back(gamma); }
     232
     233      // fission, return results to the main loop if fission is succesful
     234    } else if(1 == i) {
     235      theTempResult = (*theChannels)[1]->BreakUp(*theResidualNucleus);
     236      if(theTempResult) {
     237        size_t nsec = theTempResult->size();
     238        G4bool deletePrimary = true;
     239        for(size_t j=0; j<nsec; ++j) {
     240          if(theResidualNucleus == (*theTempResult)[j]) { deletePrimary = false; }
     241          theResult->push_back((*theTempResult)[j]);
     242        }
     243        if(deletePrimary) { delete theResidualNucleus; }
     244        delete theTempResult;
     245        return theResult;
     246      }
     247
     248      // other channels
     249    } else {
     250      theTempResult = (*theChannels)[i]->BreakUp(*theResidualNucleus);
     251      if(theTempResult) {
     252        size_t nsec = theTempResult->size();
     253        if(nsec > 0) {
     254          --nsec;
     255          for(size_t j=0; j<nsec; ++j) {
     256            theResult->push_back((*theTempResult)[j]);
     257          }
     258          // if the residual change its pointer
     259          // then delete previous residual fragment and update to the new
     260          if(theResidualNucleus != (*theTempResult)[nsec] ) {
     261            delete theResidualNucleus;
     262            theResidualNucleus = (*theTempResult)[nsec];
     263          }
     264        }
     265        delete theTempResult;
     266      }
     267    }
     268  }
     269
     270  // loop is stopped, save residual
     271  theResult->push_back(theResidualNucleus);
     272 
     273#ifdef debug
     274  G4cout << "======== Evaporation Conservation Test ===========\n"
     275         << "==================================================\n";
     276  CheckConservation(theNucleus,theResult);
     277  G4cout << "==================================================\n";
     278#endif
     279  return theResult;
     280}
     281
     282/*
     283G4FragmentVector * G4Evaporation::BreakItUp(const G4Fragment &theNucleus)
     284{
     285  G4FragmentVector * theResult = new G4FragmentVector;
    87286
    88287    // CHECK that Excitation Energy != 0
    … …  
    233432    return theResult;
    234433}
    235 
     434*/
    236435
    237436

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationDefaultGEMFactory.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4EvaporationDefaultGEMFactory.cc,v 1.1 2009/07/27 10:20:13 vnivanch Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4EvaporationDefaultGEMFactory.cc,v 1.2 2010/04/27 11:43:16 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    122122  theChannel->reserve(68);
    123123
     124  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
     125  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
     126
    124127  // JMQ 220709 standard particle evaporation channels (Z<3,A<5)
    125128  theChannel->push_back( new G4NeutronEvaporationChannel() );  // n
    … …  
    192195  theChannel->push_back( new G4Mg28GEMChannel() );     // Mg28
    193196
    194   theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
    195   theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
    196 
    197197  return theChannel;
    198198

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationFactory.cc

    r1228 r1315  
    2525//
    2626//
    27 // $Id: G4EvaporationFactory.cc,v 1.4 2006/06/29 20:10:29 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     27// $Id: G4EvaporationFactory.cc,v 1.5 2010/04/27 11:43:16 vnivanch Exp $
     28// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2929//
    3030// Hadronic Process: Nuclear De-excitations
    … …  
    3232
    3333#include "G4EvaporationFactory.hh"
    34 
    3534
    3635#include "G4NeutronEvaporationChannel.hh"
    … …  
    4443#include "G4PhotonEvaporation.hh"
    4544
     45G4EvaporationFactory::G4EvaporationFactory()
     46{}
    4647
    47 const G4EvaporationFactory &
    48 G4EvaporationFactory::operator=(const G4EvaporationFactory & )
    49 {
    50   throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationFactory::operator= meant to not be accessable.");
    51   return *this;
    52 }
    53 
    54 G4bool
    55 G4EvaporationFactory::operator==(const G4EvaporationFactory & ) const
    56 {
    57   throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationFactory::operator== meant to not be accessable.");
    58   return false;
    59 }
    60 
    61 G4bool
    62 G4EvaporationFactory::operator!=(const G4EvaporationFactory & ) const
    63 {
    64   throw G4HadronicException(__FILE__, __LINE__, "G4EvaporationFactory::operator!= meant to not be accessable.");
    65   return true;
    66 }
    67 
    68 
     48G4EvaporationFactory::~G4EvaporationFactory()
     49{}
    6950
    7051std::vector<G4VEvaporationChannel*> *
    … …  
    7556  theChannel->reserve(8);
    7657
     58  theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
     59  theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
     60
    7761  theChannel->push_back( new G4NeutronEvaporationChannel() );  // n
    7862  theChannel->push_back( new G4ProtonEvaporationChannel() );   // p
    … …  
    8266  theChannel->push_back( new G4AlphaEvaporationChannel() );    // Alpha
    8367
    84   theChannel->push_back( new G4CompetitiveFission() );         // Fission Channel
    85   theChannel->push_back( new G4PhotonEvaporation() );          // Photon Channel
    8668
    8769  return theChannel;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4EvaporationProbability.cc

    r1055 r1315  
    4141#include "G4EvaporationProbability.hh"
    4242#include "G4PairingCorrection.hh"
    43 
     43#include "G4ParticleTable.hh"
     44#include "G4IonTable.hh"
    4445
    4546

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4He3EvaporationProbability.cc

    r1055 r1315  
    221221  ecut = (ecut-a) / (p+p);
    222222  ecut2 = ecut;
    223   if (cut < 0.) ecut2 = ecut - 2.;
     223//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     224//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     225//  if (cut < 0.) ecut2 = ecut - 2.;
     226  if (cut < 0.) ecut2 = ecut;
    224227  elab = K * FragmentA / ResidualA;
    225228  sig = 0.;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4ProtonEvaporationProbability.cc

    r962 r1315  
    9191
    9292///////////////////////////////////////////////////////////////////////////////////
    93 //J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections
     93//J. M. Quesada (Dec 2007-June 2008): New inverse reaction cross sections for protons
    9494//OPT=0 Dostrovski's parameterization
    95 //OPT=1,2 Chatterjee's paramaterization
    96 //OPT=3,4 Kalbach's parameterization
     95//OPT=1 Chatterjee's parameterization
     96//OPT=2,4 Wellisch's parameterization
     97//OPT=3 Kalbach's parameterization
    9798//
    9899G4double G4ProtonEvaporationProbability::CrossSection(const  G4Fragment & fragment, const  G4double K)
    … …  
    110111  FragmentAthrd=std::pow(FragmentA,0.33333);
    111112  U=fragment.GetExcitationEnergy();
    112 
    113113
    114114  if (OPTxs==0) {std::ostringstream errOs;
    … …  
    223223//     ** p from  becchetti and greenlees (but modified with sub-barrier
    224224//     ** correction function and xp2 changed from -449)
    225 // JMQ (june 2008) : refinement of proton cross section for light systems
    226 //
     225
    227226  G4double landa, landa0, landa1, mu, mu0, mu1,nu, nu0, nu1, nu2;
    228227  G4double p, p0, p1, p2;
    … …  
    279278  ecut = (ecut-a) / (p+p);
    280279  ecut2 = ecut;
    281   if (cut < 0.) ecut2 = ecut - 2.;
     280//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     281//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     282//  if (cut < 0.) ecut2 = ecut - 2.;
     283 if (cut < 0.) ecut2 = ecut;
    282284  elab = K * FragmentA / ResidualA;
    283285  sig = 0.;
    … …  
    286288    signor2 = (ec-elab-c) / w;
    287289    signor2 = 1. + std::exp(signor2);
    288     sig = sig / signor2;
    289     // signor2 is empirical global corr'n at low elab for protons in PRECO, not enough for light targets
    290     //  refinement for proton cross section
    291     if (ResidualZ<=26)
    292       sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));     
    293                        }              //end for E<Ec
     290    sig = sig / signor2;     
     291                       }       //end for E<=Ec
    294292  else            {           //start for  E>Ec
    295293    sig = (landa*elab+mu+nu/elab) * signor;
    296 
    297     //  refinement for proton cross section
    298     if ( ResidualZ<=26 && elab <=(afit*ResidualZ+bfit)*ec)
    299            sig = sig*std::exp(-(a2*ResidualZ + b2)*(elab-(afit*ResidualZ+bfit)*ec)*(elab-(afit*ResidualZ+bfit)*ec));
    300                                                                                
    301 //
    302294    geom = 0.;
    303295

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4TritonEvaporationProbability.cc

    r1055 r1315  
    219219  ecut = (ecut-a) / (p+p);
    220220  ecut2 = ecut;
    221   if (cut < 0.) ecut2 = ecut - 2.;
     221//JMQ 290310 for avoiding unphysical increase below minimum (at ecut)
     222//ecut<0 means that there is no cut with energy axis, i.e. xs is set to 0 bellow minimum
     223//  if (cut < 0.) ecut2 = ecut - 2.;
     224  if (cut < 0.) ecut2 = ecut;
    222225  elab = K * FragmentA / ResidualA;
    223226  sig = 0.;

  • trunk/source/processes/hadronic/models/de_excitation/evaporation/src/G4VEvaporation.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 //
    27 // $Id: G4VEvaporation.cc,v 1.5 2006/06/29 20:10:49 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4VEvaporation.cc,v 1.7 2010/04/27 14:00:40 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2928//
    3029// Hadronic Process: Nuclear De-excitations
    … …  
    3231//
    3332
     33#include "G4VEvaporation.hh"
    3434
    35 #include "G4VEvaporation.hh"
    36 #include "G4HadronicException.hh"
     35G4VEvaporation::G4VEvaporation()
     36{}
    3737
    38 G4VEvaporation::G4VEvaporation(const G4VEvaporation &)
    39 {
    40  throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporation::copy_constructor meant to not be accessable");
    41 }
     38G4VEvaporation::~G4VEvaporation()
     39{}
    4240
    43 
    44 
    45 
    46 const G4VEvaporation & G4VEvaporation::operator=(const G4VEvaporation &)
    47 {
    48   throw G4HadronicException(__FILE__, __LINE__, "G4VEvaporation::operator= meant to not be accessable");
    49   return *this;
    50 }
    51 
    52 
    53 G4bool G4VEvaporation::operator==(const G4VEvaporation &) const
    54 {
    55   return false;
    56 }
    57 
    58 G4bool G4VEvaporation::operator!=(const G4VEvaporation &) const
    59 {
    60   return true;
    61 }
     41void G4VEvaporation::Initialise()
     42{}
    6243
    6344
Note: See TracChangeset for help on using the changeset viewer.