52 theProjectileFragmentation(ptr),
53 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
54 projectile3dNucleus(0),target3dNucleus(0)
74 outFile <<
"G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
75 <<
"using G4BinaryCasacde to model the interaction of a light\n"
76 <<
"nucleus with a nucleus.\n"
77 <<
"The lighter of the two nuclei is treated like a set of projectiles\n"
78 <<
"which are transported simultanously through the heavier nucleus.\n";
90 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction starts ######### " <<
G4endl;
115 if( (mom.t()-mom.mag())/
pA < 50*
MeV )
140 G4cerr <<
"G4BinaryLightIonReaction no final state for: " <<
G4endl;
146 G4cerr <<
" Target nucleus (A,Z)=("
147 << (swapped?
pA:
tA) <<
","
149 G4cerr <<
" if frequent, please submit above information as bug report"
180 std::vector<G4ReactionProduct *>::iterator iter;
193 while (std::abs(momentum.e()-pspectators.e()) > 10*
MeV)
201 G4cout <<
"Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" <<
G4endl;
204 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
209 if (++loopcount > 10 )
211 if ( momentum.vect().mag() - momentum.e()> 10*
keV )
213 G4cerr <<
"G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" <<
G4endl;
225 if ( momentum.vect().mag() - momentum.e()> 10*
keV )
228 for (iter=spectators->begin();iter!=spectators->end();iter++)
233 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
239 G4cout <<
"G4BinaryLightIonReaction.cc: mom check: " << momentum
240 <<
" 3.mag "<< momentum.vect().mag() << G4endl
241 <<
" .. pInitialState/pFinalState/spectators " <<
pInitialState <<
" "
244 G4cout <<
"G4BinaryLightIonReaction invalid final state for: " <<
G4endl;
252 G4cout <<
" if frequent, please submit above information as bug report"
254 #ifdef debug_G4BinaryLightIonReaction
256 ed <<
"G4BinaryLightIonreaction: Terminate for above error" <<
G4endl;
276 toZ.rotateZ(-1*mom.phi());
277 toZ.rotateY(-1*mom.theta());
285 G4ReactionProductVector::iterator iter;
286 #ifdef debug_BLIR_result
291 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
293 if((*iter)->GetNewlyAdded())
297 (*iter)->GetTotalEnergy(),
298 (*iter)->GetMomentum() );
300 #ifdef debug_BLIR_result
305 tmp*=toBreit.inverse();
306 tmp.setVect(-tmp.vect());
319 #ifdef debug_BLIR_result
334 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction number ends ######### " <<
G4endl;
346 const int nAttemptScale = 2500;
347 const double ErrLimit = 1.E-6;
352 G4double TotalCollisionMass = TotalCollisionMom.m();
355 for(i = 0; i < Output->size(); i++)
357 SumMom +=
G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
358 SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
362 if (SumMass > TotalCollisionMass)
return FALSE;
363 SumMass = SumMom.m2();
364 if (SumMass < 0)
return FALSE;
365 SumMass = std::sqrt(SumMass);
371 for(i = 0; i < Output->size(); i++)
375 (*Output)[i]->SetMomentum(mom.vect());
376 (*Output)[i]->SetTotalEnergy(mom.e());
386 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
389 for(i = 0; i < Output->size(); i++)
392 HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
393 G4double E = std::sqrt(HadronMom.vect().mag2() +
sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
395 (*Output)[i]->SetMomentum(HadronMom.vect());
396 (*Output)[i]->SetTotalEnergy(HadronMom.e());
400 Scale = TotalCollisionMass/Sum - 1;
402 if (std::abs(Scale) <= ErrLimit
403 || OldScale == Scale)
412 factor=
std::max(1.,std::log(std::abs(OldScale/(OldScale-Scale))));
419 G4cout <<
"G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<
G4endl;
420 G4cout <<
" Scale not unity at end of iteration loop: "<<TotalCollisionMass<<
" "<<Sum<<
" "<<Scale<<
G4endl;
421 G4cout <<
" Increase number of attempts or increase ERRLIMIT"<<
G4endl;
425 Beta = TotalCollisionMom.boostVector();
427 for(i = 0; i < Output->size(); i++)
431 (*Output)[i]->SetMomentum(mom.vect());
432 (*Output)[i]->SetTotalEnergy(mom.e());
458 if (m2Compound <
sqr(mFused) ) {
480 for(
size_t count = 0; count<cascaders->size(); count++)
482 cascaders->operator[](count)->SetNewlyAdded(
true);
516 G4LorentzVector tmpV(0,0,0,0);
517 #ifdef debug_BLIR_finalstate
518 G4LorentzVector pinitial;
520 G4LorentzVector nucleonMom(1./
pA*mom);
521 nucleonMom.setZ(nucleonMom.vect().mag());
531 nucleonPosition +=
pos;
538 initalState->push_back(it1);
539 #ifdef debug_BLIR_finalstate
545 #ifdef debug_BLIR_finalstate
546 if( result && result->size()>0)
548 G4LorentzVector presult;
549 G4ReactionProductVector::iterator iter;
551 for (iter=result->begin(); iter !=result->end(); ++iter)
553 presult +=
G4LorentzVector((*iter)->GetMomentum(),(*iter)->GetTotalEnergy());
557 <<
" final " << presult
562 if( result && result->size()==0)
576 }
while (! result && tryCount< 150);
589 G4double theStatisticalExEnergy = 0;
599 G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
601 theStatisticalExEnergy += deltaE;
604 return theStatisticalExEnergy;
613 for(i=0; i<result->size(); i++)
615 if( (*result)[i]->GetNewlyAdded() )
618 cascaders->push_back((*result)[i]);
622 pspectators +=
G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
623 spectators->push_back((*result)[i]);
673 G4ReactionProductVector::iterator ispectator;
674 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
681 G4ReactionProductVector::iterator ispectator;
682 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
684 (*ispectator)->SetNewlyAdded(
true);
685 cascaders->push_back(*ispectator);
698 G4ReactionProductVector::iterator ii;
701 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
703 (*ii)->SetNewlyAdded(
true);
705 tmp *= boost_fragments;
706 (*ii)->SetMomentum(tmp.vect());
707 (*ii)->SetTotalEnergy(tmp.e());
724 G4cout <<
"G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" <<
G4endl;
730 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
732 cascaders->push_back(*ii);
737 if ( ! EnergyIsCorrect )
743 G4cout <<
"G4BinaryLightIonReaction E/P corrections failed" <<
G4endl;
G4BinaryLightIonReaction(G4VPreCompoundModel *ptr=0)
virtual void ModelDescription(std::ostream &) const
const G4LorentzVector & GetMomentum() const
std::ostringstream G4ExceptionDescription
CLHEP::Hep3Vector G4ThreeVector
G4Fancy3DNucleus * projectile3dNucleus
CLHEP::HepLorentzRotation G4LorentzRotation
virtual G4ReactionProductVector * DeExcite(G4Fragment &aFragment)=0
virtual ~G4BinaryLightIonReaction()
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
G4ReactionProductVector * FuseNucleiAndPrompound(const G4LorentzVector &mom)
void SetNumberOfHoles(G4int valueTot, G4int valueP=0)
G4HadFinalState theResult
virtual const G4ThreeVector & GetPosition() const
void SetProjectilePotential(const G4double aPotential)
G4ReactionProductVector * BreakItUp(const G4Fragment &theInitialState)
virtual const G4ParticleDefinition * GetDefinition() const
void SetStatusChange(G4HadFinalStateStatus aS)
std::vector< G4ReactionProduct * > G4ReactionProductVector
G4BinaryCascade * theModel
G4ExcitationHandler * GetExcitationHandler() const
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
CascadeState SetState(const CascadeState new_state)
const G4ParticleDefinition * GetDefinition() const
G4double GetIonMass(G4int Z, G4int A, G4int L=0, G4int lvl=0) const
G4VPreCompoundModel * theProjectileFragmentation
void SetMomentum(const G4LorentzVector &value)
G4double GetKineticEnergy() const
void Init(G4int theA, G4int theZ)
G4double GetFermiMomentum(G4double density)
G4bool SetLighterAsProjectile(G4LorentzVector &mom, const G4LorentzRotation &toBreit)
G4double GetProjectileExcitation()
void SetNumberOfParticles(G4int value)
G4ReactionProductVector * Interact(G4LorentzVector &mom, const G4LorentzRotation &)
void DeExciteSpectatorNucleus(G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders, G4double theStatisticalExEnergy, G4LorentzVector &momentum)
virtual G4ReactionProductVector * Propagate(G4KineticTrackVector *, G4V3DNucleus *)
const G4LorentzVector & Get4Momentum() const
G4LorentzVector pInitialState
G4LorentzVector Get4Momentum() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4HadronicInteraction * FindModel(const G4String &name)
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
G4double GetTotalEnergy() const
void Init(G4int anA, G4int aZ)
static const G4double factor
G4double GetPDGMass() const
G4double GetDensity(const G4ThreeVector &aPosition) const
static G4ParticleTable * GetParticleTable()
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4HadronicInteractionRegistry * Instance()
void SetZandA_asInt(G4int Znew, G4int Anew)
G4LorentzVector SortResult(G4ReactionProductVector *result, G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders)
G4ThreeVector GetMomentum() const
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
G4double GetOuterRadius()
G4LorentzVector operator()(G4LorentzVector a, G4ReactionProduct *b)
void SetNumberOfCharged(G4int value)
G4LorentzVector pFinalState
const G4VNuclearDensity * GetNuclearDensity() const
static const double eplus
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
G4bool debug_G4BinaryLightIonReactionResults
G4Nucleon * GetNextNucleon()
void SetMomentumChange(const G4ThreeVector &aV)
G4bool EnergyAndMomentumCorrector(G4ReactionProductVector *products, G4LorentzVector &TotalCollisionMom)
static const G4double pos
static const double fermi
G4Fancy3DNucleus * target3dNucleus
G4GLOB_DLL std::ostream G4cerr
G4int GetBaryonNumber() const
CLHEP::HepLorentzVector G4LorentzVector
G4ExcitationHandler * theHandler