Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions | List of all members
G4BinaryLightIonReaction Class Reference

#include <G4BinaryLightIonReaction.hh>

Inheritance diagram for G4BinaryLightIonReaction:
Inheritance graph
[legend]
Collaboration diagram for G4BinaryLightIonReaction:
Collaboration graph
[legend]

Public Member Functions

 G4BinaryLightIonReaction (G4VPreCompoundModel *ptr=0)
 
virtual ~G4BinaryLightIonReaction ()
 
G4HadFinalState * ApplyYourself (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 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 G4String & GetModelName () 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, G4double > 		GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double > 		GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 
virtual void BuildPhysicsTable (const G4ParticleDefinition &)
 
virtual void InitialiseModel ()
 

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 ( G4VPreCompoundModel *  ptr = 0)

Definition at line 51 of file G4BinaryLightIonReaction.cc.

52 : G4HadronicInteraction("Binary Light Ion Cascade"),
53  theProjectileFragmentation(ptr),
54  pA(0),pZ(0), tA(0),tZ(0),spectatorA(0),spectatorZ(0),
55  projectile3dNucleus(0),target3dNucleus(0)
56 {
57  if(!ptr) {
58  G4HadronicInteraction* p =
59  G4HadronicInteractionRegistry::Instance()->FindModel("PRECO");
60  G4VPreCompoundModel* pre = static_cast<G4VPreCompoundModel*>(p);
61  if(!pre) { pre = new G4PreCompoundModel(); }
62  theProjectileFragmentation = pre;
63  }
64  theModel = new G4BinaryCascade(theProjectileFragmentation);
65  theHandler = theProjectileFragmentation->GetExcitationHandler();
66 
67  debug_G4BinaryLightIonReactionResults=getenv("debug_G4BinaryLightIonReactionResults")!=0;
68 }
p
const char * p
Definition: xmltok.h:285
G4VPreCompoundModel::GetExcitationHandler
G4ExcitationHandler * GetExcitationHandler() const
Definition: G4VPreCompoundModel.hh:93
G4HadronicInteraction::G4HadronicInteraction
G4HadronicInteraction(const G4String &modelName="HadronicModel")
Definition: G4HadronicInteraction.cc:42
G4HadronicInteractionRegistry::FindModel
G4HadronicInteraction * FindModel(const G4String &name)
Definition: G4HadronicInteractionRegistry.cc:117
G4HadronicInteractionRegistry::Instance
static G4HadronicInteractionRegistry * Instance()
Definition: G4HadronicInteractionRegistry.cc:36

Here is the call graph for this function:

G4BinaryLightIonReaction::~G4BinaryLightIonReaction ( )
virtual

Definition at line 70 of file G4BinaryLightIonReaction.cc.

71 {}

Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 89 of file G4BinaryLightIonReaction.cc.

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

Here is the call graph for this function:

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

Reimplemented from G4HadronicInteraction.

Definition at line 73 of file G4BinaryLightIonReaction.cc.

74 {
75  outFile << "G4Binary Light Ion Cascade is an intra-nuclear cascade model\n"
76  << "using G4BinaryCasacde to model the interaction of a light\n"
77  << "nucleus with a nucleus.\n"
78  << "The lighter of the two nuclei is treated like a set of projectiles\n"
79  << "which are transported simultanously through the heavier nucleus.\n";
80 }
void G4BinaryLightIonReaction::SetDeExcitation ( G4ExcitationHandler *  ptr)
inline

Definition at line 74 of file G4BinaryLightIonReaction.hh.

75 {
76  theProjectileFragmentation->SetExcitationHandler(ptr);
77  theHandler = ptr;
78 }
G4VPreCompoundModel::SetExcitationHandler
void SetExcitationHandler(G4ExcitationHandler *ptr)
Definition: G4VPreCompoundModel.hh:88

Here is the call graph for this function:

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 }
G4VPreCompoundModel::GetExcitationHandler
G4ExcitationHandler * GetExcitationHandler() const
Definition: G4VPreCompoundModel.hh:93

Here is the call graph for this function:

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