175 fragmentVector->clear();
183 if (verboseLevel >= 2)
185 G4cout <<
"oooooooooooooooooooooooooooooooooooooooo"
186 <<
"oooooooooooooooooooooooooooooooooooooooo"
190 G4cout <<
"Initial prefragment A=" <<A
192 <<
", excitation energy = " <<ex/
MeV <<
" MeV"
203 if (verboseLevel >= 2)
206 G4cout <<
"oooooooooooooooooooooooooooooooooooooooo"
207 <<
"oooooooooooooooooooooooooooooooooooooooo"
210 return fragmentVector;
215 lorentzVector.setE(lorentzVector.e()-ex+10.0*
eV);
219 fragmentVector->push_back(fragment);
224 fragmentVector->push_back(fragment);
226 if (verboseLevel >= 2)
228 G4cout <<
"Final fragment is in fact only a nucleon) :" <<
G4endl;
230 G4cout <<
"oooooooooooooooooooooooooooooooooooooooo"
231 <<
"oooooooooooooooooooooooooooooooooooooooo"
234 return fragmentVector;
242 if (DAabl > A) DAabl =
A;
252 G4int AF = A - DAabl;
257 G4double R = 11.8 / std::pow(AFd, 0.45);
258 G4int minZ = Z - DAabl;
259 if (minZ <= 0) minZ = 1;
267 for (
G4int ii=minZ; ii<= Z; ii++)
269 sum += std::exp(-R*std::pow(std::abs(ii - 0.486*AFd + 3.8E-04*AFd*AFd),1.5));
279 while (iz <= Z && !found)
281 found = (xi <= sig[
iz]/sum);
289 G4int DZabl = Z - ZF;
300 for (
G4int ift=0; ift<nFragTypes; ift++)
305 if (fragType[ift]->GetPDGCharge() > 0.0)
314 evapType.push_back(type);
337 evapType.erase(evapType.end()-1);
339 totalEpost += massFinalFrag;
344 if (verboseLevel >= 2)
346 G4cout <<
"Final fragment A=" <<AF
349 for (
G4int ift=0; ift<nFragTypes; ift++)
352 G4int n = std::count(evapType.begin(),evapType.end(),type);
355 <<
", number of particles emitted = " <<n <<
G4endl;
364 G4double totalEpre = massPreFrag + ex;
365 G4double excess = totalEpre - totalEpost;
370 if (produceSecondaries && evapType.size()>0)
374 SelectSecondariesByEvaporation (resultNucleus);
375 nEvap = fragmentVector->size();
376 boost = resultNucleus->
GetMomentum().findBoostToCM();
377 if (evapType.size() > 0)
378 SelectSecondariesByDefault (boost);
389 G4double p = std::sqrt(e*e-mass*mass);
392 lorentzVector.boost(-boost);
394 fragmentVector->push_back(frag);
396 delete resultNucleus;
401 if (verboseLevel >= 2)
410 G4FragmentVector::iterator iter;
411 for (iter = fragmentVector->begin(); iter != fragmentVector->end(); iter++)
422 G4cout <<
"oooooooooooooooooooooooooooooooooooooooo"
423 <<
"oooooooooooooooooooooooooooooooooooooooo"
427 return fragmentVector;
434 G4Fragment theResidualNucleus = *intermediateNucleus;
436 while (evaporate && evapType.size() != 0)
444 std::vector <G4VEvaporationChannel*> theChannels1;
445 theChannels1.clear();
446 std::vector <G4VEvaporationChannel*>::iterator i;
447 VectorOfFragmentTypes::iterator iter;
448 std::vector <VectorOfFragmentTypes::iterator> iters;
450 iter = std::find(evapType.begin(), evapType.end(),
G4Alpha::Alpha());
451 if (iter != evapType.end())
454 i = theChannels1.end() - 1;
455 (*i)->SetOPTxs(OPTxs);
456 (*i)->UseSICB(useSICB);
458 iters.push_back(iter);
460 iter = std::find(evapType.begin(), evapType.end(),
G4He3::He3());
461 if (iter != evapType.end())
464 i = theChannels1.end() - 1;
465 (*i)->SetOPTxs(OPTxs);
466 (*i)->UseSICB(useSICB);
468 iters.push_back(iter);
471 if (iter != evapType.end())
474 i = theChannels1.end() - 1;
475 (*i)->SetOPTxs(OPTxs);
476 (*i)->UseSICB(useSICB);
478 iters.push_back(iter);
481 if (iter != evapType.end())
484 i = theChannels1.end() - 1;
485 (*i)->SetOPTxs(OPTxs);
486 (*i)->UseSICB(useSICB);
488 iters.push_back(iter);
491 if (iter != evapType.end())
494 i = theChannels1.end() - 1;
495 (*i)->SetOPTxs(OPTxs);
496 (*i)->UseSICB(useSICB);
498 iters.push_back(iter);
501 if (iter != evapType.end())
504 i = theChannels1.end() - 1;
505 (*i)->SetOPTxs(OPTxs);
506 (*i)->UseSICB(useSICB);
508 iters.push_back(iter);
510 G4int nChannels = theChannels1.size();
515 std::vector<G4VEvaporationChannel*>::iterator iterEv;
516 for (iterEv=theChannels1.begin(); iterEv!=theChannels1.end(); iterEv++) {
517 totalProb += (*iterEv)->GetEmissionProbability(intermediateNucleus);
518 probEvapType[ich] = totalProb;
521 if (totalProb > 0.0) {
530 for (ii=0; ii<nChannels; ii++) {
531 if (xi < probEvapType[ii]) {
break; }
533 if (ii >= nChannels) { ii = nChannels - 1; }
535 BreakUp(*intermediateNucleus);
536 fragmentVector->push_back((*evaporationResult)[0]);
537 *intermediateNucleus = *(*evaporationResult)[1];
539 delete evaporationResult;
559 for (
unsigned i=0; i<
evapType.size(); i++)
564 G4double p = std::sqrt(e*e-mass*mass);
566 G4double sintheta = std::sqrt((1.0 - costheta)*(1.0 + costheta));
568 G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),costheta);
570 lorentzVector.boost(-boost);
588 G4cout <<
" *****************************************************************"
590 G4cout <<
" Nuclear ablation model for nuclear-nuclear interactions activated"
592 G4cout <<
" (Written by QinetiQ Ltd for the European Space Agency)"
594 G4cout <<
" !!! WARNING: This model is not well validation and should not be used for accurate simulation !!!"
596 G4cout <<
" *****************************************************************"
virtual std::vector< G4VEvaporationChannel * > * GetChannel()=0
void SelectSecondariesByDefault(G4ThreeVector)
VectorOfFragmentTypes evapType
CLHEP::Hep3Vector G4ThreeVector
G4VEvaporationFactory * theChannelFactory
G4ParticleDefinition * fragType[6]
void SelectSecondariesByEvaporation(G4Fragment *)
const G4String & GetParticleName() const
std::vector< G4VEvaporationChannel * > * theChannels
void PrintWelcomeMessage()
G4bool produceSecondaries
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
const G4LorentzVector & GetMomentum() const
std::vector< G4Fragment * > G4FragmentVector
static G4Triton * Triton()
static G4Proton * Proton()
G4double GetGroundStateMass() const
static G4Neutron * Neutron()
static const G4double A[nN]
static G4Deuteron * Deuteron()
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
G4double GetPDGCharge() const
G4FragmentVector * fragmentVector
G4double GetExcitationEnergy() const
G4int GetBaryonNumber() const
CLHEP::HepLorentzVector G4LorentzVector
G4FragmentVector * BreakItUp(const G4Fragment &theNucleus)