00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "globals.hh"
00032 #include "G4BertiniEvaporation.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4BENeutronChannel.hh"
00036 #include "G4BEProtonChannel.hh"
00037 #include "G4BEDeuteronChannel.hh"
00038 #include "G4BETritonChannel.hh"
00039 #include "G4BEHe3Channel.hh"
00040 #include "G4BEHe4Channel.hh"
00041 #include "G4BEGammaDeexcitation.hh"
00042
00043
00044
00045
00046
00047
00048 G4BertiniEvaporation::G4BertiniEvaporation()
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 }
00061
00062
00063 G4BertiniEvaporation::~G4BertiniEvaporation()
00064 {
00065 while ( !channelVector.empty() )
00066 {
00067 delete channelVector.back();
00068 channelVector.pop_back();
00069 }
00070 }
00071
00072
00073 void G4BertiniEvaporation::setVerboseLevel( const G4int verbose )
00074 {
00075 verboseLevel = verbose;
00076
00077
00078 std::vector< G4BertiniEvaporationChannel * >::iterator iChannel = channelVector.begin();
00079 for ( ; iChannel != channelVector.end() ; *iChannel++ )
00080 ( *iChannel )->setVerboseLevel( verboseLevel );
00081 }
00082
00083
00084 G4FragmentVector * G4BertiniEvaporation::BreakItUp( G4LayeredNucleus & nucleus )
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;
00095 G4ThreeVector nucleusMomentumVector;
00096 G4DynamicParticle *pEmittedParticle;
00097 std::vector< G4DynamicParticle * > secondaryParticleVector;
00098 G4FragmentVector * result = new G4FragmentVector;
00099
00100
00101 nucleusA = nucleus.GetA_asInt();
00102 nucleusZ = nucleus.GetZ_asInt();
00103 excE = nucleus.GetEnergyDeposit();
00104 nucleusMomentumVector = nucleus.GetMomentum();
00105
00106
00107 G4ThreeVector boostToLab( ( 1/G4NucleiProperties::GetNuclearMass( nucleusA, nucleusZ ) )
00108 * nucleusMomentumVector );
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
00122
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
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
00152
00153 while ( totalProbability > 0 )
00154 {
00155 G4BertiniEvaporationChannel *pSelectedChannel;
00156
00157
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
00173
00174 i = 0;
00175
00176 do
00177 {
00178 pEmittedParticle = pSelectedChannel->emit();
00179
00180
00181
00182
00183 const G4int zRes = nucleusZ - pSelectedChannel->getParticleZ();
00184 const G4int aRes = nucleusA - pSelectedChannel->getParticleA();
00185
00186 mRes = G4NucleiProperties::GetNuclearMass( aRes, zRes );
00187
00188
00189
00190
00191
00192
00193 nucleusTotalMomentum = pEmittedParticle->GetTotalMomentum();
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
00211
00212
00213
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 }
00226 else
00227 {
00228
00229
00230
00231 G4LorentzVector particleFourVector = pEmittedParticle->Get4Momentum();
00232 G4LorentzVector nucleusFourVector( - pEmittedParticle->GetMomentum(),
00233 nucleusKineticEnergy + mRes );
00234 particleFourVector.boost( boostToLab );
00235 nucleusFourVector.boost( boostToLab );
00236 G4DynamicParticle *pPartLab = new G4DynamicParticle( pEmittedParticle->GetDefinition(),
00237 particleFourVector.e(),
00238 particleFourVector.vect() );
00239 secondaryParticleVector.push_back( pPartLab );
00240
00241
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
00263 if ( nucleusA == 8 && nucleusZ == 4 )
00264 {
00265 splitBe8( excE, boostToLab, secondaryParticleVector );
00266 fillResult( secondaryParticleVector, result);
00267 return result;
00268 }
00269
00270
00271
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 }
00282 }
00283
00284
00285 G4BEGammaDeexcitation *pGammaChannel = new G4BEGammaDeexcitation;
00286 pGammaChannel->setNucleusA( nucleusA );
00287 pGammaChannel->setNucleusZ( nucleusZ );
00288 pGammaChannel->setVerboseLevel( verboseLevel );
00289 pGammaChannel->setExcitationEnergy( excE );
00290
00291
00292
00293
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
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
00327
00328 fillResult( secondaryParticleVector, result);
00329
00330 return result;
00331 }
00332
00333
00334 void G4BertiniEvaporation::splitBe8( const G4double E,
00335 const G4ThreeVector boostToLab,
00336 std::vector< G4DynamicParticle * > & secondaryParticleVector )
00337 {
00338 G4double kineticEnergy;
00339 G4double u;
00340 G4double v;
00341 G4double w;
00342 const G4double Be8DecayEnergy = 0.093 * MeV;
00343
00344 if ( E <= 0 ) throw G4HadronicException(__FILE__, __LINE__, "G4BertiniEvaporation : excitation energy < 0 " );
00345 if ( verboseLevel >= 4 ) G4cout << " Be8 split to 2 x He4" << G4endl;
00346
00347 kineticEnergy = 0.5 * ( E + Be8DecayEnergy );
00348
00349
00350 isotropicCosines( u , v , w );
00351 G4ThreeVector momentumDirectionCMS( u, v, w );
00352
00353 G4DynamicParticle *pP1Cms = new G4DynamicParticle( G4Alpha::Alpha(),
00354 momentumDirectionCMS,
00355 kineticEnergy );
00356 G4DynamicParticle *pP2Cms = new G4DynamicParticle( G4Alpha::Alpha(),
00357 -momentumDirectionCMS,
00358 kineticEnergy );
00359 G4LorentzVector fourVector1( pP1Cms->GetMomentum(),
00360 pP1Cms->GetTotalEnergy() );
00361 G4LorentzVector fourVector2( pP2Cms->GetMomentum(),
00362 pP2Cms->GetTotalEnergy() );
00363
00364
00365
00366 fourVector1.boost( boostToLab );
00367 fourVector2.boost( boostToLab );
00368 G4DynamicParticle * pP1Lab = new G4DynamicParticle( G4Alpha::Alpha(),
00369 fourVector1.vect(),
00370 fourVector1.e() );
00371 G4DynamicParticle * pP2Lab = new G4DynamicParticle( G4Alpha::Alpha(),
00372 fourVector2.vect(),
00373 fourVector2.e() );
00374 secondaryParticleVector.push_back( pP1Lab );
00375 secondaryParticleVector.push_back( pP2Lab );
00376
00377 delete pP1Cms;
00378 delete pP2Cms;
00379
00380 return;
00381 }
00382
00383
00384 void G4BertiniEvaporation::fillResult( std::vector<G4DynamicParticle *> secondaryParticleVector,
00385 G4FragmentVector * aResult )
00386 {
00387
00388 for ( size_t i = 0 ; i < secondaryParticleVector.size() ; i++ )
00389 {
00390 G4int aZ = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetPDGCharge() );
00391 G4int aA = static_cast<G4int> (secondaryParticleVector[i]->GetDefinition()->GetBaryonNumber());
00392 G4LorentzVector aMomentum = secondaryParticleVector[i]->Get4Momentum();
00393 if(aA>0)
00394 {
00395 aResult->push_back( new G4Fragment(aA, aZ, aMomentum) );
00396 }
00397 else
00398 {
00399 aResult->push_back( new G4Fragment(aMomentum, secondaryParticleVector[i]->GetDefinition()) );
00400 }
00401 }
00402 return;
00403 }
00404
00405
00406 void G4BertiniEvaporation::isotropicCosines( G4double & u, G4double & v, G4double & w )
00407 {
00408
00409 G4double CosTheta = 1.0 - 2.0 * G4UniformRand();
00410 G4double SinTheta = std::sqrt( 1.0 - CosTheta * CosTheta );
00411 G4double Phi = twopi * G4UniformRand();
00412
00413 u = std::cos( Phi ) * SinTheta;
00414 v = std::cos( Phi ) * CosTheta,
00415 w = std::sin( Phi );
00416
00417 return;
00418 }