#include <G4BertiniEvaporation.hh>
Inheritance diagram for G4BertiniEvaporation:
Public Member Functions | |
G4BertiniEvaporation () | |
~G4BertiniEvaporation () | |
virtual G4FragmentVector * | BreakItUp (const G4Fragment &theNucleus) |
G4FragmentVector * | BreakItUp (G4LayeredNucleus &nucleus) |
void | setVerboseLevel (const G4int verbose) |
Definition at line 41 of file G4BertiniEvaporation.hh.
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 }
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 }