Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
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 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 ()
 
const G4FragmentgetRecoilFragment () const
 
G4FragmentgetRecoilFragment ()
 
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 65 of file G4CollisionOutput.hh.

Constructor & Destructor Documentation

◆ G4CollisionOutput()

G4CollisionOutput::G4CollisionOutput ( )

Definition at line 77 of file G4CollisionOutput.cc.

78 : verboseLevel(0), eex_rest(0), on_shell(false) {
79 if (verboseLevel > 1)
80 G4cout << " >>> G4CollisionOutput::G4CollisionOutput" << G4endl;
81}
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout

Member Function Documentation

◆ acceptable()

G4bool G4CollisionOutput::acceptable ( ) const
inline

Definition at line 166 of file G4CollisionOutput.hh.

166{ return on_shell; };

Referenced by G4InuclCollider::collide(), and G4IntraNucleiCascader::finishCascade().

◆ add()

void G4CollisionOutput::add ( const G4CollisionOutput right)

Definition at line 106 of file G4CollisionOutput.cc.

106 {
107 addOutgoingParticles(right.outgoingParticles);
108 addOutgoingNuclei(right.outgoingNuclei);
109 theRecoilFragment = right.theRecoilFragment;
110}
void addOutgoingNuclei(const std::vector< G4InuclNuclei > &nuclea)
void addOutgoingParticles(const std::vector< G4InuclElementaryParticle > &particles)

Referenced by G4CascadeDeexcitation::collide(), G4PreCompoundDeexcitation::collide(), G4InuclCollider::collide(), G4CascadeCheckBalance::collide(), G4InuclCollider::deexcite(), G4IntraNucleiCascader::finalize(), and G4InuclCollider::rescatter().

◆ addOutgoingNuclei()

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

Definition at line 120 of file G4CollisionOutput.cc.

120 {
121 outgoingNuclei.insert(outgoingNuclei.end(), nuclea.begin(), nuclea.end());
122}

Referenced by add(), G4CascadeCheckBalance::collide(), and G4Fissioner::collide().

◆ addOutgoingNucleus()

void G4CollisionOutput::addOutgoingNucleus ( const G4InuclNuclei nuclei)
inline

Definition at line 84 of file G4CollisionOutput.hh.

84 {
85 outgoingNuclei.push_back(nuclei);
86 };

Referenced by G4EquilibriumEvaporator::collide(), and G4NonEquilibriumEvaporator::collide().

◆ addOutgoingParticle() [1/2]

void G4CollisionOutput::addOutgoingParticle ( const G4CascadParticle cparticle)

Definition at line 126 of file G4CollisionOutput.cc.

126 {
127 addOutgoingParticle(cparticle.getParticle());
128}
const G4InuclElementaryParticle & getParticle() const
void addOutgoingParticle(const G4InuclElementaryParticle &particle)

◆ addOutgoingParticle() [2/2]

◆ addOutgoingParticles() [1/3]

void G4CollisionOutput::addOutgoingParticles ( const G4ReactionProductVector rproducts)

Definition at line 137 of file G4CollisionOutput.cc.

137 {
138 if (!rproducts) return; // Sanity check, no error if null
139
140 if (verboseLevel) {
141 G4cout << " >>> G4CollisionOutput::addOutgoingParticles(G4RPVector)"
142 << G4endl;
143 }
144
145 G4ReactionProductVector::const_iterator j;
146 for (j=rproducts->begin(); j!=rproducts->end(); ++j) {
147 G4ParticleDefinition* pd = (*j)->GetDefinition();
149
150 // FIXME: Momentum returned by value; extra copying here!
151 G4LorentzVector mom((*j)->GetMomentum(), (*j)->GetTotalEnergy());
152 mom /= GeV; // Convert from GEANT4 to Bertini units
153
154 if (verboseLevel>1)
155 G4cout << " Processing " << pd->GetParticleName() << " (" << type
156 << "), momentum " << mom << " GeV" << G4endl;
157
158 // Nucleons and nuclei are jumbled together in the list
159 // NOTE: Resize and set/fill avoid unnecessary temporary copies
160 if (type) {
161 outgoingParticles.resize(numberOfOutgoingParticles()+1);
162 outgoingParticles.back().fill(mom, pd, G4InuclParticle::PreCompound);
163
164 if (verboseLevel>1) G4cout << outgoingParticles.back() << G4endl;
165 } else {
166 outgoingNuclei.resize(numberOfOutgoingNuclei()+1);
167 outgoingNuclei.back().fill(mom,pd->GetAtomicMass(),pd->GetAtomicNumber(),
169
170 if (verboseLevel>1) G4cout << outgoingNuclei.back() << G4endl;
171 }
172 }
173}
int G4int
Definition: G4Types.hh:66
G4int numberOfOutgoingParticles() const
G4int numberOfOutgoingNuclei() const
G4int GetAtomicNumber() const
G4int GetAtomicMass() const
const G4String & GetParticleName() const

◆ addOutgoingParticles() [2/3]

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

Definition at line 130 of file G4CollisionOutput.cc.

130 {
131 for (unsigned i=0; i<cparticles.size(); i++)
132 addOutgoingParticle(cparticles[i]);
133}

◆ addOutgoingParticles() [3/3]

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

◆ addRecoilFragment() [1/2]

void G4CollisionOutput::addRecoilFragment ( const G4Fragment aFragment)
inline

Definition at line 101 of file G4CollisionOutput.hh.

101 {
102 theRecoilFragment = aFragment;
103 }

◆ addRecoilFragment() [2/2]

void G4CollisionOutput::addRecoilFragment ( const G4Fragment aFragment)
inline

Definition at line 97 of file G4CollisionOutput.hh.

97 {
98 if (aFragment) addRecoilFragment(*aFragment);
99 }
void addRecoilFragment(const G4Fragment *aFragment)

Referenced by addRecoilFragment(), and G4IntraNucleiCascader::finishCascade().

◆ boostToLabFrame() [1/2]

void G4CollisionOutput::boostToLabFrame ( const G4LorentzConvertor convertor)

Definition at line 294 of file G4CollisionOutput.cc.

294 {
295 if (verboseLevel > 1)
296 G4cout << " >>> G4CollisionOutput::boostToLabFrame" << G4endl;
297
298 if (!outgoingParticles.empty()) {
299 particleIterator ipart = outgoingParticles.begin();
300 for(; ipart != outgoingParticles.end(); ipart++) {
301 ipart->setMomentum(boostToLabFrame(ipart->getMomentum(), convertor));
302 }
303
304 std::sort(outgoingParticles.begin(), outgoingParticles.end(),
306 }
307
308 if (!outgoingNuclei.empty()) {
309 nucleiIterator inuc = outgoingNuclei.begin();
310 for (; inuc != outgoingNuclei.end(); inuc++) {
311 inuc->setMomentum(boostToLabFrame(inuc->getMomentum(), convertor));
312 }
313 }
314
315 if (theRecoilFragment.GetA() > 0.) {
316 theRecoilFragment.SetMomentum(boostToLabFrame(theRecoilFragment.GetMomentum(), convertor));
317 }
318}
std::vector< G4InuclElementaryParticle >::iterator particleIterator
Definition: G4BigBanger.cc:60
std::vector< G4InuclNuclei >::const_iterator nucleiIterator
void boostToLabFrame(const G4LorentzConvertor &convertor)
const G4LorentzVector & GetMomentum() const
Definition: G4Fragment.hh:251
G4double GetA() const
Definition: G4Fragment.hh:283
void SetMomentum(const G4LorentzVector &value)
Definition: G4Fragment.hh:256

Referenced by boostToLabFrame(), G4InuclCollider::collide(), and G4EquilibriumEvaporator::collide().

◆ boostToLabFrame() [2/2]

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

Definition at line 322 of file G4CollisionOutput.cc.

323 {
324 if (convertor.reflectionNeeded()) mom.setZ(-mom.z());
325 mom = convertor.rotate(mom);
326 mom = convertor.backToTheLab(mom);
327
328 return mom;
329}
G4bool reflectionNeeded() const
G4LorentzVector backToTheLab(const G4LorentzVector &mom) const
G4LorentzVector rotate(const G4LorentzVector &mom) const

◆ getOutgoingNuclei() [1/2]

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

Definition at line 139 of file G4CollisionOutput.hh.

139{ return outgoingNuclei; };

◆ getOutgoingNuclei() [2/2]

◆ getOutgoingParticles() [1/2]

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

Definition at line 129 of file G4CollisionOutput.hh.

129 {
130 return outgoingParticles;
131 };

◆ getOutgoingParticles() [2/2]

◆ getRecoilFragment() [1/2]

G4Fragment & G4CollisionOutput::getRecoilFragment ( )
inline

Definition at line 143 of file G4CollisionOutput.hh.

143{ return theRecoilFragment; }

◆ getRecoilFragment() [2/2]

const G4Fragment & G4CollisionOutput::getRecoilFragment ( ) const
inline

Definition at line 141 of file G4CollisionOutput.hh.

141{ return theRecoilFragment; }

Referenced by G4InuclCollider::collide(), and G4InuclCollider::rescatter().

◆ getRemainingExitationEnergy()

double G4CollisionOutput::getRemainingExitationEnergy ( ) const
inline

Definition at line 165 of file G4CollisionOutput.hh.

165{ return eex_rest; };

◆ getTotalBaryonNumber()

G4int G4CollisionOutput::getTotalBaryonNumber ( ) const

Definition at line 245 of file G4CollisionOutput.cc.

245 {
246 if (verboseLevel > 1)
247 G4cout << " >>> G4CollisionOutput::getTotalBaryonNumber" << G4endl;
248
249 G4int baryon = 0;
250 G4int i(0);
251 for(i=0; i < G4int(outgoingParticles.size()); i++) {
252 baryon += outgoingParticles[i].baryon();
253 }
254 for(i=0; i < G4int(outgoingNuclei.size()); i++) {
255 baryon += G4int(outgoingNuclei[i].getA());
256 }
257 baryon += theRecoilFragment.GetA_asInt();
258
259 return baryon;
260}
G4int GetA_asInt() const
Definition: G4Fragment.hh:218

Referenced by G4CascadeCheckBalance::collide().

◆ getTotalCharge()

G4int G4CollisionOutput::getTotalCharge ( ) const

Definition at line 228 of file G4CollisionOutput.cc.

228 {
229 if (verboseLevel > 1)
230 G4cout << " >>> G4CollisionOutput::getTotalCharge" << G4endl;
231
232 G4int charge = 0;
233 G4int i(0);
234 for(i=0; i < G4int(outgoingParticles.size()); i++) {
235 charge += G4int(outgoingParticles[i].getCharge());
236 }
237 for(i=0; i < G4int(outgoingNuclei.size()); i++) {
238 charge += G4int(outgoingNuclei[i].getCharge());
239 }
240 charge += theRecoilFragment.GetZ_asInt();
241
242 return charge;
243}
G4int GetZ_asInt() const
Definition: G4Fragment.hh:223

Referenced by G4CascadeCheckBalance::collide().

◆ getTotalOutputMomentum()

G4LorentzVector G4CollisionOutput::getTotalOutputMomentum ( ) const

Definition at line 211 of file G4CollisionOutput.cc.

211 {
212 if (verboseLevel > 1)
213 G4cout << " >>> G4CollisionOutput::getTotalOutputMomentum" << G4endl;
214
215 G4LorentzVector tot_mom;
216 G4int i(0);
217 for(i=0; i < G4int(outgoingParticles.size()); i++) {
218 tot_mom += outgoingParticles[i].getMomentum();
219 }
220 for(i=0; i < G4int(outgoingNuclei.size()); i++) {
221 tot_mom += outgoingNuclei[i].getMomentum();
222 }
223 tot_mom += theRecoilFragment.GetMomentum()/GeV; // Need Bertini units!
224
225 return tot_mom;
226}

Referenced by G4CascadeCheckBalance::collide(), and setOnShell().

◆ getTotalStrangeness()

G4int G4CollisionOutput::getTotalStrangeness ( ) const

Definition at line 262 of file G4CollisionOutput.cc.

262 {
263 if (verboseLevel > 1)
264 G4cout << " >>> G4CollisionOutput::getTotalStrangeness" << G4endl;
265
266 G4int strange = 0;
267 G4int i(0);
268 for(i=0; i < G4int(outgoingParticles.size()); i++) {
269 strange += outgoingParticles[i].getStrangeness();
270 }
271
272 return strange;
273}

Referenced by G4CascadeCheckBalance::collide().

◆ numberOfOutgoingNuclei()

G4int G4CollisionOutput::numberOfOutgoingNuclei ( ) const
inline

◆ numberOfOutgoingParticles()

◆ operator=()

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

Definition at line 84 of file G4CollisionOutput.cc.

85{
86 if (this != &right) {
87 verboseLevel = right.verboseLevel;
88 outgoingParticles = right.outgoingParticles;
89 outgoingNuclei = right.outgoingNuclei;
90 theRecoilFragment = right.theRecoilFragment;
91 eex_rest = right.eex_rest;
92 on_shell = right.on_shell;
93 }
94 return *this;
95}

◆ printCollisionOutput()

void G4CollisionOutput::printCollisionOutput ( std::ostream &  os = G4cout) const

Definition at line 276 of file G4CollisionOutput.cc.

276 {
277 os << " Output: " << G4endl
278 << " Outgoing Particles: " << outgoingParticles.size() << G4endl;
279
280 G4int i(0);
281 for(i=0; i < G4int(outgoingParticles.size()); i++)
282 os << outgoingParticles[i] << G4endl;
283
284 os << " Outgoing Nuclei: " << outgoingNuclei.size() << G4endl;
285 for(i=0; i < G4int(outgoingNuclei.size()); i++)
286 os << outgoingNuclei[i] << G4endl;
287
288 if (theRecoilFragment.GetA() > 0) os << theRecoilFragment << G4endl;
289}

Referenced by G4CascadeInterface::ApplyYourself(), G4EvaporationInuclCollider::collide(), G4CascadeCoalescence::FindClusters(), G4IntraNucleiCascader::finishCascade(), G4NucleiModel::generateParticleFate(), G4CascadeInterface::Propagate(), setOnShell(), G4CascadeInterface::throwNonConservationFailure(), and G4CascadeColliderBase::validateOutput().

◆ removeOutgoingNucleus() [1/3]

void G4CollisionOutput::removeOutgoingNucleus ( const G4InuclNuclei nuclei)

Definition at line 196 of file G4CollisionOutput.cc.

196 {
197 nucleiIterator pos =
198 std::find(outgoingNuclei.begin(), outgoingNuclei.end(), nuclei);
199 if (pos != outgoingNuclei.end()) outgoingNuclei.erase(pos);
200}

◆ removeOutgoingNucleus() [2/3]

void G4CollisionOutput::removeOutgoingNucleus ( const G4InuclNuclei nuclei)
inline

Definition at line 115 of file G4CollisionOutput.hh.

115 {
116 if (nuclei) removeOutgoingNucleus(*nuclei);
117 }
void removeOutgoingNucleus(G4int index)

◆ removeOutgoingNucleus() [3/3]

void G4CollisionOutput::removeOutgoingNucleus ( G4int  index)

Definition at line 183 of file G4CollisionOutput.cc.

183 {
184 if (index >= 0 && index < numberOfOutgoingNuclei())
185 outgoingNuclei.erase(outgoingNuclei.begin()+(size_t)index);
186}

Referenced by removeOutgoingNucleus().

◆ removeOutgoingParticle() [1/3]

void G4CollisionOutput::removeOutgoingParticle ( const G4InuclElementaryParticle particle)

Definition at line 190 of file G4CollisionOutput.cc.

190 {
191 particleIterator pos =
192 std::find(outgoingParticles.begin(), outgoingParticles.end(), particle);
193 if (pos != outgoingParticles.end()) outgoingParticles.erase(pos);
194}

◆ removeOutgoingParticle() [2/3]

void G4CollisionOutput::removeOutgoingParticle ( const G4InuclElementaryParticle particle)
inline

Definition at line 109 of file G4CollisionOutput.hh.

109 {
110 if (particle) removeOutgoingParticle(*particle);
111 }
void removeOutgoingParticle(G4int index)

◆ removeOutgoingParticle() [3/3]

void G4CollisionOutput::removeOutgoingParticle ( G4int  index)

Definition at line 178 of file G4CollisionOutput.cc.

178 {
179 if (index >= 0 && index < numberOfOutgoingParticles())
180 outgoingParticles.erase(outgoingParticles.begin()+(size_t)index);
181}

Referenced by removeOutgoingParticle().

◆ removeRecoilFragment()

void G4CollisionOutput::removeRecoilFragment ( )

Definition at line 204 of file G4CollisionOutput.cc.

204 {
205 static const G4Fragment emptyFragment; // Default ctor is all zeros
206 theRecoilFragment = emptyFragment;
207}

Referenced by G4InuclCollider::collide(), G4InuclCollider::rescatter(), and reset().

◆ reset()

◆ rotateEvent()

void G4CollisionOutput::rotateEvent ( const G4LorentzRotation rotate)

Definition at line 333 of file G4CollisionOutput.cc.

333 {
334 if (verboseLevel > 1)
335 G4cout << " >>> G4CollisionOutput::rotateEvent" << G4endl;
336
337 particleIterator ipart = outgoingParticles.begin();
338 for(; ipart != outgoingParticles.end(); ipart++)
339 ipart->setMomentum(ipart->getMomentum()*=rotate);
340
341 nucleiIterator inuc = outgoingNuclei.begin();
342 for (; inuc != outgoingNuclei.end(); inuc++)
343 inuc->setMomentum(inuc->getMomentum()*=rotate);
344
345 if (theRecoilFragment.GetA() > 0.) {
346 G4LorentzVector mom = theRecoilFragment.GetMomentum(); // Need copy
347 theRecoilFragment.SetMomentum(mom*=rotate);
348 }
349
350}

Referenced by G4CascadeInterface::ApplyYourself().

◆ setOnShell()

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

Definition at line 378 of file G4CollisionOutput.cc.

379 {
380 if (verboseLevel > 1)
381 G4cout << " >>> G4CollisionOutput::setOnShell" << G4endl;
382
383 const G4double accuracy = 0.00001; // momentum concerves at the level of 10 keV
384
385 on_shell = false;
386
387 G4LorentzVector ini_mom = bullet->getMomentum();
388 G4LorentzVector momt = target->getMomentum();
390 if(verboseLevel > 2){
391 G4cout << " bullet momentum = " << ini_mom.e() <<", "<< ini_mom.x() <<", "<< ini_mom.y()<<", "<< ini_mom.z()<<G4endl;
392 G4cout << " target momentum = " << momt.e()<<", "<< momt.x()<<", "<< momt.y()<<", "<< momt.z()<<G4endl;
393 G4cout << " Fstate momentum = " << out_mom.e()<<", "<< out_mom.x()<<", "<< out_mom.y()<<", "<< out_mom.z()<<G4endl;
394 }
395
396 ini_mom += momt;
397 G4LorentzVector mom_non_cons = ini_mom - out_mom;
398
399 G4double pnc = mom_non_cons.rho();
400 G4double enc = mom_non_cons.e();
401
403
404 if(verboseLevel > 2) {
406 G4cout << " momentum non conservation: " << G4endl
407 << " e " << enc << " p " << pnc << G4endl
408 << " remaining exitation " << eex_rest << G4endl;
409 }
410
411 if (std::fabs(enc) <= accuracy && pnc <= accuracy) {
412 on_shell = true;
413 return;
414 }
415
416 // Adjust "last" particle's four-momentum to balance event
417 // ONLY adjust particles with sufficient e or p to remain physical!
418
419 if (verboseLevel > 2) G4cout << " re-balancing four-momenta" << G4endl;
420
421 G4int npart = outgoingParticles.size();
422 G4int nnuc = outgoingNuclei.size();
423
424 if (npart > 0) {
425 for (G4int ip=npart-1; ip>=0; ip--) {
426 if (outgoingParticles[ip].getKineticEnergy()+enc > 0.) {
427 G4LorentzVector last_mom = outgoingParticles[ip].getMomentum();
428 last_mom += mom_non_cons;
429 outgoingParticles[ip].setMomentum(last_mom);
430 break;
431 }
432 }
433 } else if (nnuc > 0) {
434 for (G4int in=nnuc-1; in>=0; in--) {
435 if (outgoingNuclei[in].getKineticEnergy()+enc > 0.) {
436 G4LorentzVector last_mom = outgoingNuclei[in].getMomentum();
437 last_mom += mom_non_cons;
438 outgoingNuclei[in].setMomentum(last_mom);
439 break;
440 }
441 }
442 } else if (theRecoilFragment.GetA() > 0.) {
443 // NOTE: G4Fragment does not use Bertini units!
444 G4LorentzVector last_mom = theRecoilFragment.GetMomentum()/GeV;
445 if ((last_mom.e()-last_mom.m())+enc > 0.) {
446 last_mom += mom_non_cons;
447 theRecoilFragment.SetMomentum(last_mom*GeV);
448 }
449 }
450
451 // Recompute momentum non-conservation parameters
452 out_mom = getTotalOutputMomentum();
453 mom_non_cons = ini_mom - out_mom;
454 pnc = mom_non_cons.rho();
455 enc = mom_non_cons.e();
456
457 if(verboseLevel > 2){
459 G4cout << " momentum non conservation after (1): " << G4endl
460 << " e " << enc << " p " << pnc << G4endl;
461 }
462
463 // Can energy be balanced just with nuclear excitation?
464 G4bool need_hard_tuning = true;
465
466 G4double encMeV = mom_non_cons.e() / GeV; // Excitation below is in MeV
467 if (theRecoilFragment.GetA() > 0.) {
468 G4double eex = theRecoilFragment.GetExcitationEnergy();
469 if (eex > 0.0 && eex + encMeV >= 0.0) {
470 // NOTE: G4Fragment doesn't have function to change excitation energy
471 // ==> theRecoilFragment.SetExcitationEnergy(eex+encMeV);
472
473 G4LorentzVector fragMom = theRecoilFragment.GetMomentum();
474 G4double newMass = fragMom.m() + encMeV; // .m() includes eex
475 fragMom.setVectM(fragMom.vect(), newMass);
476 need_hard_tuning = false;
477 }
478 } else if (outgoingNuclei.size() > 0) {
479 for (G4int i=0; i < G4int(outgoingNuclei.size()); i++) {
480 G4double eex = outgoingNuclei[i].getExitationEnergy();
481
482 if(eex > 0.0 && eex + encMeV >= 0.0) {
483 outgoingNuclei[i].setExitationEnergy(eex+encMeV);
484 need_hard_tuning = false;
485 break;
486 }
487 }
488 if (need_hard_tuning && encMeV > 0.) {
489 outgoingNuclei[0].setExitationEnergy(encMeV);
490 need_hard_tuning = false;
491 }
492 }
493
494 if (!need_hard_tuning) {
495 on_shell = true;
496 return;
497 }
498
499 // Momentum (hard) tuning required for energy conservation
500 if (verboseLevel > 2)
501 G4cout << " trying hard (particle-pair) tuning" << G4endl;
502
503 std::pair<std::pair<G4int, G4int>, G4int> tune_par = selectPairToTune(mom_non_cons.e());
504 std::pair<G4int, G4int> tune_particles = tune_par.first;
505 G4int mom_ind = tune_par.second;
506
507 if(verboseLevel > 2) {
508 G4cout << " p1 " << tune_particles.first << " p2 " << tune_particles.second
509 << " ind " << mom_ind << G4endl;
510 }
511
512 G4bool tuning_possible =
513 (tune_particles.first >= 0 && tune_particles.second >= 0 &&
514 mom_ind >= G4LorentzVector::X);
515
516 if (!tuning_possible) {
517 if (verboseLevel > 2) G4cout << " tuning impossible " << G4endl;
518 return;
519 }
520
521 G4LorentzVector mom1 = outgoingParticles[tune_particles.first].getMomentum();
522 G4LorentzVector mom2 = outgoingParticles[tune_particles.second].getMomentum();
523 G4double newE12 = mom1.e() + mom2.e() + mom_non_cons.e();
524 G4double R = 0.5 * (newE12 * newE12 + mom2.e() * mom2.e() - mom1.e() * mom1.e()) / newE12;
525 G4double Q = -(mom1[mom_ind] + mom2[mom_ind]) / newE12;
526 G4double UDQ = 1.0 / (Q * Q - 1.0);
527 G4double W = (R * Q + mom2[mom_ind]) * UDQ;
528 G4double V = (mom2.e() * mom2.e() - R * R) * UDQ;
529 G4double DET = W * W + V;
530
531 if (DET < 0.0) {
532 if (verboseLevel > 2) G4cout << " DET < 0 " << G4endl;
533 return;
534 }
535
536 // Tuning allowed only for non-negative determinant
537 G4double x1 = -(W + std::sqrt(DET));
538 G4double x2 = -(W - std::sqrt(DET));
539
540 // choose the appropriate solution
541 G4bool xset = false;
542 G4double x = 0.0;
543
544 if(mom_non_cons.e() > 0.0) { // x has to be > 0.0
545 if(x1 > 0.0) {
546 if(R + Q * x1 >= 0.0) {
547 x = x1;
548 xset = true;
549 };
550 };
551 if(!xset && x2 > 0.0) {
552 if(R + Q * x2 >= 0.0) {
553 x = x2;
554 xset = true;
555 };
556 };
557 } else {
558 if(x1 < 0.0) {
559 if(R + Q * x1 >= 0.) {
560 x = x1;
561 xset = true;
562 };
563 };
564 if(!xset && x2 < 0.0) {
565 if(R + Q * x2 >= 0.0) {
566 x = x2;
567 xset = true;
568 };
569 };
570 } // if(mom_non_cons.e() > 0.0)
571
572 if(!xset) {
573 if(verboseLevel > 2)
574 G4cout << " no appropriate solution found " << G4endl;
575 return;
576 } // if(xset)
577
578 // retune momentums
579 mom1[mom_ind] += x;
580 mom2[mom_ind] -= x;
581 outgoingParticles[tune_particles.first ].setMomentum(mom1);
582 outgoingParticles[tune_particles.second].setMomentum(mom2);
583 out_mom = getTotalOutputMomentum();
584 std::sort(outgoingParticles.begin(), outgoingParticles.end(), G4ParticleLargerEkin());
585 mom_non_cons = ini_mom - out_mom;
586 pnc = mom_non_cons.rho();
587 enc = mom_non_cons.e();
588
589 on_shell = (std::fabs(enc) < accuracy || pnc < accuracy);
590
591 if(verboseLevel > 2) {
592 G4cout << " momentum non conservation tuning: " << G4endl
593 << " e " << enc << " p " << pnc
594 << (on_shell?" success":" FAILURE") << G4endl;
595 }
596}
double G4double
Definition: G4Types.hh:64
bool G4bool
Definition: G4Types.hh:67
void setVectM(const Hep3Vector &spatial, double mass)
Hep3Vector vect() const
G4LorentzVector getTotalOutputMomentum() const
void printCollisionOutput(std::ostream &os=G4cout) const
G4double GetExcitationEnergy() const
Definition: G4Fragment.hh:235
G4LorentzVector getMomentum() const

Referenced by G4InuclCollider::collide(), and G4IntraNucleiCascader::finishCascade().

◆ setRemainingExitationEnergy()

void G4CollisionOutput::setRemainingExitationEnergy ( )

Definition at line 601 of file G4CollisionOutput.cc.

601 {
602 eex_rest = theRecoilFragment.GetExcitationEnergy() / GeV;
603
604 for(G4int i=0; i < G4int(outgoingNuclei.size()); i++)
605 eex_rest += outgoingNuclei[i].getExitationEnergyInGeV();
606}

Referenced by setOnShell().

◆ setVerboseLevel()

void G4CollisionOutput::setVerboseLevel ( G4int  verbose)
inline

◆ trivialise()

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

Definition at line 353 of file G4CollisionOutput.cc.

354 {
355 if (verboseLevel > 1)
356 G4cout << " >>> G4CollisionOutput::trivialize" << G4endl;
357
358 reset(); // Discard existing output, replace with bullet/target
359
360 if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
361 outgoingNuclei.push_back(*nuclei_target);
362 } else {
363 G4InuclElementaryParticle* particle =
364 dynamic_cast<G4InuclElementaryParticle*>(target);
365 outgoingParticles.push_back(*particle);
366 }
367
368 if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
369 outgoingNuclei.push_back(*nuclei_bullet);
370 } else {
371 G4InuclElementaryParticle* particle =
372 dynamic_cast<G4InuclElementaryParticle*>(bullet);
373 outgoingParticles.push_back(*particle);
374 }
375}

Referenced by G4InuclCollider::collide(), and G4IntraNucleiCascader::finalize().


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