Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4BinaryLightIonReaction Class Reference

#include <G4BinaryLightIonReaction.hh>

+ Inheritance diagram for G4BinaryLightIonReaction:

Public Member Functions

 G4BinaryLightIonReaction (G4VPreCompoundModel *ptr=0)
 
virtual ~G4BinaryLightIonReaction ()
 
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &theNucleus)
 
void SetPrecompound (G4VPreCompoundModel *ptr)
 
void SetDeExcitation (G4ExcitationHandler *ptr)
 
virtual void ModelDescription (std::ostream &) const
 
- Public Member Functions inherited from G4HadronicInteraction
 G4HadronicInteraction (const G4String &modelName="HadronicModel")
 
virtual ~G4HadronicInteraction ()
 
virtual G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
virtual G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
 
G4double GetMinEnergy () const
 
G4double GetMinEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMinEnergy (G4double anEnergy)
 
void SetMinEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMinEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4double GetMaxEnergy () const
 
G4double GetMaxEnergy (const G4Material *aMaterial, const G4Element *anElement) const
 
void SetMaxEnergy (const G4double anEnergy)
 
void SetMaxEnergy (G4double anEnergy, const G4Element *anElement)
 
void SetMaxEnergy (G4double anEnergy, const G4Material *aMaterial)
 
G4int GetVerboseLevel () const
 
void SetVerboseLevel (G4int value)
 
const G4StringGetModelName () const
 
void DeActivateFor (const G4Material *aMaterial)
 
void ActivateFor (const G4Material *aMaterial)
 
void DeActivateFor (const G4Element *anElement)
 
void ActivateFor (const G4Element *anElement)
 
G4bool IsBlocked (const G4Material *aMaterial) const
 
G4bool IsBlocked (const G4Element *anElement) const
 
void SetRecoilEnergyThreshold (G4double val)
 
G4double GetRecoilEnergyThreshold () const
 
virtual const std::pair< G4double, G4doubleGetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double, G4doubleGetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void ModelDescription (std::ostream &outFile) const
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 
 G4HadronicInteraction (const G4HadronicInteraction &right)=delete
 
const G4HadronicInteractionoperator= (const G4HadronicInteraction &right)=delete
 
G4bool operator== (const G4HadronicInteraction &right) const =delete
 
G4bool operator!= (const G4HadronicInteraction &right) const =delete
 

Additional Inherited Members

- Protected Member Functions inherited from G4HadronicInteraction
void SetModelName (const G4String &nam)
 
G4bool IsBlocked () const
 
void Block ()
 
- Protected Attributes inherited from G4HadronicInteraction
G4HadFinalState theParticleChange
 
G4int verboseLevel
 
G4double theMinEnergy
 
G4double theMaxEnergy
 
G4bool isBlocked
 

Detailed Description

Definition at line 34 of file G4BinaryLightIonReaction.hh.

Constructor & Destructor Documentation

◆ G4BinaryLightIonReaction()

G4BinaryLightIonReaction::G4BinaryLightIonReaction ( G4VPreCompoundModel ptr = 0)

Definition at line 54 of file G4BinaryLightIonReaction.cc.

55: G4HadronicInteraction("Binary Light Ion Cascade"),
56 theProjectileFragmentation(ptr),
57 pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
58 projectile3dNucleus(0),target3dNucleus(0)
59{
60 if(!ptr) {
63 G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
64 if(!pre) { pre = new G4PreCompoundModel(); }
65 theProjectileFragmentation = pre;
66 }
67 theModel = new G4BinaryCascade(theProjectileFragmentation);
68 theHandler = theProjectileFragmentation->GetExcitationHandler();
69 theBLIR_ID = G4PhysicsModelCatalog::GetModelID("model_G4BinaryLightIonReaction");
70 debug_G4BinaryLightIonReactionResults=std::getenv("debug_G4BinaryLightIonReactionResults")!=0;
71}
G4HadronicInteraction * FindModel(const G4String &name)
static G4HadronicInteractionRegistry * Instance()
static G4int GetModelID(const G4int modelIndex)
G4ExcitationHandler * GetExcitationHandler() const

◆ ~G4BinaryLightIonReaction()

G4BinaryLightIonReaction::~G4BinaryLightIonReaction ( )
virtual

Definition at line 73 of file G4BinaryLightIonReaction.cc.

74{}

Member Function Documentation

◆ ApplyYourself()

G4HadFinalState * G4BinaryLightIonReaction::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus theNucleus 
)
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 91 of file G4BinaryLightIonReaction.cc.

93{
94 if(std::getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction starts ######### " << G4endl;
95 G4ping debug("debug_G4BinaryLightIonReaction");
96 pA=aTrack.GetDefinition()->GetBaryonNumber();
97 pZ=G4lrint(aTrack.GetDefinition()->GetPDGCharge()/eplus);
98 tA=targetNucleus.GetA_asInt();
99 tZ=targetNucleus.GetZ_asInt();
100 G4double timePrimary = aTrack.GetGlobalTime();
101 G4LorentzVector mom(aTrack.Get4Momentum());
102 //G4cout << "proj mom : " << mom << G4endl;
103 G4LorentzRotation toBreit(mom.boostVector());
104
105 G4bool swapped=SetLighterAsProjectile(mom, toBreit);
106 //G4cout << "after swap, swapped? / mom " << swapped << " / " << mom <<G4endl;
107 G4ReactionProductVector * result = 0;
108 G4ReactionProductVector * cascaders=0; //new G4ReactionProductVector;
109// G4double m_nucl(0); // to check energy balance
110
111
112 // G4double m1=G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(pZ,pA);
113 // G4cout << "Entering the decision point "
114 // << (mom.t()-mom.mag())/pA << " "
115 // << pA<<" "<< pZ<<" "
116 // << tA<<" "<< tZ<<G4endl
117 // << " "<<mom.t()-mom.mag()<<" "
118 // << mom.t()- m1<<G4endl;
119 if( (mom.t()-mom.mag())/pA < 50*MeV )
120 {
121 // G4cout << "Using pre-compound only, E= "<<mom.t()-mom.mag()<<G4endl;
122 // m_nucl = mom.mag();
123 cascaders=FuseNucleiAndPrompound(mom);
124 if( !cascaders )
125 {
126
127 // abort!! happens for too low energy for nuclei to fuse
128
129 theResult.Clear();
130 theResult.SetStatusChange(isAlive);
131 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
132 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
133 return &theResult;
134 }
135 }
136 else
137 {
138 result=Interact(mom,toBreit);
139
140 if(! result )
141 {
142 // abort!!
143
144 G4cerr << "G4BinaryLightIonReaction no final state for: " << G4endl;
145 G4cerr << " Primary " << aTrack.GetDefinition()
146 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
147 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
148 << ", kinetic energy " << aTrack.GetKineticEnergy()
149 << G4endl;
150 G4cerr << " Target nucleus (A,Z)=("
151 << (swapped?pA:tA) << ","
152 << (swapped?pZ:tZ) << ")" << G4endl;
153 G4cerr << " if frequent, please submit above information as bug report"
154 << G4endl << G4endl;
155
156 theResult.Clear();
157 theResult.SetStatusChange(isAlive);
158 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
159 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
160 return &theResult;
161 }
162
163 // Calculate excitation energy,
164 G4double theStatisticalExEnergy = GetProjectileExcitation();
165
166
167 pInitialState = mom;
168 //G4cout << "BLIC: pInitialState from aTrack : " << pInitialState;
169 pInitialState.setT(pInitialState.getT() +
171 //G4cout << "BLIC: target nucleus added : " << pInitialState << G4endl;
172
173 delete target3dNucleus;target3dNucleus=0;
174 delete projectile3dNucleus;projectile3dNucleus=0;
175
177
178 cascaders = new G4ReactionProductVector;
179
180 G4LorentzVector pspectators=SortResult(result,spectators,cascaders);
181 // this also sets spectatorA and spectatorZ
182
183 // pFinalState=std::accumulate(cascaders->begin(),cascaders->end(),pFinalState,ReactionProduct4Mom);
184
185 std::vector<G4ReactionProduct *>::iterator iter;
186
187 // G4cout << "pInitialState, pFinalState / pspectators"<< pInitialState << " / " << pFinalState << " / " << pspectators << G4endl;
188 // if ( spectA-spectatorA !=0 || spectZ-spectatorZ !=0)
189 // {
190 // G4cout << "spect Nucl != spectators: nucl a,z; spect a,z" <<
191 // spectatorA <<" "<< spectatorZ <<" ; " << spectA <<" "<< spectZ << G4endl;
192 // }
193 delete result;
194 result=0;
195 G4LorentzVector momentum(pInitialState-pFinalState);
196 G4int loopcount(0);
197 //G4cout << "BLIC: momentum, pspectators : " << momentum << " / " << pspectators << G4endl;
198 while (std::abs(momentum.e()-pspectators.e()) > 10*MeV) /* Loop checking, 31.08.2015, G.Folger */
199 // see if on loopcount
200 {
201 G4LorentzVector pCorrect(pInitialState-pspectators);
202 //G4cout << "BLIC:: BIC nonconservation? (pInitialState-pFinalState) / spectators :" << momentum << " / " << pspectators << "pCorrect "<< pCorrect<< G4endl;
203 // Correct outgoing casacde particles.... to have momentum of (initial state - spectators)
204 G4bool EnergyIsCorrect=EnergyAndMomentumCorrector(cascaders, pCorrect);
205 if ( ! EnergyIsCorrect && debug_G4BinaryLightIonReactionResults)
206 {
207 G4cout << "Warning - G4BinaryLightIonReaction E/P correction for cascaders failed" << G4endl;
208 }
209 pFinalState=G4LorentzVector(0,0,0,0);
210 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
211 {
212 pFinalState += G4LorentzVector( (*iter)->GetMomentum(), (*iter)->GetTotalEnergy() );
213 }
214 momentum=pInitialState-pFinalState;
215 if (++loopcount > 10 )
216 {
217 if ( momentum.vect().mag() - momentum.e()> 10*keV )
218 {
219 G4cerr << "G4BinaryLightIonReaction.cc: Cannot correct 4-momentum of cascade particles" << G4endl;
220 throw G4HadronicException(__FILE__, __LINE__, "G4BinaryCasacde::ApplyCollision()");
221 } else {
222 break;
223 }
224
225 }
226 }
227
228 if (spectatorA > 0 )
229 {
230 // check spectator momentum
231 if ( momentum.vect().mag() - momentum.e()< 10*keV )
232 {
233 // DeExciteSpectatorNucleus() also handles also case of A=1, Z=0,1
234 DeExciteSpectatorNucleus(spectators, cascaders, theStatisticalExEnergy, momentum);
235
236 } else { // momentum non-conservation --> fail
237 for (iter=spectators->begin();iter!=spectators->end();iter++)
238 {
239 delete *iter;
240 }
241 delete spectators;
242 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
243 {
244 delete *iter;
245 }
246 delete cascaders;
247
248 G4cout << "G4BinaryLightIonReaction.cc: mom check: " << momentum
249 << " 3.mag "<< momentum.vect().mag() << G4endl
250 << " .. pInitialState/pFinalState/spectators " << pInitialState <<" "
251 << pFinalState << " " << pspectators << G4endl
252 << " .. A,Z " << spectatorA <<" "<< spectatorZ << G4endl;
253 G4cout << "G4BinaryLightIonReaction invalid final state for: " << G4endl;
254 G4cout << " Primary " << aTrack.GetDefinition()
255 << ", (A,Z)=(" << aTrack.GetDefinition()->GetBaryonNumber()
256 << "," << aTrack.GetDefinition()->GetPDGCharge()/eplus << ") "
257 << ", kinetic energy " << aTrack.GetKineticEnergy()
258 << G4endl;
259 G4cout << " Target nucleus (A,Z)=(" << targetNucleus.GetA_asInt()
260 << "," << targetNucleus.GetZ_asInt() << ")" << G4endl;
261 G4cout << " if frequent, please submit above information as bug report"
262 << G4endl << G4endl;
263#ifdef debug_G4BinaryLightIonReaction
265 ed << "G4BinaryLightIonreaction: Terminate for above error" << G4endl;
266 G4Exception("G4BinaryLightIonreaction::ApplyYourSelf()", "BLIC001", FatalException,
267 ed);
268
269#endif
270 theResult.Clear();
271 theResult.SetStatusChange(isAlive);
272 theResult.SetEnergyChange(aTrack.GetKineticEnergy());
273 theResult.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
274 return &theResult;
275 }
276 } else { // no spectators
277 delete spectators;
278 }
279 }
280 // Rotate to lab
282 toZ.rotateZ(-1*mom.phi());
283 toZ.rotateY(-1*mom.theta());
284 G4LorentzRotation toLab(toZ.inverse());
285
286 // Fill the particle change, while rotating. Boost from projectile breit-frame in case we swapped.
287 // theResult.Clear();
288 theResult.Clear();
289 theResult.SetStatusChange(stopAndKill);
290 G4LorentzVector ptot(0);
291 G4ReactionProductVector::iterator iter;
292 #ifdef debug_BLIR_result
293 G4LorentzVector p_raw;
294 #endif
295 //G4int i=0;
296
297 for(iter=cascaders->begin(); iter!=cascaders->end(); iter++)
298 {
299 if((*iter)->GetNewlyAdded())
300 {
301 G4DynamicParticle * aNewDP =
302 new G4DynamicParticle((*iter)->GetDefinition(),
303 (*iter)->GetTotalEnergy(),
304 (*iter)->GetMomentum() );
305 G4LorentzVector tmp = aNewDP->Get4Momentum();
306 #ifdef debug_BLIR_result
307 p_raw+= tmp;
308 #endif
309 if(swapped)
310 {
311 tmp*=toBreit.inverse();
312 tmp.setVect(-tmp.vect());
313 }
314 tmp *= toLab;
315 aNewDP->Set4Momentum(tmp);
316 G4HadSecondary aNew = G4HadSecondary(aNewDP);
317 G4double time = 0; //(*iter)->GetCreationTime();
318 //if(time < 0.0) { time = 0.0; }
319 aNew.SetTime(timePrimary + time);
320 //aNew.SetCreatorModelID((*iter)->GetCreatorModelID()); //AR-02Aug2021 : For some reasons, it does NOT work!
321 aNew.SetCreatorModelID(theBLIR_ID);
322
323 theResult.AddSecondary(aNew);
324 ptot += tmp;
325 //G4cout << "BLIC: Secondary " << aNew->GetDefinition()->GetParticleName()
326 // <<" "<< aNew->GetMomentum()<<" "<< aNew->GetTotalEnergy() << G4endl;
327 }
328 delete *iter;
329 }
330 delete cascaders;
331
332#ifdef debug_BLIR_result
333 //G4cout << "Result analysis, secondaries " << theResult.GetNumberOfSecondaries() << G4endl;
334 //G4cout << "p_tot_raw " << p_raw << " sum p final " << ptot << G4endl;
336 GetIonMass(targetNucleus.GetZ_asInt(),targetNucleus.GetA_asInt());
337 // delete? tZ=targetNucleus.GetZ_asInt();
338
339 //G4cout << "BLIC Energy conservation initial/primary/nucleus/final/delta(init-final) "
340 // << aTrack.GetTotalEnergy() + m_nucl <<" "<< aTrack.GetTotalEnergy() <<" "<< m_nucl <<" "<<ptot.e()
341 // <<" "<< aTrack.GetTotalEnergy() + m_nucl - ptot.e() << G4endl;
342 G4cout << "BLIC momentum conservation " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl)
343 << " ptot " << ptot << " delta " << aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot
344 << " 3mom.mag() " << (aTrack.Get4Momentum()+ G4LorentzVector(m_nucl) - ptot).vect().mag() << G4endl;
345#endif
346
347 if(std::getenv("BLICDEBUG") ) G4cerr << " ######### Binary Light Ion Reaction number ends ######### " << G4endl;
348
349 return &theResult;
350}
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
@ isAlive
@ stopAndKill
CLHEP::HepLorentzVector G4LorentzVector
std::vector< G4ReactionProduct * > G4ReactionProductVector
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
Hep3Vector unit() const
HepLorentzRotation & rotateY(double delta)
HepLorentzRotation & rotateZ(double delta)
HepLorentzRotation inverse() const
Hep3Vector vect() const
void setVect(const Hep3Vector &)
G4LorentzVector Get4Momentum() const
void Set4Momentum(const G4LorentzVector &momentum)
void SetStatusChange(G4HadFinalStateStatus aS)
void AddSecondary(G4DynamicParticle *aP, G4int mod=-1)
void SetEnergyChange(G4double anEnergy)
void SetMomentumChange(const G4ThreeVector &aV)
const G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
const G4LorentzVector & Get4Momentum() const
G4double GetGlobalTime() const
void SetTime(G4double aT)
void SetCreatorModelID(G4int id)
G4double GetIonMass(G4int Z, G4int A, G4int nL=0, G4int lvl=0) const
Definition: G4IonTable.cc:1517
G4double GetPDGCharge() const
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
Definition: G4ping.hh:35
int G4lrint(double ad)
Definition: templates.hh:134

◆ ModelDescription()

void G4BinaryLightIonReaction::ModelDescription ( std::ostream &  outFile) const
virtual

Reimplemented from G4HadronicInteraction.

Definition at line 76 of file G4BinaryLightIonReaction.cc.

77{
78 outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
79 << "using G4BinaryCasacde to model the interaction of a light\n"
80 << "nucleus with a nucleus.\n"
81 << "The lighter of the two nuclei is treated like a set of projectiles\n"
82 << "which are transported simultaneously through the heavier nucleus.\n";
83}

◆ SetDeExcitation()

void G4BinaryLightIonReaction::SetDeExcitation ( G4ExcitationHandler ptr)
inline

Definition at line 74 of file G4BinaryLightIonReaction.hh.

75{
76 theProjectileFragmentation->SetExcitationHandler(ptr);
77 theHandler = ptr;
78}
void SetExcitationHandler(G4ExcitationHandler *ptr)

◆ SetPrecompound()

void G4BinaryLightIonReaction::SetPrecompound ( G4VPreCompoundModel ptr)
inline

Definition at line 69 of file G4BinaryLightIonReaction.hh.

70{
71 if(ptr) { theProjectileFragmentation = ptr; }
72 theHandler = theProjectileFragmentation->GetExcitationHandler();
73}

The documentation for this class was generated from the following files: