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
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00078
00079 #include <iomanip>
00080 #include <numeric>
00081
00082 #include "G4WilsonAblationModel.hh"
00083 #include "G4PhysicalConstants.hh"
00084 #include "G4SystemOfUnits.hh"
00085 #include "Randomize.hh"
00086 #include "G4ParticleTable.hh"
00087 #include "G4IonTable.hh"
00088 #include "G4Alpha.hh"
00089 #include "G4He3.hh"
00090 #include "G4Triton.hh"
00091 #include "G4Deuteron.hh"
00092 #include "G4Proton.hh"
00093 #include "G4Neutron.hh"
00094 #include "G4AlphaEvaporationChannel.hh"
00095 #include "G4He3EvaporationChannel.hh"
00096 #include "G4TritonEvaporationChannel.hh"
00097 #include "G4DeuteronEvaporationChannel.hh"
00098 #include "G4ProtonEvaporationChannel.hh"
00099 #include "G4NeutronEvaporationChannel.hh"
00100 #include "G4LorentzVector.hh"
00101 #include "G4VEvaporationChannel.hh"
00102
00104
00105 G4WilsonAblationModel::G4WilsonAblationModel()
00106 {
00107
00108
00109
00110
00111 PrintWelcomeMessage();
00112
00113
00114
00115
00116 verboseLevel = 0;
00117
00118
00119
00120
00121
00122 B = 10.0 * MeV;
00123
00124
00125
00126
00127
00128 produceSecondaries = true;
00129
00130
00131
00132
00133
00134 nFragTypes = 6;
00135 fragType[0] = G4Alpha::Alpha();
00136 fragType[1] = G4He3::He3();
00137 fragType[2] = G4Triton::Triton();
00138 fragType[3] = G4Deuteron::Deuteron();
00139 fragType[4] = G4Proton::Proton();
00140 fragType[5] = G4Neutron::Neutron();
00141
00142
00143
00144
00145 verboseLevel = 0;
00146 theChannelFactory = new G4EvaporationFactory();
00147 theChannels = theChannelFactory->GetChannel();
00148
00149
00150
00151
00152
00153 OPTxs = 3;
00154 useSICB = false;
00155 fragmentVector = 0;
00156 }
00158
00159 G4WilsonAblationModel::~G4WilsonAblationModel()
00160 {
00161 if (theChannels != 0) theChannels = 0;
00162 if (theChannelFactory != 0) delete theChannelFactory;
00163 }
00165
00166 G4FragmentVector *G4WilsonAblationModel::BreakItUp
00167 (const G4Fragment &theNucleus)
00168 {
00169
00170
00171
00172
00173
00174 fragmentVector = new G4FragmentVector;
00175 fragmentVector->clear();
00176
00177
00178
00179
00180 G4int A = theNucleus.GetA_asInt();
00181 G4int Z = theNucleus.GetZ_asInt();
00182 G4double ex = theNucleus.GetExcitationEnergy();
00183 if (verboseLevel >= 2)
00184 {
00185 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
00186 <<"oooooooooooooooooooooooooooooooooooooooo"
00187 <<G4endl;
00188 G4cout.precision(6);
00189 G4cout <<"IN G4WilsonAblationModel" <<G4endl;
00190 G4cout <<"Initial prefragment A=" <<A
00191 <<", Z=" <<Z
00192 <<", excitation energy = " <<ex/MeV <<" MeV"
00193 <<G4endl;
00194 }
00195
00196
00197
00198
00199
00200
00201 if (A == 0)
00202 {
00203 if (verboseLevel >= 2)
00204 {
00205 G4cout <<"No nucleus to decay" <<G4endl;
00206 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
00207 <<"oooooooooooooooooooooooooooooooooooooooo"
00208 <<G4endl;
00209 }
00210 return fragmentVector;
00211 }
00212 else if (A == 1)
00213 {
00214 G4LorentzVector lorentzVector = theNucleus.GetMomentum();
00215 lorentzVector.setE(lorentzVector.e()-ex+10.0*eV);
00216 if (Z == 0)
00217 {
00218 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Neutron::Neutron());
00219 fragmentVector->push_back(fragment);
00220 }
00221 else
00222 {
00223 G4Fragment *fragment = new G4Fragment(lorentzVector,G4Proton::Proton());
00224 fragmentVector->push_back(fragment);
00225 }
00226 if (verboseLevel >= 2)
00227 {
00228 G4cout <<"Final fragment is in fact only a nucleon) :" <<G4endl;
00229 G4cout <<(*fragmentVector)[0] <<G4endl;
00230 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
00231 <<"oooooooooooooooooooooooooooooooooooooooo"
00232 <<G4endl;
00233 }
00234 return fragmentVector;
00235 }
00236
00237
00238
00239
00240
00241 G4int DAabl = (G4int) (ex / B);
00242 if (DAabl > A) DAabl = A;
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252 G4int AF = A - DAabl;
00253 G4int ZF = 0;
00254 if (AF > 0)
00255 {
00256 G4double AFd = (G4double) AF;
00257 G4double R = 11.8 / std::pow(AFd, 0.45);
00258 G4int minZ = Z - DAabl;
00259 if (minZ <= 0) minZ = 1;
00260
00261
00262
00263
00264
00265 G4double sig[100];
00266 G4double sum = 0.0;
00267 for (G4int ii=minZ; ii<= Z; ii++)
00268 {
00269 sum += std::exp(-R*std::pow(std::abs(ii - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
00270 sig[ii] = sum;
00271 }
00272
00273
00274
00275
00276 G4double xi = G4UniformRand();
00277 G4int iz = minZ;
00278 G4bool found = false;
00279 while (iz <= Z && !found)
00280 {
00281 found = (xi <= sig[iz]/sum);
00282 if (!found) iz++;
00283 }
00284 if (iz > Z)
00285 ZF = Z;
00286 else
00287 ZF = iz;
00288 }
00289 G4int DZabl = Z - ZF;
00290
00291
00292
00293
00294
00295
00296
00297
00298 G4double totalEpost = 0.0;
00299 evapType.clear();
00300 for (G4int ift=0; ift<nFragTypes; ift++)
00301 {
00302 G4ParticleDefinition *type = fragType[ift];
00303 G4double n = std::floor((G4double) DAabl / type->GetBaryonNumber() + 1.0E-10);
00304 G4double n1 = 1.0E+10;
00305 if (fragType[ift]->GetPDGCharge() > 0.0)
00306 n1 = std::floor((G4double) DZabl / type->GetPDGCharge() + 1.0E-10);
00307 if (n > n1) n = n1;
00308 if (n > 0.0)
00309 {
00310 G4double mass = type->GetPDGMass();
00311 for (G4int j=0; j<(G4int) n; j++)
00312 {
00313 totalEpost += mass;
00314 evapType.push_back(type);
00315 }
00316 DAabl -= (G4int) (n * type->GetBaryonNumber() + 1.0E-10);
00317 DZabl -= (G4int) (n * type->GetPDGCharge() + 1.0E-10);
00318 }
00319 }
00320
00321
00322
00323
00324
00325
00326
00327
00328 G4double massFinalFrag = 0.0;
00329 if (AF > 0)
00330 massFinalFrag = G4ParticleTable::GetParticleTable()->GetIonTable()->
00331 GetIonMass(ZF,AF);
00332 else
00333 {
00334 G4ParticleDefinition *type = evapType[evapType.size()-1];
00335 AF = type->GetBaryonNumber();
00336 ZF = (G4int) (type->GetPDGCharge() + 1.0E-10);
00337 evapType.erase(evapType.end()-1);
00338 }
00339 totalEpost += massFinalFrag;
00340
00341
00342
00343
00344 if (verboseLevel >= 2)
00345 {
00346 G4cout <<"Final fragment A=" <<AF
00347 <<", Z=" <<ZF
00348 <<G4endl;
00349 for (G4int ift=0; ift<nFragTypes; ift++)
00350 {
00351 G4ParticleDefinition *type = fragType[ift];
00352 G4int n = std::count(evapType.begin(),evapType.end(),type);
00353 if (n > 0)
00354 G4cout <<"Particle type: " <<std::setw(10) <<type->GetParticleName()
00355 <<", number of particles emitted = " <<n <<G4endl;
00356 }
00357 }
00358
00359
00360
00361
00362
00363 G4double massPreFrag = theNucleus.GetGroundStateMass();
00364 G4double totalEpre = massPreFrag + ex;
00365 G4double excess = totalEpre - totalEpost;
00366
00367 G4Fragment *resultNucleus = new G4Fragment(A, Z, theNucleus.GetMomentum());
00368 G4ThreeVector boost(0.0,0.0,0.0);
00369 G4int nEvap = 0;
00370 if (produceSecondaries && evapType.size()>0)
00371 {
00372 if (excess > 0.0)
00373 {
00374 SelectSecondariesByEvaporation (resultNucleus);
00375 nEvap = fragmentVector->size();
00376 boost = resultNucleus->GetMomentum().findBoostToCM();
00377 if (evapType.size() > 0)
00378 SelectSecondariesByDefault (boost);
00379 }
00380 else
00381 SelectSecondariesByDefault(G4ThreeVector(0.0,0.0,0.0));
00382 }
00383
00384 if (AF > 0)
00385 {
00386 G4double mass = G4ParticleTable::GetParticleTable()->GetIonTable()->
00387 GetIonMass(ZF,AF);
00388 G4double e = mass + 10.0*eV;
00389 G4double p = std::sqrt(e*e-mass*mass);
00390 G4ThreeVector direction(0.0,0.0,1.0);
00391 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
00392 lorentzVector.boost(-boost);
00393 G4Fragment* frag = new G4Fragment(AF, ZF, lorentzVector);
00394 fragmentVector->push_back(frag);
00395 }
00396 delete resultNucleus;
00397
00398
00399
00400
00401 if (verboseLevel >= 2)
00402 {
00403 if (nEvap > 0)
00404 {
00405 G4cout <<"----------------------" <<G4endl;
00406 G4cout <<"Evaporated particles :" <<G4endl;
00407 G4cout <<"----------------------" <<G4endl;
00408 }
00409 G4int ie = 0;
00410 G4FragmentVector::iterator iter;
00411 for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++)
00412 {
00413 if (ie == nEvap)
00414 {
00415
00416 G4cout <<"---------------------------------" <<G4endl;
00417 G4cout <<"Particles from default emission :" <<G4endl;
00418 G4cout <<"---------------------------------" <<G4endl;
00419 }
00420 G4cout <<*iter <<G4endl;
00421 }
00422 G4cout <<"oooooooooooooooooooooooooooooooooooooooo"
00423 <<"oooooooooooooooooooooooooooooooooooooooo"
00424 <<G4endl;
00425 }
00426
00427 return fragmentVector;
00428 }
00430
00431 void G4WilsonAblationModel::SelectSecondariesByEvaporation
00432 (G4Fragment *intermediateNucleus)
00433 {
00434 G4Fragment theResidualNucleus = *intermediateNucleus;
00435 G4bool evaporate = true;
00436 while (evaporate && evapType.size() != 0)
00437 {
00438
00439
00440
00441
00442
00443
00444 std::vector <G4VEvaporationChannel*> theChannels1;
00445 theChannels1.clear();
00446 std::vector <G4VEvaporationChannel*>::iterator i;
00447 VectorOfFragmentTypes::iterator iter;
00448 std::vector <VectorOfFragmentTypes::iterator> iters;
00449 iters.clear();
00450 iter = std::find(evapType.begin(), evapType.end(), G4Alpha::Alpha());
00451 if (iter != evapType.end())
00452 {
00453 theChannels1.push_back(new G4AlphaEvaporationChannel);
00454 i = theChannels1.end() - 1;
00455 (*i)->SetOPTxs(OPTxs);
00456 (*i)->UseSICB(useSICB);
00457
00458 iters.push_back(iter);
00459 }
00460 iter = std::find(evapType.begin(), evapType.end(), G4He3::He3());
00461 if (iter != evapType.end())
00462 {
00463 theChannels1.push_back(new G4He3EvaporationChannel);
00464 i = theChannels1.end() - 1;
00465 (*i)->SetOPTxs(OPTxs);
00466 (*i)->UseSICB(useSICB);
00467
00468 iters.push_back(iter);
00469 }
00470 iter = std::find(evapType.begin(), evapType.end(), G4Triton::Triton());
00471 if (iter != evapType.end())
00472 {
00473 theChannels1.push_back(new G4TritonEvaporationChannel);
00474 i = theChannels1.end() - 1;
00475 (*i)->SetOPTxs(OPTxs);
00476 (*i)->UseSICB(useSICB);
00477
00478 iters.push_back(iter);
00479 }
00480 iter = std::find(evapType.begin(), evapType.end(), G4Deuteron::Deuteron());
00481 if (iter != evapType.end())
00482 {
00483 theChannels1.push_back(new G4DeuteronEvaporationChannel);
00484 i = theChannels1.end() - 1;
00485 (*i)->SetOPTxs(OPTxs);
00486 (*i)->UseSICB(useSICB);
00487
00488 iters.push_back(iter);
00489 }
00490 iter = std::find(evapType.begin(), evapType.end(), G4Proton::Proton());
00491 if (iter != evapType.end())
00492 {
00493 theChannels1.push_back(new G4ProtonEvaporationChannel);
00494 i = theChannels1.end() - 1;
00495 (*i)->SetOPTxs(OPTxs);
00496 (*i)->UseSICB(useSICB);
00497
00498 iters.push_back(iter);
00499 }
00500 iter = std::find(evapType.begin(), evapType.end(), G4Neutron::Neutron());
00501 if (iter != evapType.end())
00502 {
00503 theChannels1.push_back(new G4NeutronEvaporationChannel);
00504 i = theChannels1.end() - 1;
00505 (*i)->SetOPTxs(OPTxs);
00506 (*i)->UseSICB(useSICB);
00507
00508 iters.push_back(iter);
00509 }
00510 G4int nChannels = theChannels1.size();
00511
00512 G4double totalProb = 0.0;
00513 G4int ich = 0;
00514 G4double probEvapType[6] = {0.0};
00515 std::vector<G4VEvaporationChannel*>::iterator iterEv;
00516 for (iterEv=theChannels1.begin(); iterEv!=theChannels1.end(); iterEv++) {
00517 totalProb += (*iterEv)->GetEmissionProbability(intermediateNucleus);
00518 probEvapType[ich] = totalProb;
00519 ++ich;
00520 }
00521 if (totalProb > 0.0) {
00522
00523
00524
00525
00526
00527
00528 G4double xi = totalProb*G4UniformRand();
00529 G4int ii = 0;
00530 for (ii=0; ii<nChannels; ii++) {
00531 if (xi < probEvapType[ii]) { break; }
00532 }
00533 if (ii >= nChannels) { ii = nChannels - 1; }
00534 G4FragmentVector *evaporationResult = theChannels1[ii]->
00535 BreakUp(*intermediateNucleus);
00536 fragmentVector->push_back((*evaporationResult)[0]);
00537 *intermediateNucleus = *(*evaporationResult)[1];
00538
00539 delete evaporationResult;
00540
00541 }
00542 else
00543 {
00544
00545
00546
00547
00548
00549 evaporate = false;
00550 }
00551 }
00552
00553 return;
00554 }
00556
00557 void G4WilsonAblationModel::SelectSecondariesByDefault (G4ThreeVector boost)
00558 {
00559 for (unsigned i=0; i<evapType.size(); i++)
00560 {
00561 G4ParticleDefinition *type = evapType[i];
00562 G4double mass = type->GetPDGMass();
00563 G4double e = mass + 10.0*eV;
00564 G4double p = std::sqrt(e*e-mass*mass);
00565 G4double costheta = 2.0*G4UniformRand() - 1.0;
00566 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
00567 G4double phi = twopi * G4UniformRand() * rad;
00568 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
00569 G4LorentzVector lorentzVector = G4LorentzVector(direction*p, e);
00570 lorentzVector.boost(-boost);
00571
00572
00573
00574
00575 G4int A = type->GetBaryonNumber();
00576 G4int Z = (G4int) (type->GetPDGCharge() + 1.0E-10);
00577 G4Fragment *fragment =
00578 new G4Fragment(A, Z, lorentzVector);
00579
00580 fragmentVector->push_back(fragment);
00581 }
00582 }
00584
00585 void G4WilsonAblationModel::PrintWelcomeMessage ()
00586 {
00587 G4cout <<G4endl;
00588 G4cout <<" *****************************************************************"
00589 <<G4endl;
00590 G4cout <<" Nuclear ablation model for nuclear-nuclear interactions activated"
00591 <<G4endl;
00592 G4cout <<" (Written by QinetiQ Ltd for the European Space Agency)"
00593 <<G4endl;
00594 G4cout <<" *****************************************************************"
00595 <<G4endl;
00596 G4cout << G4endl;
00597
00598 return;
00599 }
00601