G4BertiniEvaporation Class Reference

#include <G4BertiniEvaporation.hh>

Inheritance diagram for G4BertiniEvaporation:

G4VEvaporation

Public Member Functions

 G4BertiniEvaporation ()
 ~G4BertiniEvaporation ()
virtual G4FragmentVectorBreakItUp (const G4Fragment &theNucleus)
G4FragmentVectorBreakItUp (G4LayeredNucleus &nucleus)
void setVerboseLevel (const G4int verbose)

Detailed Description

Definition at line 41 of file G4BertiniEvaporation.hh.


Constructor & Destructor Documentation

G4BertiniEvaporation::G4BertiniEvaporation (  ) 

Definition at line 48 of file G4BertiniEvaporation.cc.

References G4cout, and G4endl.

00049 {
00050     G4cout << "Info G4BertiniEvaporation: This is still very fresh code."<< G4endl;
00051     G4cout << "     G4BertiniEvaporation: feed-back for improvement is very wellcome."<< G4endl;
00052     verboseLevel = 0;
00053  
00054     channelVector.push_back( new G4BENeutronChannel );
00055     channelVector.push_back( new G4BEProtonChannel  );
00056     channelVector.push_back( new G4BEDeuteronChannel);
00057     channelVector.push_back( new G4BETritonChannel  );
00058     channelVector.push_back( new G4BEHe3Channel     );
00059     channelVector.push_back( new G4BEHe4Channel     );
00060 }

G4BertiniEvaporation::~G4BertiniEvaporation (  ) 

Definition at line 63 of file G4BertiniEvaporation.cc.

00064 {
00065   while ( !channelVector.empty() )
00066     {
00067       delete channelVector.back();
00068       channelVector.pop_back();
00069     }
00070 }


Member Function Documentation

G4FragmentVector * G4BertiniEvaporation::BreakItUp ( G4LayeredNucleus nucleus  ) 

Definition at line 84 of file G4BertiniEvaporation.cc.

References G4BEGammaDeexcitation::emit(), G4BertiniEvaporationChannel::emit(), G4cout, G4endl, G4UniformRand, G4DynamicParticle::Get4Momentum(), G4Nucleus::GetA_asInt(), G4DynamicParticle::GetDefinition(), G4Nucleus::GetEnergyDeposit(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4DynamicParticle::GetMomentum(), G4LayeredNucleus::GetMomentum(), G4BertiniEvaporationChannel::getName(), G4NucleiProperties::GetNuclearMass(), G4BertiniEvaporationChannel::getParticleA(), G4BertiniEvaporationChannel::getParticleZ(), G4BertiniEvaporationChannel::getQ(), G4DynamicParticle::GetTotalEnergy(), G4DynamicParticle::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4BEGammaDeexcitation::setExcitationEnergy(), G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4BEGammaDeexcitation::setNucleusA(), G4BEGammaDeexcitation::setNucleusZ(), G4BertiniEvaporationChannel::setProbability(), and G4BEGammaDeexcitation::setVerboseLevel().

00085 {
00086   G4int nucleusA;
00087   G4int nucleusZ;
00088   G4int i;
00089   G4double totalProbability;
00090   G4double excE;
00091   G4double newExcitation;
00092   G4double nucleusTotalMomentum;
00093   G4double nucleusKineticEnergy;
00094   G4double mRes; // Mass of residual nucleus.
00095   G4ThreeVector nucleusMomentumVector;
00096   G4DynamicParticle *pEmittedParticle;
00097   std::vector< G4DynamicParticle * > secondaryParticleVector;
00098   G4FragmentVector * result = new G4FragmentVector;
00099   
00100   // Read properties of the nucleus.
00101   nucleusA = nucleus.GetA_asInt();
00102   nucleusZ = nucleus.GetZ_asInt();
00103   excE                  = nucleus.GetEnergyDeposit();
00104   nucleusMomentumVector = nucleus.GetMomentum();
00105 
00106   // Move to CMS frame, save initial velocity of the nucleus to boostToLab vector.
00107   G4ThreeVector boostToLab( ( 1/G4NucleiProperties::GetNuclearMass( nucleusA, nucleusZ ) ) 
00108                             * nucleusMomentumVector ); // xx mass ok?
00109 
00110   if ( verboseLevel >= 10 )
00111     G4cout << " G4BertiniEvaporation : initial kinematics : boostToLab vector = " << boostToLab << G4endl
00112            << "                     excitation energy  : " << excE<< G4endl;
00113 
00114   if ( nucleusA == 8 && nucleusZ == 4 )
00115     {
00116       splitBe8( excE, boostToLab, secondaryParticleVector );
00117       fillResult( secondaryParticleVector, result);
00118       return result;
00119     }
00120   
00121   // Initialize evaporation channels and calculate sum of emission
00122   // probabilities.
00123   std::vector< G4BertiniEvaporationChannel * >::iterator iChannel = channelVector.begin();  
00124   totalProbability = 0;
00125   for ( ; iChannel != channelVector.end() ; *iChannel++ )                  
00126      {
00127        ( *iChannel )->setNucleusA( nucleusA );
00128        ( *iChannel )->setNucleusZ( nucleusZ );
00129        ( *iChannel )->setExcitationEnergy( excE );
00130        ( *iChannel )->setPairingCorrection( 1 );
00131        ( *iChannel )->calculateProbability();
00132        totalProbability += ( *iChannel )->getProbability();
00133      }
00134 
00135   // If no evaporation process is possible, try without pairing energy.
00136   if ( totalProbability == 0 )
00137     {
00138       if ( verboseLevel >= 4 )
00139         G4cout << "G4BertiniEvaporation : no emission possible with pairing correction, trying without it" << G4endl;
00140       for ( iChannel = channelVector.begin() ; iChannel != channelVector.end() ; *iChannel++ )
00141         {
00142           ( *iChannel )->setPairingCorrection(0);
00143           ( *iChannel )->calculateProbability();
00144           totalProbability += ( *iChannel )->getProbability();
00145         }
00146       if ( verboseLevel >= 4 )
00147         G4cout << "                       probability without correction " << totalProbability << G4endl;
00148 
00149     }
00150 
00151   // Normal evaporation cycle follows.
00152   // Particles are evaporated till all probabilities are zero.
00153   while ( totalProbability > 0 )
00154     {
00155       G4BertiniEvaporationChannel *pSelectedChannel;
00156 
00157       // Sample active emission channel.
00158       G4double sampledProbability = G4UniformRand() * totalProbability;
00159       if ( verboseLevel >= 10 ) G4cout << "          RandomProb : " << sampledProbability << G4endl;
00160       for ( iChannel = channelVector.begin() ; iChannel != channelVector.end() ; *iChannel++ )
00161         {
00162           sampledProbability -= ( *iChannel )->getProbability();
00163           if ( sampledProbability < 0 ) break;
00164         }
00165       pSelectedChannel = ( *iChannel );
00166       if ( iChannel == channelVector.end() )
00167         throw G4HadronicException(__FILE__, __LINE__,  "G4BertiniEvaporation : Error while sampling evaporation particle" );
00168 
00169       if ( verboseLevel >= 4 ) 
00170         G4cout << "G4BertiniEvaporation : particle " << pSelectedChannel->getName() << " selected " << G4endl;
00171 
00172       // Max 10 tries to get a physically acceptable particle energy
00173       // in this emission channel.
00174       i = 0;
00175 
00176       do
00177         {
00178           pEmittedParticle = pSelectedChannel->emit();
00179           // This loop checks that particle with too large energy is not emitted.
00180           // CMS frame is considered in this loop. Nonrelativistic treatment. xxx
00181 
00182 
00183           const G4int zRes = nucleusZ - pSelectedChannel->getParticleZ(); 
00184           const G4int aRes  = nucleusA - pSelectedChannel->getParticleA(); 
00185           // const G4double eBind = G4NucleiProperties::GetBindingEnergy( aRes, zRes );  // Binding energy of the nucleus.
00186           mRes  = G4NucleiProperties::GetNuclearMass( aRes, zRes ); // Mass of the target nucleus
00187           //      In HETC88:
00188           //      eBind = Z * (-0.78244) + A * 8.36755 - cameron ( A , Z );
00189           //      mRes =  zRes * 938.79304 + ( aRes - zRes ) * 939.57548 - eBind; 
00190           //      where cameron ( A, Z ) is the mass excess by Cameron, see Canadian Journal of Physics,
00191           //      vol. 35, 1957, p.1022
00192 
00193           nucleusTotalMomentum = pEmittedParticle->GetTotalMomentum(); // CMS frame
00194           nucleusKineticEnergy = std::pow( nucleusTotalMomentum, 2 ) / ( 2 * mRes );
00195           newExcitation = excE - pEmittedParticle->GetKineticEnergy() - nucleusKineticEnergy - pSelectedChannel->getQ();
00196 
00197           if ( verboseLevel >= 10)
00198             G4cout << "G4BertiniEvaporation : Kinematics " << G4endl
00199                    << "                    part kinetic E in CMS = " 
00200                    << pEmittedParticle->GetKineticEnergy() << G4endl
00201                    << "                    new excitation E      =  " 
00202                    << newExcitation << G4endl;
00203 
00204           i++;
00205           if ( !( newExcitation > 0 ) && i < 10) delete pEmittedParticle;
00206         } while ( !( newExcitation > 0 )  &&  i < 10 );
00207 
00208       if ( i >= 10 ) 
00209         {
00210           // No appropriate particle energy found.
00211           // Set probability of this channel to zero 
00212           // and try to sample another particle on the 
00213           // next round.
00214           delete pEmittedParticle;
00215 
00216           if ( verboseLevel >= 4 )
00217             G4cout << "G4BertiniEvaporation : No appropriate energy for particle " 
00218                    << pSelectedChannel->getName() << " found." << G4endl;
00219 
00220           pSelectedChannel->setProbability( 0 );
00221 
00222           totalProbability = 0;
00223           for ( ; iChannel != channelVector.end() ; *iChannel++ )                  
00224             totalProbability += ( *iChannel )->getProbability();
00225         } // Treatment of physically unacceptable particle ends.
00226       else
00227         {
00228           // Good particle found. 
00229 
00230           // Transform particle properties to lab frame and save it.
00231           G4LorentzVector particleFourVector = pEmittedParticle->Get4Momentum();
00232           G4LorentzVector nucleusFourVector(  - pEmittedParticle->GetMomentum(), // CMS Frame.
00233                                               nucleusKineticEnergy + mRes ); // Total energy.
00234           particleFourVector.boost( boostToLab );
00235           nucleusFourVector.boost(  boostToLab );
00236           G4DynamicParticle *pPartLab  = new G4DynamicParticle( pEmittedParticle->GetDefinition(),
00237                                                                 particleFourVector.e(), // Total energy
00238                                                                 particleFourVector.vect() ); // momentum vector
00239           secondaryParticleVector.push_back( pPartLab );
00240 
00241           // Update residual nucleus and boostToLab vector.
00242           nucleusA = nucleusA - pSelectedChannel->getParticleA();
00243           nucleusZ = nucleusZ - pSelectedChannel->getParticleZ();
00244           excE = newExcitation;
00245           boostToLab = G4ThreeVector ( ( 1/mRes ) * nucleusFourVector.vect() );
00246 
00247           if ( verboseLevel >= 10 )
00248             G4cout << "  particle mom in cms " << pEmittedParticle->GetMomentum() << G4endl 
00249                    << "  particle mom in cms " << pEmittedParticle->GetTotalMomentum() << G4endl 
00250                    << "          Etot in cms " << pEmittedParticle->GetTotalEnergy() << G4endl
00251                    << "          Ekin in cms " << pEmittedParticle->GetKineticEnergy() << G4endl
00252                    << "        particle mass " << pEmittedParticle->GetMass() << G4endl
00253                    << "          boost  vect " << boostToLab << G4endl
00254                    << "          boosted 4v  " << particleFourVector << G4endl
00255                    << "          nucleus 4v  " << nucleusFourVector << G4endl
00256                    << "          nucl cm mom " << nucleusTotalMomentum << G4endl
00257                    << " particle k.e. in lab " << pPartLab->GetKineticEnergy() << G4endl
00258                    << "     new boost vector " << boostToLab << G4endl;
00259 
00260           delete pEmittedParticle;
00261 
00262           // If the residual nucleus is Be8, split it.
00263           if ( nucleusA == 8 && nucleusZ == 4 )
00264             {
00265               splitBe8( excE, boostToLab, secondaryParticleVector );
00266               fillResult( secondaryParticleVector, result);
00267               return result;
00268             }
00269 
00270           // Update evaporation channels and 
00271           // total emission probability.
00272           totalProbability = 0;
00273           for ( iChannel = channelVector.begin() ; iChannel != channelVector.end() ; *iChannel++ )  
00274             {
00275               ( *iChannel )->setNucleusA( nucleusA ); 
00276               ( *iChannel )->setNucleusZ( nucleusZ );
00277               ( *iChannel )->setExcitationEnergy( excE ); 
00278               ( *iChannel )->calculateProbability();
00279               totalProbability += ( *iChannel )->getProbability();
00280             }
00281         } // Treatment of physically acceptable particle ends.
00282     } // End of evaporation cycle.
00283 
00284   // Gamma de-excitation.
00285   G4BEGammaDeexcitation *pGammaChannel = new G4BEGammaDeexcitation;
00286   pGammaChannel->setNucleusA( nucleusA );
00287   pGammaChannel->setNucleusZ( nucleusZ );
00288   pGammaChannel->setVerboseLevel( verboseLevel );
00289   pGammaChannel->setExcitationEnergy( excE );
00290 
00291   // Emit equally sampled photons until all the excitation energy is
00292   // used; the last photon gets the remaining de-excitation energy.
00293   // Change of momentum of the nucleus is neglected.
00294   G4double gammaEnergy;
00295   G4double totalDeExcEnergy = 0;
00296 
00297   while ( excE > 0 )
00298     {
00299       pEmittedParticle = pGammaChannel->emit();
00300       gammaEnergy = pEmittedParticle->GetKineticEnergy();
00301       
00302       if ( ( totalDeExcEnergy + gammaEnergy ) > excE ) 
00303         {
00304           // All the remaining energy is used here.
00305           gammaEnergy = excE - totalDeExcEnergy;
00306           excE = 0;
00307         }
00308 
00309       totalDeExcEnergy += gammaEnergy; 
00310 
00311       if ( verboseLevel >= 10 )
00312         G4cout << " G4BertiniEvaporation : gamma de-excitation, g of " << gammaEnergy << " MeV " << G4endl;
00313       
00314       G4LorentzVector gammaFourVector( pEmittedParticle->GetMomentum(),
00315                                        pEmittedParticle->GetKineticEnergy() );
00316       gammaFourVector.boost( boostToLab );
00317       pEmittedParticle->SetKineticEnergy( gammaFourVector.e() );
00318       pEmittedParticle->SetMomentumDirection( gammaFourVector.vect().unit() );
00319 
00320       secondaryParticleVector.push_back( pEmittedParticle );
00321 
00322       } 
00323 
00324   delete pGammaChannel;
00325 
00326   // Residual nucleus is not returned.
00327 
00328   fillResult( secondaryParticleVector, result);
00329   
00330   return result;
00331 }

virtual G4FragmentVector* G4BertiniEvaporation::BreakItUp ( const G4Fragment theNucleus  )  [inline, virtual]

Implements G4VEvaporation.

Definition at line 48 of file G4BertiniEvaporation.hh.

References G4Fragment::GetA(), G4Fragment::GetExcitationEnergy(), and G4Fragment::GetZ().

00049   {
00050     G4LayeredNucleus aNuc( theNucleus.GetA(), theNucleus.GetZ() );
00051     aNuc.AddExcitationEnergy(theNucleus.GetExcitationEnergy());
00052     return BreakItUp(aNuc);
00053   }

void G4BertiniEvaporation::setVerboseLevel ( const G4int  verbose  ) 

Definition at line 73 of file G4BertiniEvaporation.cc.

00074 {
00075   verboseLevel = verbose;
00076 
00077   // Update verbose level to all evaporation channels.
00078   std::vector< G4BertiniEvaporationChannel * >::iterator iChannel = channelVector.begin();  
00079   for ( ; iChannel != channelVector.end() ; *iChannel++ )                  
00080     ( *iChannel )->setVerboseLevel( verboseLevel );
00081 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:51:31 2013 for Geant4 by  doxygen 1.4.7