51 theProjectileFragmentation(ptr),
52 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
53 projectile3dNucleus(0),target3dNucleus(0)
73 outFile <<
"G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
74 <<
"using G4BinaryCasacde to model the interaction of a light\n"
75 <<
"nucleus with a nucleus.\n"
76 <<
"The lighter of the two nuclei is treated like a set of projectiles\n"
77 <<
"which are transported simultanously through the heavier nucleus.\n";
89 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction starts ######### " <<
G4endl;
114 if( (mom.t()-mom.mag())/
pA < 50*
MeV )
139 G4cerr <<
"G4BinaryLightIonReaction no final state for: " <<
G4endl;
145 G4cerr <<
" Target nucleus (A,Z)=("
146 << (swapped?
pA:
tA) <<
","
148 G4cerr <<
" if frequent, please submit above information as bug report"
178 std::vector<G4ReactionProduct *>::iterator iter;
191 while (std::abs(momentum.e()-pspectators.e()) > 10*
MeV)
199 G4cout <<
"Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" <<
G4endl;
202 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
207 if (++loopcount > 10 )
209 if ( momentum.vect().mag() - momentum.e()> 10*
keV )
211 G4cerr <<
"G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" <<
G4endl;
223 if ( momentum.vect().mag() - momentum.e()> 10*
keV )
226 for (iter=spectators->begin();iter!=spectators->end();iter++)
231 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
237 G4cout <<
"G4BinaryLightIonReaction.cc: mom check: " << momentum
238 <<
" 3.mag "<< momentum.vect().mag() << G4endl
239 <<
" .. pInitialState/pFinalState/spectators " <<
pInitialState <<
" "
242 G4cout <<
"G4BinaryLightIonReaction invalid final state for: " <<
G4endl;
250 G4cout <<
" if frequent, please submit above information as bug report"
268 toZ.rotateZ(-1*mom.phi());
269 toZ.rotateY(-1*mom.theta());
277 G4ReactionProductVector::iterator iter;
278 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
280 if((*iter)->GetNewlyAdded())
284 (*iter)->GetTotalEnergy(),
285 (*iter)->GetMomentum() );
289 tmp*=toBreit.inverse();
290 tmp.setVect(-tmp.vect());
306 #ifdef debug_BLIR_result
308 G4cout <<
" Energy conservation initial/primary/nucleus/final/delta(init-final) "
313 if(getenv(
"BLICDEBUG") )
G4cerr <<
" ######### Binary Light Ion Reaction number ends ######### " <<
G4endl;
325 const int nAttemptScale = 2500;
326 const double ErrLimit = 1.E-6;
331 G4double TotalCollisionMass = TotalCollisionMom.m();
334 for(i = 0; i < Output->size(); i++)
336 SumMom +=
G4LorentzVector((*Output)[i]->GetMomentum(),(*Output)[i]->GetTotalEnergy());
337 SumMass += (*Output)[i]->GetDefinition()->GetPDGMass();
341 if (SumMass > TotalCollisionMass)
return FALSE;
342 SumMass = SumMom.m2();
343 if (SumMass < 0)
return FALSE;
344 SumMass = std::sqrt(SumMass);
350 for(i = 0; i < Output->size(); i++)
354 (*Output)[i]->SetMomentum(mom.vect());
355 (*Output)[i]->SetTotalEnergy(mom.e());
365 for(cAttempt = 0; cAttempt < nAttemptScale; cAttempt++)
368 for(i = 0; i < Output->size(); i++)
371 HadronMom.setVect(HadronMom.vect()+ factor*Scale*HadronMom.vect());
372 G4double E = std::sqrt(HadronMom.vect().mag2() +
sqr((*Output)[i]->GetDefinition()->GetPDGMass()));
374 (*Output)[i]->SetMomentum(HadronMom.vect());
375 (*Output)[i]->SetTotalEnergy(HadronMom.e());
379 Scale = TotalCollisionMass/Sum - 1;
381 if (std::abs(Scale) <= ErrLimit
382 || OldScale == Scale)
391 factor=
std::max(1.,std::log(std::abs(OldScale/(OldScale-Scale))));
398 G4cout <<
"G4G4BinaryLightIonReaction::EnergyAndMomentumCorrector - Warning"<<
G4endl;
399 G4cout <<
" Scale not unity at end of iteration loop: "<<TotalCollisionMass<<
" "<<Sum<<
" "<<Scale<<
G4endl;
400 G4cout <<
" Increase number of attempts or increase ERRLIMIT"<<
G4endl;
404 Beta = TotalCollisionMom.boostVector();
406 for(i = 0; i < Output->size(); i++)
410 (*Output)[i]->SetMomentum(mom.vect());
411 (*Output)[i]->SetTotalEnergy(mom.e());
437 if (m2Compound <
sqr(mFused) ) {
460 for(
size_t count = 0; count<cascaders->size(); count++)
462 cascaders->operator[](count)->SetNewlyAdded(
true);
496 G4LorentzVector tmpV(0,0,0,0);
497 G4LorentzVector nucleonMom(1./
pA*mom);
498 nucleonMom.setZ(nucleonMom.vect().mag());
508 nucleonPosition +=
pos;
515 initalState->push_back(it1);
519 if( result && result->size()==0)
533 }
while (! result && tryCount< 150);
547 G4double theStatisticalExEnergy = 0;
563 G4double localFermiEnergy = std::sqrt(nucMass*nucMass + localPfermi*localPfermi) - nucMass;
565 theStatisticalExEnergy += deltaE;
568 return theStatisticalExEnergy;
577 for(i=0; i<result->size(); i++)
579 if( (*result)[i]->GetNewlyAdded() )
582 cascaders->push_back((*result)[i]);
586 pspectators +=
G4LorentzVector( (*result)[i]->GetMomentum(), (*result)[i]->GetTotalEnergy() );
587 spectators->push_back((*result)[i]);
636 G4ReactionProductVector::iterator ispectator;
637 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
644 G4ReactionProductVector::iterator ispectator;
645 for (ispectator=spectators->begin();ispectator!=spectators->end();ispectator++)
647 (*ispectator)->SetNewlyAdded(
true);
648 cascaders->push_back(*ispectator);
649 pFragments+=
G4LorentzVector((*ispectator)->GetMomentum(),(*ispectator)->GetTotalEnergy());
660 G4ReactionProductVector::iterator ii;
663 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
665 (*ii)->SetNewlyAdded(
true);
667 tmp *= boost_fragments;
668 (*ii)->SetMomentum(tmp.vect());
669 (*ii)->SetTotalEnergy(tmp.e());
685 G4cout <<
"G4BinaryLightIonReaction E/P correction for nucleus failed, will try to correct overall" <<
G4endl;
691 for(ii=proFrag->begin(); ii!=proFrag->end(); ii++)
693 cascaders->push_back(*ii);
698 if ( ! EnergyIsCorrect )
704 G4cout <<
"G4BinaryLightIonReaction E/P corrections failed" <<
G4endl;
G4BinaryLightIonReaction(G4VPreCompoundModel *ptr=0)
virtual void ModelDescription(std::ostream &) const
const G4LorentzVector & GetMomentum() const
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) const
void SetA(G4double value)
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
G4HadronicInteraction * FindModel(const G4String &name)
void Set4Momentum(const G4LorentzVector &momentum)
void SetEnergyChange(G4double anEnergy)
G4double GetTotalEnergy() const
void Init(G4int anA, G4int aZ)
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()
G4LorentzVector SortResult(G4ReactionProductVector *result, G4ReactionProductVector *spectators, G4ReactionProductVector *cascaders)
G4ThreeVector GetMomentum() const
virtual G4ParticleDefinition * GetDefinition() const
void SetZ(G4double value)
G4double GetOuterRadius()
G4LorentzVector operator()(G4LorentzVector a, G4ReactionProduct *b)
void SetNumberOfCharged(G4int value)
G4LorentzVector pFinalState
const G4VNuclearDensity * GetNuclearDensity() const
static const double eplus
G4double GetPDGCharge() const
G4bool debug_G4BinaryLightIonReactionResults
G4Nucleon * GetNextNucleon()
void SetMomentumChange(const G4ThreeVector &aV)
G4bool EnergyAndMomentumCorrector(G4ReactionProductVector *products, G4LorentzVector &TotalCollisionMom)
G4int GetNumberOfSecondaries() const
static const G4double pos
void AddSecondary(G4DynamicParticle *aP)
static const double fermi
G4Fancy3DNucleus * target3dNucleus
G4GLOB_DLL std::ostream G4cerr
G4int GetBaryonNumber() const
G4double GetTotalEnergy() const
CLHEP::HepLorentzVector G4LorentzVector
G4ExcitationHandler * theHandler