Geant4.10
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Public Member Functions
G4BinaryLightIonReaction Class Reference

#include <G4BinaryLightIonReaction.hh>

Inheritance diagram for G4BinaryLightIonReaction:
G4HadronicInteraction

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 G4double SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A)
 
virtual G4bool IsApplicable (const G4HadProjectile &, G4Nucleus &)
 
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)
 
const G4HadronicInteractionGetMyPointer () const
 
virtual G4int GetVerboseLevel () const
 
virtual 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
 
G4bool operator== (const G4HadronicInteraction &right) const
 
G4bool operator!= (const G4HadronicInteraction &right) const
 
virtual const std::pair
< G4double, G4double
GetFatalEnergyCheckLevels () const
 
virtual std::pair< G4double,
G4double
GetEnergyMomentumCheckLevels () const
 
void SetEnergyMomentumCheckLevels (G4double relativeLevel, G4double absoluteLevel)
 

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 49 of file G4BinaryLightIonReaction.cc.

References G4HadronicInteractionRegistry::FindModel(), G4VPreCompoundModel::GetExcitationHandler(), and G4HadronicInteractionRegistry::Instance().

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

Definition at line 68 of file G4BinaryLightIonReaction.cc.

69 {}

Member Function Documentation

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

Implements G4HadronicInteraction.

Definition at line 87 of file G4BinaryLightIonReaction.cc.

References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), debug, CLHEP::HepLorentzVector::e(), python.hepunit::eplus, G4cerr, G4cout, G4endl, G4lrint(), G4HadProjectile::Get4Momentum(), G4DynamicParticle::Get4Momentum(), G4Nucleus::GetA_asInt(), G4ParticleDefinition::GetBaryonNumber(), G4HadProjectile::GetDefinition(), G4IonTable::GetIonMass(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4HadFinalState::GetNumberOfSecondaries(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGCharge(), CLHEP::HepLorentzVector::getT(), G4HadProjectile::GetTotalEnergy(), G4Nucleus::GetZ_asInt(), CLHEP::HepLorentzRotation::inverse(), isAlive, python.hepunit::keV, CLHEP::Hep3Vector::mag(), python.hepunit::MeV, CLHEP::HepLorentzRotation::rotateY(), CLHEP::HepLorentzRotation::rotateZ(), G4DynamicParticle::Set4Momentum(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetMomentumChange(), G4HadFinalState::SetStatusChange(), CLHEP::HepLorentzVector::setT(), CLHEP::HepLorentzVector::setVect(), stopAndKill, CLHEP::Hep3Vector::unit(), and CLHEP::HepLorentzVector::vect().

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

Reimplemented from G4HadronicInteraction.

Definition at line 71 of file G4BinaryLightIonReaction.cc.

72 {
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";
78 }
std::ofstream outFile
Definition: GammaRayTel.cc:68
void G4BinaryLightIonReaction::SetDeExcitation ( G4ExcitationHandler ptr)
inline

Definition at line 74 of file G4BinaryLightIonReaction.hh.

References G4VPreCompoundModel::SetExcitationHandler().

75 {
76  theProjectileFragmentation->SetExcitationHandler(ptr);
77  theHandler = ptr;
78 }
void SetExcitationHandler(G4ExcitationHandler *ptr)
void G4BinaryLightIonReaction::SetPrecompound ( G4VPreCompoundModel ptr)
inline

Definition at line 69 of file G4BinaryLightIonReaction.hh.

References G4VPreCompoundModel::GetExcitationHandler().

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

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