Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4CollisionOutput Class Reference

#include <G4CollisionOutput.hh>

Public Member Functions

 G4CollisionOutput ()
 
G4CollisionOutputoperator= (const G4CollisionOutput &right)
 
void setVerboseLevel (G4int verbose)
 
void reset ()
 
void add (const G4CollisionOutput &right)
 
void addOutgoingParticle (const G4InuclElementaryParticle &particle)
 
void addOutgoingParticles (const std::vector< G4InuclElementaryParticle > &particles)
 
void addOutgoingNucleus (const G4InuclNuclei &nuclei)
 
void addOutgoingNuclei (const std::vector< G4InuclNuclei > &nuclea)
 
void addOutgoingParticle (const G4CascadParticle &cparticle)
 
void addOutgoingParticles (const std::vector< G4CascadParticle > &cparticles)
 
void addOutgoingParticles (const G4ReactionProductVector *rproducts)
 
void addRecoilFragment (const G4Fragment *aFragment)
 
void addRecoilFragment (const G4Fragment &aFragment)
 
void removeOutgoingParticle (G4int index)
 
void removeOutgoingParticle (const G4InuclElementaryParticle &particle)
 
void removeOutgoingParticle (const G4InuclElementaryParticle *particle)
 
void removeOutgoingNucleus (G4int index)
 
void removeOutgoingNucleus (const G4InuclNuclei &nuclei)
 
void removeOutgoingNucleus (const G4InuclNuclei *nuclei)
 
void removeRecoilFragment (G4int index=-1)
 
G4int numberOfOutgoingParticles () const
 
const std::vector
< G4InuclElementaryParticle > & 
getOutgoingParticles () const
 
std::vector
< G4InuclElementaryParticle > & 
getOutgoingParticles ()
 
G4int numberOfOutgoingNuclei () const
 
const std::vector
< G4InuclNuclei > & 
getOutgoingNuclei () const
 
std::vector< G4InuclNuclei > & getOutgoingNuclei ()
 
G4int numberOfFragments () const
 
const G4FragmentgetRecoilFragment (G4int index=0) const
 
const std::vector< G4Fragment > & getRecoilFragments () const
 
std::vector< G4Fragment > & getRecoilFragments ()
 
G4LorentzVector getTotalOutputMomentum () const
 
G4int getTotalCharge () const
 
G4int getTotalBaryonNumber () const
 
G4int getTotalStrangeness () const
 
void printCollisionOutput (std::ostream &os=G4cout) const
 
void boostToLabFrame (const G4LorentzConvertor &convertor)
 
G4LorentzVector boostToLabFrame (G4LorentzVector mom, const G4LorentzConvertor &convertor) const
 
void rotateEvent (const G4LorentzRotation &rotate)
 
void trivialise (G4InuclParticle *bullet, G4InuclParticle *target)
 
void setOnShell (G4InuclParticle *bullet, G4InuclParticle *target)
 
void setRemainingExitationEnergy ()
 
double getRemainingExitationEnergy () const
 
G4bool acceptable () const
 

Detailed Description

Definition at line 67 of file G4CollisionOutput.hh.

Constructor & Destructor Documentation

G4CollisionOutput::G4CollisionOutput ( )

Definition at line 82 of file G4CollisionOutput.cc.

83  : verboseLevel(0), eex_rest(0), on_shell(false) {
84  if (verboseLevel > 1)
85  G4cout << " >>> G4CollisionOutput::G4CollisionOutput" << G4endl;
86 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Member Function Documentation

G4bool G4CollisionOutput::acceptable ( ) const
inline

Definition at line 174 of file G4CollisionOutput.hh.

174 { return on_shell; };

Here is the caller graph for this function:

void G4CollisionOutput::add ( const G4CollisionOutput right)

Definition at line 121 of file G4CollisionOutput.cc.

121  {
122  addOutgoingParticles(right.outgoingParticles);
123  addOutgoingNuclei(right.outgoingNuclei);
124  recoilFragments = right.recoilFragments; // Replace, don't combine
125  eex_rest = 0.;
126  on_shell = false;
127 }
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::addOutgoingNuclei ( const std::vector< G4InuclNuclei > &  nuclea)

Definition at line 137 of file G4CollisionOutput.cc.

137  {
138  outgoingNuclei.insert(outgoingNuclei.end(), nuclea.begin(), nuclea.end());
139 }

Here is the caller graph for this function:

void G4CollisionOutput::addOutgoingNucleus ( const G4InuclNuclei nuclei)
inline

Definition at line 86 of file G4CollisionOutput.hh.

86  {
87  outgoingNuclei.push_back(nuclei);
88  };

Here is the caller graph for this function:

void G4CollisionOutput::addOutgoingParticle ( const G4InuclElementaryParticle particle)
inline

Definition at line 80 of file G4CollisionOutput.hh.

80  {
81  outgoingParticles.push_back(particle);
82  }

Here is the caller graph for this function:

void G4CollisionOutput::addOutgoingParticle ( const G4CascadParticle cparticle)

Definition at line 143 of file G4CollisionOutput.cc.

143  {
144  addOutgoingParticle(cparticle.getParticle());
145 }
void addOutgoingParticle(const G4InuclElementaryParticle &particle)
const G4InuclElementaryParticle & getParticle() const

Here is the call graph for this function:

void G4CollisionOutput::addOutgoingParticles ( const std::vector< G4InuclElementaryParticle > &  particles)

Definition at line 132 of file G4CollisionOutput.cc.

132  {
133  outgoingParticles.insert(outgoingParticles.end(),
134  particles.begin(), particles.end());
135 }

Here is the caller graph for this function:

void G4CollisionOutput::addOutgoingParticles ( const std::vector< G4CascadParticle > &  cparticles)

Definition at line 147 of file G4CollisionOutput.cc.

147  {
148  for (unsigned i=0; i<cparticles.size(); i++)
149  addOutgoingParticle(cparticles[i]);
150 }
void addOutgoingParticle(const G4InuclElementaryParticle &particle)

Here is the call graph for this function:

void G4CollisionOutput::addOutgoingParticles ( const G4ReactionProductVector rproducts)

Definition at line 154 of file G4CollisionOutput.cc.

154  {
155  if (!rproducts) return; // Sanity check, no error if null
156 
157  if (verboseLevel) {
158  G4cout << " >>> G4CollisionOutput::addOutgoingParticles(G4RPVector)"
159  << G4endl;
160  }
161 
162  G4ReactionProductVector::const_iterator j;
163  for (j=rproducts->begin(); j!=rproducts->end(); ++j) {
164  const G4ParticleDefinition* pd = (*j)->GetDefinition();
166 
167  // FIXME: Momentum returned by value; extra copying here!
168  G4LorentzVector mom((*j)->GetMomentum(), (*j)->GetTotalEnergy());
169  mom /= GeV; // Convert from GEANT4 to Bertini units
170 
171  if (verboseLevel>1)
172  G4cout << " Processing " << pd->GetParticleName() << " (" << type
173  << "), momentum " << mom << " GeV" << G4endl;
174 
175  // Nucleons and nuclei are jumbled together in the list
176  // NOTE: Resize and set/fill avoid unnecessary temporary copies
177  if (type) {
178  outgoingParticles.resize(numberOfOutgoingParticles()+1);
179  outgoingParticles.back().fill(mom, pd, G4InuclParticle::PreCompound);
180 
181  if (verboseLevel>1) G4cout << outgoingParticles.back() << G4endl;
182  } else {
183  outgoingNuclei.resize(numberOfOutgoingNuclei()+1);
184  outgoingNuclei.back().fill(mom,pd->GetAtomicMass(),pd->GetAtomicNumber(),
186 
187  if (verboseLevel>1) G4cout << outgoingNuclei.back() << G4endl;
188  }
189  }
190 }
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4int GetAtomicNumber() const
G4GLOB_DLL std::ostream G4cout
G4int numberOfOutgoingParticles() const
G4int GetAtomicMass() const
G4int numberOfOutgoingNuclei() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

void G4CollisionOutput::addRecoilFragment ( const G4Fragment aFragment)
inline

Definition at line 99 of file G4CollisionOutput.hh.

99  {
100  if (aFragment) addRecoilFragment(*aFragment);
101  }
void addRecoilFragment(const G4Fragment *aFragment)

Here is the caller graph for this function:

void G4CollisionOutput::addRecoilFragment ( const G4Fragment aFragment)
inline

Definition at line 103 of file G4CollisionOutput.hh.

103  {
104  recoilFragments.push_back(aFragment);
105  }
void G4CollisionOutput::boostToLabFrame ( const G4LorentzConvertor convertor)

Definition at line 319 of file G4CollisionOutput.cc.

319  {
320  if (verboseLevel > 1)
321  G4cout << " >>> G4CollisionOutput::boostToLabFrame" << G4endl;
322 
323  particleIterator ipart = outgoingParticles.begin();
324  for(; ipart != outgoingParticles.end(); ipart++) {
325  ipart->setMomentum(boostToLabFrame(ipart->getMomentum(), convertor));
326  }
327 
328  std::sort(outgoingParticles.begin(), outgoingParticles.end(),
330 
331  nucleiIterator inuc = outgoingNuclei.begin();
332  for (; inuc != outgoingNuclei.end(); inuc++) {
333  inuc->setMomentum(boostToLabFrame(inuc->getMomentum(), convertor));
334  }
335 
336  // NOTE: Fragment momentum must be converted to and from Bertini units
337  G4LorentzVector fmom;
338  fragmentIterator ifrag = recoilFragments.begin();
339  for (; ifrag != recoilFragments.end(); ifrag++) {
340  fmom = ifrag->GetMomentum() / GeV;
341  ifrag->SetMomentum(boostToLabFrame(fmom, convertor)*GeV);
342  }
343 }
G4GLOB_DLL std::ostream G4cout
void boostToLabFrame(const G4LorentzConvertor &convertor)
static constexpr double GeV
Definition: G4SIunits.hh:217
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
std::vector< G4Fragment >::iterator fragmentIterator
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

G4LorentzVector G4CollisionOutput::boostToLabFrame ( G4LorentzVector  mom,
const G4LorentzConvertor convertor 
) const

Definition at line 347 of file G4CollisionOutput.cc.

348  {
349  if (convertor.reflectionNeeded()) mom.setZ(-mom.z());
350  mom = convertor.rotate(mom);
351  mom = convertor.backToTheLab(mom);
352 
353  return mom;
354 }
G4LorentzVector rotate(const G4LorentzVector &mom) const
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4bool reflectionNeeded() const

Here is the call graph for this function:

const std::vector<G4InuclNuclei>& G4CollisionOutput::getOutgoingNuclei ( ) const
inline

Definition at line 137 of file G4CollisionOutput.hh.

137  {
138  return outgoingNuclei;
139  };

Here is the caller graph for this function:

std::vector<G4InuclNuclei>& G4CollisionOutput::getOutgoingNuclei ( )
inline

Definition at line 141 of file G4CollisionOutput.hh.

141 { return outgoingNuclei; };
const std::vector<G4InuclElementaryParticle>& G4CollisionOutput::getOutgoingParticles ( ) const
inline

Definition at line 127 of file G4CollisionOutput.hh.

127  {
128  return outgoingParticles;
129  };

Here is the caller graph for this function:

std::vector<G4InuclElementaryParticle>& G4CollisionOutput::getOutgoingParticles ( )
inline

Definition at line 131 of file G4CollisionOutput.hh.

131  {
132  return outgoingParticles;
133  };
const G4Fragment & G4CollisionOutput::getRecoilFragment ( G4int  index = 0) const

Definition at line 113 of file G4CollisionOutput.cc.

113  {
114  return ( (index >= 0 && index < numberOfFragments())
115  ? recoilFragments[index] : emptyFragment);
116 }
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

const std::vector<G4Fragment>& G4CollisionOutput::getRecoilFragments ( ) const
inline

Definition at line 147 of file G4CollisionOutput.hh.

147  {
148  return recoilFragments;
149  };
std::vector<G4Fragment>& G4CollisionOutput::getRecoilFragments ( )
inline

Definition at line 151 of file G4CollisionOutput.hh.

151 { return recoilFragments; };
double G4CollisionOutput::getRemainingExitationEnergy ( ) const
inline

Definition at line 173 of file G4CollisionOutput.hh.

173 { return eex_rest; };
G4int G4CollisionOutput::getTotalBaryonNumber ( ) const

Definition at line 268 of file G4CollisionOutput.cc.

268  {
269  if (verboseLevel > 1)
270  G4cout << " >>> G4CollisionOutput::getTotalBaryonNumber" << G4endl;
271 
272  G4int baryon = 0;
273  G4int i(0);
274  for(i=0; i < numberOfOutgoingParticles(); i++) {
275  baryon += outgoingParticles[i].baryon();
276  }
277  for(i=0; i < numberOfOutgoingNuclei(); i++) {
278  baryon += G4int(outgoingNuclei[i].getA());
279  }
280  for(i=0; i < numberOfFragments(); i++) {
281  baryon += recoilFragments[i].GetA_asInt();
282  }
283 
284  return baryon;
285 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4int numberOfOutgoingParticles() const
G4int numberOfOutgoingNuclei() const
#define G4endl
Definition: G4ios.hh:61
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4CollisionOutput::getTotalCharge ( ) const

Definition at line 249 of file G4CollisionOutput.cc.

249  {
250  if (verboseLevel > 1)
251  G4cout << " >>> G4CollisionOutput::getTotalCharge" << G4endl;
252 
253  G4int charge = 0;
254  G4int i(0);
255  for(i=0; i < numberOfOutgoingParticles(); i++) {
256  charge += G4int(outgoingParticles[i].getCharge());
257  }
258  for(i=0; i < numberOfOutgoingNuclei(); i++) {
259  charge += G4int(outgoingNuclei[i].getCharge());
260  }
261  for(i=0; i < numberOfFragments(); i++) {
262  charge += recoilFragments[i].GetZ_asInt();
263  }
264 
265  return charge;
266 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4int numberOfOutgoingParticles() const
G4int numberOfOutgoingNuclei() const
#define G4endl
Definition: G4ios.hh:61
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4LorentzVector G4CollisionOutput::getTotalOutputMomentum ( ) const

Definition at line 230 of file G4CollisionOutput.cc.

230  {
231  if (verboseLevel > 1)
232  G4cout << " >>> G4CollisionOutput::getTotalOutputMomentum" << G4endl;
233 
234  G4LorentzVector tot_mom;
235  G4int i(0);
236  for(i=0; i < numberOfOutgoingParticles(); i++) {
237  tot_mom += outgoingParticles[i].getMomentum();
238  }
239  for(i=0; i < numberOfOutgoingNuclei(); i++) {
240  tot_mom += outgoingNuclei[i].getMomentum();
241  }
242  for(i=0; i < numberOfFragments(); i++) {
243  tot_mom += recoilFragments[i].GetMomentum()/GeV; // Need Bertini units!
244  }
245 
246  return tot_mom;
247 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4int numberOfOutgoingParticles() const
G4int numberOfOutgoingNuclei() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4CollisionOutput::getTotalStrangeness ( ) const

Definition at line 287 of file G4CollisionOutput.cc.

287  {
288  if (verboseLevel > 1)
289  G4cout << " >>> G4CollisionOutput::getTotalStrangeness" << G4endl;
290 
291  G4int strange = 0;
292  G4int i(0);
293  for(i=0; i < numberOfOutgoingParticles(); i++) {
294  strange += outgoingParticles[i].getStrangeness();
295  }
296 
297  return strange;
298 }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
G4int numberOfOutgoingParticles() const
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4CollisionOutput::numberOfFragments ( ) const
inline

Definition at line 143 of file G4CollisionOutput.hh.

143 { return recoilFragments.size(); }

Here is the caller graph for this function:

G4int G4CollisionOutput::numberOfOutgoingNuclei ( ) const
inline

Definition at line 135 of file G4CollisionOutput.hh.

135 { return outgoingNuclei.size(); };

Here is the caller graph for this function:

G4int G4CollisionOutput::numberOfOutgoingParticles ( ) const
inline

Definition at line 125 of file G4CollisionOutput.hh.

125 { return outgoingParticles.size(); }

Here is the caller graph for this function:

G4CollisionOutput & G4CollisionOutput::operator= ( const G4CollisionOutput right)

Definition at line 89 of file G4CollisionOutput.cc.

90 {
91  if (this != &right) {
92  verboseLevel = right.verboseLevel;
93  outgoingParticles = right.outgoingParticles;
94  outgoingNuclei = right.outgoingNuclei;
95  recoilFragments = right.recoilFragments;
96  eex_rest = right.eex_rest;
97  on_shell = right.on_shell;
98  }
99  return *this;
100 }
void G4CollisionOutput::printCollisionOutput ( std::ostream &  os = G4cout) const

Definition at line 301 of file G4CollisionOutput.cc.

301  {
302  os << " Output: " << G4endl
303  << " Outgoing Particles: " << numberOfOutgoingParticles() << G4endl;
304 
305  G4int i(0);
306  for(i=0; i < numberOfOutgoingParticles(); i++)
307  os << outgoingParticles[i] << G4endl;
308 
309  os << " Outgoing Nuclei: " << numberOfOutgoingNuclei() << G4endl;
310  for(i=0; i < numberOfOutgoingNuclei(); i++)
311  os << outgoingNuclei[i] << G4endl;
312  for(i=0; i < numberOfFragments(); i++)
313  os << recoilFragments[i] << G4endl;
314 }
int G4int
Definition: G4Types.hh:78
G4int numberOfOutgoingParticles() const
G4int numberOfOutgoingNuclei() const
#define G4endl
Definition: G4ios.hh:61
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::removeOutgoingNucleus ( G4int  index)

Definition at line 200 of file G4CollisionOutput.cc.

200  {
201  if (index >= 0 && index < numberOfOutgoingNuclei())
202  outgoingNuclei.erase(outgoingNuclei.begin()+(size_t)index);
203 }
G4int numberOfOutgoingNuclei() const

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::removeOutgoingNucleus ( const G4InuclNuclei nuclei)

Definition at line 213 of file G4CollisionOutput.cc.

213  {
215  std::find(outgoingNuclei.begin(), outgoingNuclei.end(), nuclei);
216  if (pos != outgoingNuclei.end()) outgoingNuclei.erase(pos);
217 }
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
static const G4double pos
void G4CollisionOutput::removeOutgoingNucleus ( const G4InuclNuclei nuclei)
inline

Definition at line 117 of file G4CollisionOutput.hh.

117  {
118  if (nuclei) removeOutgoingNucleus(*nuclei);
119  }
void removeOutgoingNucleus(G4int index)

Here is the call graph for this function:

void G4CollisionOutput::removeOutgoingParticle ( G4int  index)

Definition at line 195 of file G4CollisionOutput.cc.

195  {
196  if (index >= 0 && index < numberOfOutgoingParticles())
197  outgoingParticles.erase(outgoingParticles.begin()+(size_t)index);
198 }
G4int numberOfOutgoingParticles() const

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::removeOutgoingParticle ( const G4InuclElementaryParticle particle)

Definition at line 207 of file G4CollisionOutput.cc.

207  {
209  std::find(outgoingParticles.begin(), outgoingParticles.end(), particle);
210  if (pos != outgoingParticles.end()) outgoingParticles.erase(pos);
211 }
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
static const G4double pos
void G4CollisionOutput::removeOutgoingParticle ( const G4InuclElementaryParticle particle)
inline

Definition at line 111 of file G4CollisionOutput.hh.

111  {
112  if (particle) removeOutgoingParticle(*particle);
113  }
void removeOutgoingParticle(G4int index)

Here is the call graph for this function:

void G4CollisionOutput::removeRecoilFragment ( G4int  index = -1)

Definition at line 221 of file G4CollisionOutput.cc.

221  {
222  if (index < 0) recoilFragments.clear();
223  else if (index < numberOfFragments())
224  recoilFragments.erase(recoilFragments.begin()+(size_t)index);
225 }
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::reset ( )

Definition at line 102 of file G4CollisionOutput.cc.

102  {
103  outgoingNuclei.clear();
104  outgoingParticles.clear();
105  recoilFragments.clear();
106  eex_rest = 0.;
107  on_shell = false;
108 }

Here is the caller graph for this function:

void G4CollisionOutput::rotateEvent ( const G4LorentzRotation rotate)

Definition at line 358 of file G4CollisionOutput.cc.

358  {
359  if (verboseLevel > 1)
360  G4cout << " >>> G4CollisionOutput::rotateEvent" << G4endl;
361 
362  particleIterator ipart = outgoingParticles.begin();
363  for(; ipart != outgoingParticles.end(); ipart++)
364  ipart->setMomentum(ipart->getMomentum()*=rotate);
365 
366  nucleiIterator inuc = outgoingNuclei.begin();
367  for (; inuc != outgoingNuclei.end(); inuc++)
368  inuc->setMomentum(inuc->getMomentum()*=rotate);
369 
370  fragmentIterator ifrag = recoilFragments.begin();
371  for (; ifrag != recoilFragments.end(); ifrag++) {
372  G4LorentzVector mom = ifrag->GetMomentum(); // Need copy
373  ifrag->SetMomentum(mom*=rotate);
374  }
375 }
G4GLOB_DLL std::ostream G4cout
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:65
std::vector< G4Fragment >::iterator fragmentIterator
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

void G4CollisionOutput::setOnShell ( G4InuclParticle bullet,
G4InuclParticle target 
)

Definition at line 403 of file G4CollisionOutput.cc.

404  {
405  if (verboseLevel > 1)
406  G4cout << " >>> G4CollisionOutput::setOnShell" << G4endl;
407 
408  const G4double accuracy = 0.00001; // momentum concerves at the level of 10 keV
409 
410  on_shell = false;
411 
412  G4LorentzVector ini_mom = bullet->getMomentum();
413  G4LorentzVector momt = target->getMomentum();
415  if(verboseLevel > 2){
416  G4cout << " bullet momentum = " << ini_mom.e() <<", "<< ini_mom.x() <<", "<< ini_mom.y()<<", "<< ini_mom.z()<<G4endl;
417  G4cout << " target momentum = " << momt.e()<<", "<< momt.x()<<", "<< momt.y()<<", "<< momt.z()<<G4endl;
418  G4cout << " Fstate momentum = " << out_mom.e()<<", "<< out_mom.x()<<", "<< out_mom.y()<<", "<< out_mom.z()<<G4endl;
419  }
420 
421  ini_mom += momt;
422 
423  mom_non_cons = ini_mom - out_mom;
424  G4double pnc = mom_non_cons.rho();
425  G4double enc = mom_non_cons.e();
426 
428 
429  if(verboseLevel > 2) {
431  G4cout << " momentum non conservation: " << G4endl
432  << " e " << enc << " p " << pnc << G4endl
433  << " remaining exitation " << eex_rest << G4endl;
434  }
435 
436  if (std::fabs(enc) <= accuracy && pnc <= accuracy) {
437  on_shell = true;
438  return;
439  }
440 
441  // Adjust "last" particle's four-momentum to balance event
442  // ONLY adjust particles with sufficient e or p to remain physical!
443 
444  if (verboseLevel > 2) G4cout << " re-balancing four-momenta" << G4endl;
445 
447  G4int nnuc = numberOfOutgoingNuclei();
448  G4int nfrag = numberOfFragments();
449 
450  G4LorentzVector last_mom; // Buffer to reduce memory churn
451 
452  if (npart > 0) {
453  for (G4int ip=npart-1; ip>=0; ip--) {
454  if (outgoingParticles[ip].getKineticEnergy()+enc > 0.) {
455  last_mom = outgoingParticles[ip].getMomentum();
456  last_mom += mom_non_cons;
457  outgoingParticles[ip].setMomentum(last_mom);
458  break;
459  }
460  }
461  } else if (nnuc > 0) {
462  for (G4int in=nnuc-1; in>=0; in--) {
463  if (outgoingNuclei[in].getKineticEnergy()+enc > 0.) {
464  last_mom = outgoingNuclei[in].getMomentum();
465  last_mom += mom_non_cons;
466  outgoingNuclei[in].setMomentum(last_mom);
467  break;
468  }
469  }
470  } else if (nfrag > 0) {
471  for (G4int ifr=nfrag-1; ifr>=0; ifr--) {
472  // NOTE: G4Fragment does not use Bertini units!
473  last_mom = recoilFragments[ifr].GetMomentum()/GeV;
474  if ((last_mom.e()-last_mom.m())+enc > 0.) {
475  last_mom += mom_non_cons;
476  recoilFragments[ifr].SetMomentum(last_mom*GeV);
477  break;
478  }
479  }
480  }
481 
482  // Recompute momentum non-conservation parameters
483  out_mom = getTotalOutputMomentum();
484  mom_non_cons = ini_mom - out_mom;
485  pnc = mom_non_cons.rho();
486  enc = mom_non_cons.e();
487 
488  if(verboseLevel > 2){
490  G4cout << " momentum non conservation after (1): " << G4endl
491  << " e " << enc << " p " << pnc << G4endl;
492  }
493 
494  // Can energy be balanced just with nuclear excitation?
495  G4bool need_hard_tuning = true;
496 
497  G4double encMeV = mom_non_cons.e() / GeV; // Excitation below is in MeV
498  if (nfrag > 0) {
499  for (G4int i=0; i < nfrag; i++) {
500  G4double eex = recoilFragments[0].GetExcitationEnergy();
501  if (eex > 0.0 && eex + encMeV >= 0.0) {
502  // NOTE: G4Fragment doesn't have function to change excitation energy
503  // ==> theRecoilFragment.SetExcitationEnergy(eex+encMeV);
504 
505  G4LorentzVector fragMom = recoilFragments[i].GetMomentum();
506  G4double newMass = fragMom.m() + encMeV; // .m() includes eex
507  fragMom.setVectM(fragMom.vect(), newMass);
508  need_hard_tuning = false;
509  break;
510  }
511  }
512  } else if (nnuc > 0) {
513  for (G4int i=0; i < nnuc; i++) {
514  G4double eex = outgoingNuclei[i].getExitationEnergy();
515  if (eex > 0.0 && eex + encMeV >= 0.0) {
516  outgoingNuclei[i].setExitationEnergy(eex+encMeV);
517  need_hard_tuning = false;
518  break;
519  }
520  }
521  if (need_hard_tuning && encMeV > 0.) {
522  outgoingNuclei[0].setExitationEnergy(encMeV);
523  need_hard_tuning = false;
524  }
525  }
526 
527  if (!need_hard_tuning) {
528  on_shell = true;
529  return;
530  }
531 
532  // Momentum (hard) tuning required for energy conservation
533  if (verboseLevel > 2)
534  G4cout << " trying hard (particle-pair) tuning" << G4endl;
535 
536  /*****
537  // Hard tuning of quasielastic particle against nucleus
538  if (npart == 1) {
539  if (verboseLevel > 2)
540  G4cout << " tuning particle 0 against recoil nucleus" << G4endl;
541 
542  G4LorentzVector mom1 = outgoingParticles[0].getMomentum();
543  G4LorentzVector mom2 = recoilFragments[0].GetMomentum()/GeV;
544 
545  // Preserve momentum direction of outgoing particle
546  G4ThreeVector pdir = mom1.vect().unit();
547 
548  // Two-body kinematics (nucleon against nucleus) in overall CM system
549  G4double ecmsq = ini_mom.m2();
550  G4double a = 0.5 * (ecmsq - mom1.m2() - mom2.m2());
551  G4double pmod = std::sqrt((a*a - mom1.m2()*mom2.m2()) / ecmsq);
552 
553  mom1.setVectM(pdir*pmod, mom1.m());
554  mom2.setVectM(-pdir*pmod, mom2.m());
555 
556  // Boost CM momenta back into lab frame, assign back to particles
557  mom1.boost(-ini_mom.boostVector());
558  mom2.boost(-ini_mom.boostVector());
559 
560  outgoingParticles[0].setMomentum(mom1);
561  recoilFragments[0].SetMomentum(mom2*GeV);
562  } else { // Hard tuning using pair of outgoing particles
563  *****/
564  std::pair<std::pair<G4int, G4int>, G4int> tune_par = selectPairToTune(enc);
565  std::pair<G4int, G4int> tune_particles = tune_par.first;
566  G4int mom_ind = tune_par.second;
567 
568  G4bool tuning_possible =
569  (tune_particles.first >= 0 && tune_particles.second >= 0 &&
570  mom_ind >= G4LorentzVector::X);
571 
572  if (!tuning_possible) {
573  if (verboseLevel > 2) G4cout << " tuning impossible " << G4endl;
574  return;
575  }
576 
577  if(verboseLevel > 2) {
578  G4cout << " p1 " << tune_particles.first << " p2 " << tune_particles.second
579  << " ind " << mom_ind << G4endl;
580  }
581 
582  G4LorentzVector mom1 = outgoingParticles[tune_particles.first].getMomentum();
583  G4LorentzVector mom2 = outgoingParticles[tune_particles.second].getMomentum();
584 
585  if (tuneSelectedPair(mom1, mom2, mom_ind)) {
586  outgoingParticles[tune_particles.first ].setMomentum(mom1);
587  outgoingParticles[tune_particles.second].setMomentum(mom2);
588  }
589  else return;
590  //***** }
591 
592  out_mom = getTotalOutputMomentum();
593  std::sort(outgoingParticles.begin(), outgoingParticles.end(), G4ParticleLargerEkin());
594  mom_non_cons = ini_mom - out_mom;
595  pnc = mom_non_cons.rho();
596  enc = mom_non_cons.e();
597 
598  on_shell = (std::fabs(enc) < accuracy || pnc < accuracy);
599 
600  if(verboseLevel > 2) {
601  G4cout << " momentum non conservation tuning: " << G4endl
602  << " e " << enc << " p " << pnc
603  << (on_shell?" success":" FAILURE") << G4endl;
604  }
605 }
G4LorentzVector getMomentum() const
void printCollisionOutput(std::ostream &os=G4cout) const
int G4int
Definition: G4Types.hh:78
G4LorentzVector getTotalOutputMomentum() const
void setVectM(const Hep3Vector &spatial, double mass)
Hep3Vector vect() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4int numberOfOutgoingParticles() const
double rho() const
G4int numberOfOutgoingNuclei() const
static constexpr double GeV
Definition: G4SIunits.hh:217
#define G4endl
Definition: G4ios.hh:61
G4int numberOfFragments() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::setRemainingExitationEnergy ( )

Definition at line 610 of file G4CollisionOutput.cc.

610  {
611  eex_rest = 0.;
612  G4int i(0);
613  for (i=0; i < numberOfOutgoingNuclei(); i++)
614  eex_rest += outgoingNuclei[i].getExitationEnergyInGeV();
615  for (i=0; i < numberOfFragments(); i++)
616  eex_rest += recoilFragments[i].GetExcitationEnergy() / GeV;
617 }
int G4int
Definition: G4Types.hh:78
G4int numberOfOutgoingNuclei() const
static constexpr double GeV
Definition: G4SIunits.hh:217
G4int numberOfFragments() const

Here is the call graph for this function:

Here is the caller graph for this function:

void G4CollisionOutput::setVerboseLevel ( G4int  verbose)
inline

Definition at line 72 of file G4CollisionOutput.hh.

72 { verboseLevel = verbose; };

Here is the caller graph for this function:

void G4CollisionOutput::trivialise ( G4InuclParticle bullet,
G4InuclParticle target 
)

Definition at line 378 of file G4CollisionOutput.cc.

379  {
380  if (verboseLevel > 1)
381  G4cout << " >>> G4CollisionOutput::trivialize" << G4endl;
382 
383  reset(); // Discard existing output, replace with bullet/target
384 
385  if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
386  outgoingNuclei.push_back(*nuclei_target);
387  } else {
388  G4InuclElementaryParticle* particle =
389  dynamic_cast<G4InuclElementaryParticle*>(target);
390  outgoingParticles.push_back(*particle);
391  }
392 
393  if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
394  outgoingNuclei.push_back(*nuclei_bullet);
395  } else {
396  G4InuclElementaryParticle* particle =
397  dynamic_cast<G4InuclElementaryParticle*>(bullet);
398  outgoingParticles.push_back(*particle);
399  }
400 }
const XML_Char * target
Definition: expat.h:268
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

Here is the caller graph for this function:


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