47 #ifdef G4MULTITHREADED
59 theProjectileFragmentation(ptr),
60 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
61 projectile3dNucleus(0),target3dNucleus(0)
73 #ifdef G4MULTITHREADED
78 #ifdef G4MULTITHREADED
93 outFile <<
"G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
94 <<
"using G4BinaryCasacde to model the interaction of a light\n"
95 <<
"nucleus with a nucleus.\n"
96 <<
"The lighter of the two nuclei is treated like a set of projectiles\n"
97 <<
"which are transported simultanously through the heavier nucleus.\n";
109 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction starts ######### " <<
G4endl;
134 if( (mom.t()-mom.mag())/
pA < 50*
MeV )
159 G4cerr <<
"G4BinaryLightIonReaction no final state for: " <<
G4endl;
165 G4cerr <<
" Target nucleus (A,Z)=("
166 << (swapped?
pA:
tA) <<
","
168 G4cerr <<
" if frequent, please submit above information as bug report"
200 std::vector<G4ReactionProduct *>::iterator iter;
213 while (std::abs(momentum.
e()-pspectators.
e()) > 10*
MeV)
222 G4cout <<
"Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" <<
G4endl;
225 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
230 if (++loopcount > 10 )
232 if ( momentum.
vect().
mag() - momentum.
e()> 10*
keV )
234 G4cerr <<
"G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" <<
G4endl;
246 if ( momentum.
vect().
mag() - momentum.
e()< 10*
keV )
252 for (iter=spectators->begin();iter!=spectators->end();iter++)
257 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
263 G4cout <<
"G4BinaryLightIonReaction.cc: mom check: " << momentum
264 <<
" 3.mag "<< momentum.
vect().
mag() << G4endl
265 <<
" .. pInitialState/pFinalState/spectators " <<
pInitialState <<
" "
268 G4cout <<
"G4BinaryLightIonReaction invalid final state for: " <<
G4endl;
276 G4cout <<
" if frequent, please submit above information as bug report"
278 #ifdef debug_G4BinaryLightIonReaction
280 ed <<
"G4BinaryLightIonreaction: Terminate for above error" <<
G4endl;
306 G4ReactionProductVector::iterator iter;
307 #ifdef debug_BLIR_result
312 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
314 if((*iter)->GetNewlyAdded())
318 (*iter)->GetTotalEnergy(),
319 (*iter)->GetMomentum() );
321 #ifdef debug_BLIR_result
326 tmp*=toBreit.inverse();
327 tmp.setVect(-tmp.vect());
334 aNew.
SetTime(timePrimary + time);
346 #ifdef debug_BLIR_result
361 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction number ends ######### " <<
G4endl;
373 const int nAttemptScale = 2500;
374 const double ErrLimit = 1.E-6;
379 G4double TotalCollisionMass = TotalCollisionMom.
m();
382 for(i = 0; i < Output->size(); i++)
384 SumMom +=
G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
385 SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
389 if (SumMass > TotalCollisionMass)
return FALSE;
390 SumMass = SumMom.
m2();
391 if (SumMass < 0)
return FALSE;
392 SumMass = std::sqrt(SumMass);
398 for(i = 0; i < Output->size(); i++)
402 (*Output)[i]->SetMomentum(mom.
vect());
403 (*Output)[i]->SetTotalEnergy(mom.
e());
413 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
416 for(i = 0; i < Output->size(); i++)
420 G4double E = std::sqrt(HadronMom.
vect().
mag2() +
sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
422 (*Output)[i]->SetMomentum(HadronMom.
vect());
423 (*Output)[i]->SetTotalEnergy(HadronMom.
e());
427 Scale = TotalCollisionMass/Sum - 1;
429 if (std::abs(Scale) <= ErrLimit
430 || OldScale ==
Scale)
439 factor=
std::max(1.,
G4Log(std::abs(OldScale/(OldScale-Scale))));
446 G4cout <<
"G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<
G4endl;
447 G4cout <<
" Scale not unity at end of iteration loop: "<<TotalCollisionMass<<
" "<<Sum<<
" "<<Scale<<
G4endl;
448 G4cout <<
" Increase number of attempts or increase ERRLIMIT"<<
G4endl;
454 for(i = 0; i < Output->size(); i++)
458 (*Output)[i]->SetMomentum(mom.
vect());
459 (*Output)[i]->SetTotalEnergy(mom.
e());
485 if (m2Compound <
sqr(mFused) ) {
507 for(
size_t count = 0; count<cascaders->size(); count++)
509 cascaders->operator[](count)->SetNewlyAdded(
true);
543 G4LorentzVector tmpV(0,0,0,0);
544 #ifdef debug_BLIR_finalstate
545 G4LorentzVector pinitial;
547 G4LorentzVector nucleonMom(1./
pA*mom);
558 nucleonPosition +=
pos;
565 initalState->push_back(it1);
566 #ifdef debug_BLIR_finalstate
572 #ifdef debug_BLIR_finalstate
573 if( result && result->size()>0)
575 G4LorentzVector presult;
576 G4ReactionProductVector::iterator iter;
578 for (iter=result->begin(); iter !=result->end(); ++iter)
580 presult +=
G4LorentzVector((*iter)->GetMomentum(),(*iter)->GetTotalEnergy());
584 <<
" final " << presult
589 if( result && result->size()==0)
603 }
while (! result && tryCount< 150);
611 G4double theStatisticalExEnergy = 0;
621 G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
623 theStatisticalExEnergy += deltaE;
626 return theStatisticalExEnergy;
635 for(i=0; i<result->size(); i++)
637 if( (*result)[i]->GetNewlyAdded() )
640 cascaders->push_back((*result)[i]);
644 pspectators +=
G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
645 spectators->push_back((*result)[i]);
695 G4ReactionProductVector::iterator ispectator;
696 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
703 G4ReactionProductVector::iterator ispectator;
704 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
706 (*ispectator)->SetNewlyAdded(
true);
707 cascaders->push_back(*ispectator);
720 G4ReactionProductVector::iterator ii;
723 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
725 (*ii)->SetNewlyAdded(
true);
727 tmp *= boost_fragments;
728 (*ii)->SetMomentum(tmp.vect());
729 (*ii)->SetTotalEnergy(tmp.e());
746 G4cout <<
"G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" <<
G4endl;
752 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
754 cascaders->push_back(*ii);
759 if ( ! EnergyIsCorrect )
765 G4cout <<
"G4BinaryLightIonReaction E/P corrections failed" <<
G4endl;
void SetProjectilePotential(const G4double aPotential)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void Init(G4int anA, G4int aZ)
CLHEP::Hep3Vector G4ThreeVector
G4LorentzVector pFinalState
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
std::ostringstream G4ExceptionDescription
const G4LorentzVector & Get4Momentum() const
std::vector< ExP01TrackerHit * > a
static const G4double pos
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
G4Fancy3DNucleus * target3dNucleus
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
G4double GetFermiMomentum(G4double density)
G4LorentzVector operator()(G4LorentzVector a, G4ReactionProduct *b)
void SetMomentumChange(const G4ThreeVector &aV)
static constexpr double keV
G4Nucleon * GetNextNucleon()
G4ExcitationHandler * GetExcitationHandler() const
G4double GetDensity(const G4ThreeVector &aPosition) const
void SetZandA_asInt(G4int Znew, G4int Anew)
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
G4IonTable * GetIonTable() const
G4double GetPDGCharge() const
G4ReactionProductVector * Interact(G4LorentzVector &mom, const G4LorentzRotation &)
void setVect(const Hep3Vector &)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4bool debug_G4BinaryLightIonReactionResults
G4double GetPDGMass() const
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *, G4V3DNucleus *)
const G4LorentzVector & GetMomentum() const
G4double G4Log(G4double x)
virtual ~G4BinaryLightIonReaction()
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
const G4ThreeVector & GetPosition() const
void SetEnergyChange(G4double anEnergy)
CLHEP::HepLorentzRotation G4LorentzRotation
G4double GetProjectileExcitation()
G4double GetTotalEnergy() const
static constexpr double fermi
#define G4MUTEX_INITIALIZER
G4double GetGlobalTime() const
void DeExciteSpectatorNucleus(G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders, G4double theStatisticalExEnergy, G4LorentzVector &momentum)
CascadeState SetState(const CascadeState new_state)
HepLorentzRotation & rotateZ(double delta)
G4BinaryCascade * theModel
G4int GetBaryonNumber() const
static G4HadronicInteractionRegistry * Instance()
G4double GetOuterRadius()
void SetCreatorModelType(G4int idx)
void Init(G4int theA, G4int theZ)
void Set4Momentum(const G4LorentzVector &momentum)
G4LorentzVector Get4Momentum() const
HepLorentzRotation inverse() const
void SetMomentum(const G4LorentzVector &value)
G4ThreeVector GetMomentum() const
void SetNumberOfCharged(G4int value)
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4GLOB_DLL std::ostream G4cerr
#define G4MUTEXUNLOCK(mutex)
const G4LorentzVector & Get4Momentum() const
G4double G4ParticleHPJENDLHEData::G4double result
const G4VNuclearDensity * GetNuclearDensity() const
G4double GetKineticEnergy() const
static constexpr double eplus
virtual const G4ParticleDefinition * GetDefinition() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
const G4ParticleDefinition * GetDefinition() const
G4VPreCompoundModel * theProjectileFragmentation
G4BinaryLightIonReaction(G4VPreCompoundModel *ptr=0)
virtual void ModelDescription(std::ostream &) const
#define G4MUTEXLOCK(mutex)
G4HadFinalState theResult
void SetNumberOfParticles(G4int value)
G4GLOB_DLL std::ostream G4cout
Hep3Vector boostVector() const
G4Fancy3DNucleus * projectile3dNucleus
void SetTime(G4double aT)
CLHEP::HepLorentzVector G4LorentzVector
HepLorentzRotation & rotateY(double delta)
G4ReactionProductVector * FuseNucleiAndPrompound(const G4LorentzVector &mom)
G4HadronicInteraction * FindModel(const G4String &name)
G4LorentzVector pInitialState
simulatedPeak Scale(1/simulationNormalisationFactor)
G4bool EnergyAndMomentumCorrector(G4ReactionProductVector *products, G4LorentzVector &TotalCollisionMom)
G4LorentzVector SortResult(G4ReactionProductVector *result, G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders)
G4ExcitationHandler * theHandler
G4bool SetLighterAsProjectile(G4LorentzVector &mom, const G4LorentzRotation &toBreit)
void SetStatusChange(G4HadFinalStateStatus aS)
static G4int Register(const G4String &)