ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNucleus.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLNucleus.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 /*
39  * G4INCLNucleus.cc
40  *
41  * \date Jun 5, 2009
42  * \author Pekka Kaitaniemi
43  */
44 
45 #ifndef G4INCLNucleus_hh
46 #define G4INCLNucleus_hh 1
47 
48 #include "G4INCLGlobals.hh"
49 #include "G4INCLLogger.hh"
50 #include "G4INCLParticle.hh"
51 #include "G4INCLIAvatar.hh"
52 #include "G4INCLNucleus.hh"
53 #include "G4INCLKinematicsUtils.hh"
54 #include "G4INCLDecayAvatar.hh"
56 #include "G4INCLCluster.hh"
57 #include "G4INCLClusterDecay.hh"
58 #include "G4INCLDeJongSpin.hh"
59 #include "G4INCLParticleSpecies.hh"
60 #include "G4INCLParticleTable.hh"
61 #include <iterator>
62 #include <cstdlib>
63 #include <sstream>
64 // #include <cassert>
65 
66 namespace G4INCL {
67 
68  Nucleus::Nucleus(G4int mass, G4int charge, G4int strangess, Config const * const conf, const G4double universeRadius)
69  : Cluster(charge,mass,strangess,true),
70  theInitialZ(charge), theInitialA(mass), theInitialS(strangess),
71  theNpInitial(0), theNnInitial(0),
72  theNpionplusInitial(0), theNpionminusInitial(0),
73  theNkaonplusInitial(0), theNkaonminusInitial(0),
74  initialInternalEnergy(0.),
75  incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
76  initialCenterOfMass(0.,0.,0.),
77  remnant(true),
78  initialEnergy(0.),
79  tryCN(false),
80  theUniverseRadius(universeRadius),
81  isNucleusNucleus(false),
82  theProjectileRemnant(NULL),
83  theDensity(NULL),
84  thePotential(NULL)
85  {
86  PotentialType potentialType;
87  G4bool pionPotential;
88  if(conf) {
89  potentialType = conf->getPotentialType();
90  pionPotential = conf->getPionPotential();
91  } else { // By default we don't use energy dependent
92  // potential. This is convenient for some tests.
93  potentialType = IsospinPotential;
94  pionPotential = true;
95  }
96 
97  thePotential = NuclearPotential::createPotential(potentialType, theA, theZ, pionPotential);
98 
101 
103 
105  theParticleSampler->setDensity(theDensity);
106 
107  if(theUniverseRadius<0)
108  theUniverseRadius = theDensity->getMaximumRadius();
109  theStore = new Store(conf);
110  }
111 
113  delete theStore;
115  /* We don't delete the potential and the density here any more -- Factories
116  * are caching them
117  delete thePotential;
118  delete theDensity;*/
119  }
120 
122  // Reset the variables connected with the projectile remnant
123  delete theProjectileRemnant;
124  theProjectileRemnant = NULL;
125 
127  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
129  }
131  particles.clear();
134  }
135 
137  if(!finalstate) // do nothing if no final state was returned
138  return;
139 
140  G4double totalEnergy = 0.0;
141 
142  FinalStateValidity const validity = finalstate->getValidity();
143  if(validity == ValidFS) {
144 
145  ParticleList const &created = finalstate->getCreatedParticles();
146  for(ParticleIter iter=created.begin(), e=created.end(); iter!=e; ++iter) {
147  theStore->add((*iter));
148  if(!(*iter)->isOutOfWell()) {
149  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
150  }
151  }
152 
153  ParticleList const &deleted = finalstate->getDestroyedParticles();
154  for(ParticleIter iter=deleted.begin(), e=deleted.end(); iter!=e; ++iter) {
156  }
157 
158  ParticleList const &modified = finalstate->getModifiedParticles();
159  for(ParticleIter iter=modified.begin(), e=modified.end(); iter!=e; ++iter) {
161  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
162  }
163 
164  ParticleList const &out = finalstate->getOutgoingParticles();
165  for(ParticleIter iter=out.begin(), e=out.end(); iter!=e; ++iter) {
166  if((*iter)->isCluster()) {
167  Cluster *clusterOut = dynamic_cast<Cluster*>((*iter));
168 // assert(clusterOut);
169 #ifdef INCLXX_IN_GEANT4_MODE
170  if(!clusterOut)
171  continue;
172 #endif
173  ParticleList const &components = clusterOut->getParticles();
174  for(ParticleIter in=components.begin(), end=components.end(); in!=end; ++in)
176  } else {
178  }
179  totalEnergy += (*iter)->getEnergy(); // No potential here because the particle is gone
180  theA -= (*iter)->getA();
181  theZ -= (*iter)->getZ();
182  theS -= (*iter)->getS();
183  theStore->addToOutgoing(*iter);
184  (*iter)->setEmissionTime(theStore->getBook().getCurrentTime());
185  }
186 
187  ParticleList const &entering = finalstate->getEnteringParticles();
188  for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
189  insertParticle(*iter);
190  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
191  }
192 
193  // actually perform the removal of the scheduled avatars
195  } else if(validity == ParticleBelowFermiFS || validity == ParticleBelowZeroFS) {
196  INCL_DEBUG("A Particle is entering below the Fermi sea:" << '\n' << finalstate->print() << '\n');
197  tryCN = true;
198  ParticleList const &entering = finalstate->getEnteringParticles();
199  for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
200  insertParticle(*iter);
201  }
202  }
203 
204  if(validity==ValidFS &&
205  std::abs(totalEnergy - finalstate->getTotalEnergyBeforeInteraction()) > 0.1) {
206  INCL_ERROR("Energy nonconservation! Energy at the beginning of the event = "
207  << finalstate->getTotalEnergyBeforeInteraction()
208  <<" and after interaction = "
209  << totalEnergy << '\n'
210  << finalstate->print());
211  }
212  }
213 
215  INCL_WARN("Useless Nucleus::propagateParticles -method called." << '\n');
216  }
217 
219  G4double totalEnergy = 0.0;
220  ParticleList const &inside = theStore->getParticles();
221  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
222  if((*p)->isNucleon()) // Ugly: we should calculate everything using total energies!
223  totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
224  else if((*p)->isResonance())
225  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::effectiveNucleonMass;
226  else if((*p)->isHyperon())
227  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::getRealMass((*p)->getType());
228  else
229  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
230  }
231 
232  return totalEnergy;
233  }
234 
236  // If the remnant consists of only one nucleon, we need to apply a special
237  // procedure to put it on mass shell.
238  if(theA==1) {
239  emitInsidePions();
241  remnant=false;
242  return;
243  }
244 
245  // Compute the recoil momentum and angular momentum
248 
249  ParticleList const &outgoing = theStore->getOutgoingParticles();
250  for(ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p) {
251  theMomentum -= (*p)->getMomentum();
252  theSpin -= (*p)->getAngularMomentum();
253  }
257  }
258 
259  // Subtract orbital angular momentum
262 
265  remnant=true;
266  }
267 
269  ThreeVector cm(0.,0.,0.);
270  G4double totalMass = 0.0;
271  ParticleList const &inside = theStore->getParticles();
272  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
273  const G4double mass = (*p)->getMass();
274  cm += (*p)->getPosition() * mass;
275  totalMass += mass;
276  }
277  cm /= totalMass;
278  return cm;
279  }
280 
282  const G4double totalEnergy = computeTotalEnergy();
283  const G4double separationEnergies = computeSeparationEnergyBalance();
284 
285  return totalEnergy - initialInternalEnergy - separationEnergies;
286  }
287 
288  std::string Nucleus::print()
289  {
290  std::stringstream ss;
291  ss << "Particles in the nucleus:" << '\n'
292  << "Inside:" << '\n';
293  G4int counter = 1;
294  ParticleList const &inside = theStore->getParticles();
295  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
296  ss << "index = " << counter << '\n'
297  << (*p)->print();
298  counter++;
299  }
300  ss <<"Outgoing:" << '\n';
301  ParticleList const &outgoing = theStore->getOutgoingParticles();
302  for(ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p)
303  ss << (*p)->print();
304 
305  return ss.str();
306  }
307 
310  ParticleList deltas;
311  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
312  if((*i)->isDelta()) deltas.push_back((*i));
313  }
314  if(deltas.empty()) return false;
315 
316  for(ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
317  INCL_DEBUG("Decay outgoing delta particle:" << '\n'
318  << (*i)->print() << '\n');
319  const ThreeVector beta = -(*i)->boostVector();
320  const G4double deltaMass = (*i)->getMass();
321 
322  // Set the delta momentum to zero and sample the decay in the CM frame.
323  // This makes life simpler if we are using real particle masses.
324  (*i)->setMomentum(ThreeVector());
325  (*i)->setEnergy((*i)->getMass());
326 
327  // Use a DecayAvatar
328  IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
329  FinalState *fs = decay->getFinalState();
330  Particle * const pion = fs->getCreatedParticles().front();
331  Particle * const nucleon = fs->getModifiedParticles().front();
332 
333  // Adjust the decay momentum if we are using the real masses
334  const G4double decayMomentum = KinematicsUtils::momentumInCM(deltaMass,
335  nucleon->getTableMass(),
336  pion->getTableMass());
337  ThreeVector newMomentum = pion->getMomentum();
338  newMomentum *= decayMomentum / newMomentum.mag();
339 
340  pion->setTableMass();
341  pion->setMomentum(newMomentum);
342  pion->adjustEnergyFromMomentum();
343  pion->setEmissionTime(nucleon->getEmissionTime());
344  pion->boost(beta);
346 
347  nucleon->setTableMass();
348  nucleon->setMomentum(-newMomentum);
349  nucleon->adjustEnergyFromMomentum();
350  nucleon->boost(beta);
351 
352  theStore->addToOutgoing(pion);
353 
354  delete fs;
355  delete decay;
356  }
357 
358  return true;
359  }
360 
362  /* If there is a pion potential, do nothing (deltas will be counted as
363  * excitation energy).
364  * If, however, the remnant is unphysical (Z<0 or Z>A), force the deltas to
365  * decay and get rid of all the pions. In case you're wondering, you can
366  * end up with Z<0 or Z>A if the remnant contains more pi- than protons or
367  * more pi+ than neutrons, respectively.
368  */
369  const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
370  if(thePotential->hasPionPotential() && !unphysicalRemnant)
371  return false;
372 
373  // Build a list of deltas (avoid modifying the list you are iterating on).
374  ParticleList const &inside = theStore->getParticles();
375  ParticleList deltas;
376  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
377  if((*i)->isDelta()) deltas.push_back((*i));
378 
379  // Loop over the deltas, make them decay
380  for(ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
381  INCL_DEBUG("Decay inside delta particle:" << '\n'
382  << (*i)->print() << '\n');
383  // Create a forced-decay avatar. Note the last boolean parameter. Note
384  // also that if the remnant is unphysical we more or less explicitly give
385  // up energy conservation and CDPP by passing a NULL pointer for the
386  // nucleus.
387  IAvatar *decay;
388  if(unphysicalRemnant) {
389  INCL_WARN("Forcing delta decay inside an unphysical remnant (A=" << theA
390  << ", Z=" << theZ << "). Might lead to energy-violation warnings."
391  << '\n');
392  decay = new DecayAvatar((*i), 0.0, NULL, true);
393  } else
394  decay = new DecayAvatar((*i), 0.0, this, true);
395  FinalState *fs = decay->getFinalState();
396 
397  // The pion can be ejected only if we managed to satisfy energy
398  // conservation and if pion emission does not lead to negative excitation
399  // energies.
400  if(fs->getValidity()==ValidFS) {
401  // Apply the final state to the nucleus
402  applyFinalState(fs);
403  }
404  delete fs;
405  delete decay;
406  }
407 
408  // If the remnant is unphysical, emit all the pions
409  if(unphysicalRemnant) {
410  INCL_DEBUG("Remnant is unphysical: Z=" << theZ << ", A=" << theA << ", emitting all the pions" << '\n');
411  emitInsidePions();
412  }
413 
414  return true;
415  }
416 
418 
419  /* Transform each strange particles into a lambda
420  * Every Kaon (KPlus and KZero) are emited
421  */
422  const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
423  if(unphysicalRemnant){
425  INCL_WARN("Remnant is unphysical: Z=" << theZ << ", A=" << theA << ", too much strange particles? -> all emit" << '\n');
426  return false;
427  }
428 
429  /* Build a list of particles with a strangeness == -1 except Lambda,
430  * and two other one for proton and neutron
431  */
432  ParticleList const &inside = theStore->getParticles();
433  ParticleList stranges;
434  ParticleList protons;
435  ParticleList neutrons;
436  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i){
437  if((*i)->isSigma() || (*i)->isAntiKaon()) stranges.push_back((*i));
438  else if((*i)->isNucleon() && (*i)->getZ() == 1) protons.push_back((*i));
439  else if((*i)->isNucleon() && (*i)->getZ() == 0) neutrons.push_back((*i));
440  }
441 
442  if((stranges.size() > protons.size()) || (stranges.size() > neutrons.size())){
443  INCL_WARN("Remnant is unphysical: Nproton=" << protons.size() << ", Nneutron=" << neutrons.size() << ", Strange particles : " << stranges.size() << '\n');
445  return false;
446  }
447 
448  // Loop over the strange particles, make them absorbe
449  ParticleIter protonIter = protons.begin();
450  ParticleIter neutronIter = neutrons.begin();
451  for(ParticleIter i=stranges.begin(), e=stranges.end(); i!=e; ++i) {
452  INCL_DEBUG("Absorbe inside strange particles:" << '\n'
453  << (*i)->print() << '\n');
454  IAvatar *decay;
455  if((*i)->getType() == SigmaMinus){
456  decay = new DecayAvatar((*protonIter), (*i), 0.0, this, true);
457  ++protonIter;
458  }
459  else if((*i)->getType() == SigmaPlus){
460  decay = new DecayAvatar((*neutronIter), (*i), 0.0, this, true);
461  ++neutronIter;
462  }
463  else if(Random::shoot()*(protons.size() + neutrons.size()) < protons.size()){
464  decay = new DecayAvatar((*protonIter), (*i), 0.0, this, true);
465  ++protonIter;
466  }
467  else {
468  decay = new DecayAvatar((*neutronIter), (*i), 0.0, this, true);
469  ++neutronIter;
470  }
471  FinalState *fs = decay->getFinalState();
472  applyFinalState(fs);
473  delete fs;
474  delete decay;
475  }
476 
477  return true;
478  }
479 
482  ParticleList pionResonances;
483  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
484 // if((*i)->isEta() || (*i)->isOmega()) pionResonances.push_back((*i));
485  if(((*i)->isEta() && timeThreshold > ParticleTable::getWidth(Eta)) || ((*i)->isOmega() && timeThreshold > ParticleTable::getWidth(Omega))) pionResonances.push_back((*i));
486  }
487  if(pionResonances.empty()) return false;
488 
489  for(ParticleIter i=pionResonances.begin(), e=pionResonances.end(); i!=e; ++i) {
490  INCL_DEBUG("Decay outgoing pionResonances particle:" << '\n'
491  << (*i)->print() << '\n');
492  const ThreeVector beta = -(*i)->boostVector();
493  const G4double pionResonanceMass = (*i)->getMass();
494 
495  // Set the pionResonance momentum to zero and sample the decay in the CM frame.
496  // This makes life simpler if we are using real particle masses.
497  (*i)->setMomentum(ThreeVector());
498  (*i)->setEnergy((*i)->getMass());
499 
500  // Use a DecayAvatar
501  IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
502  FinalState *fs = decay->getFinalState();
503 
504  Particle * const theModifiedParticle = fs->getModifiedParticles().front();
505  ParticleList const &created = fs->getCreatedParticles();
506  Particle * const theCreatedParticle1 = created.front();
507 
508  if (created.size() == 1) {
509 
510  // Adjust the decay momentum if we are using the real masses
511  const G4double decayMomentum = KinematicsUtils::momentumInCM(pionResonanceMass,theModifiedParticle->getTableMass(),theCreatedParticle1->getTableMass());
512  ThreeVector newMomentum = theCreatedParticle1->getMomentum();
513  newMomentum *= decayMomentum / newMomentum.mag();
514 
515  theCreatedParticle1->setTableMass();
516  theCreatedParticle1->setMomentum(newMomentum);
517  theCreatedParticle1->adjustEnergyFromMomentum();
518  //theCreatedParticle1->setEmissionTime(nucleon->getEmissionTime());
519  theCreatedParticle1->boost(beta);
520  theCreatedParticle1->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
521 
522  theModifiedParticle->setTableMass();
523  theModifiedParticle->setMomentum(-newMomentum);
524  theModifiedParticle->adjustEnergyFromMomentum();
525  theModifiedParticle->boost(beta);
526 
527  theStore->addToOutgoing(theCreatedParticle1);
528  }
529  else if (created.size() == 2) {
530  Particle * const theCreatedParticle2 = created.back();
531 
532  theCreatedParticle1->boost(beta);
533  theCreatedParticle1->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
534  theCreatedParticle2->boost(beta);
535  theCreatedParticle2->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
536  theModifiedParticle->boost(beta);
537 
538  theStore->addToOutgoing(theCreatedParticle1);
539  theStore->addToOutgoing(theCreatedParticle2);
540  }
541  else {
542  INCL_ERROR("Wrong number (< 2) of created particles during the decay of a pion resonance");
543  }
544  delete fs;
545  delete decay;
546  }
547 
548  return true;
549  }
550 
553  ParticleList neutralsigma;
554  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
555  if((*i)->getType() == SigmaZero && timeThreshold > ParticleTable::getWidth(SigmaZero)) neutralsigma.push_back((*i));
556  }
557  if(neutralsigma.empty()) return false;
558 
559  for(ParticleIter i=neutralsigma.begin(), e=neutralsigma.end(); i!=e; ++i) {
560  INCL_DEBUG("Decay outgoing neutral sigma:" << '\n'
561  << (*i)->print() << '\n');
562  const ThreeVector beta = -(*i)->boostVector();
563  const G4double neutralsigmaMass = (*i)->getMass();
564 
565  // Set the neutral sigma momentum to zero and sample the decay in the CM frame.
566  // This makes life simpler if we are using real particle masses.
567  (*i)->setMomentum(ThreeVector());
568  (*i)->setEnergy((*i)->getMass());
569 
570  // Use a DecayAvatar
571  IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
572  FinalState *fs = decay->getFinalState();
573 
574  Particle * const theModifiedParticle = fs->getModifiedParticles().front();
575  ParticleList const &created = fs->getCreatedParticles();
576  Particle * const theCreatedParticle = created.front();
577 
578  if (created.size() == 1) {
579 
580  // Adjust the decay momentum if we are using the real masses
581  const G4double decayMomentum = KinematicsUtils::momentumInCM(neutralsigmaMass,theModifiedParticle->getTableMass(),theCreatedParticle->getTableMass());
582  ThreeVector newMomentum = theCreatedParticle->getMomentum();
583  newMomentum *= decayMomentum / newMomentum.mag();
584 
585  theCreatedParticle->setTableMass();
586  theCreatedParticle->setMomentum(newMomentum);
587  theCreatedParticle->adjustEnergyFromMomentum();
588  theCreatedParticle->boost(beta);
589  theCreatedParticle->setBiasCollisionVector(theModifiedParticle->getBiasCollisionVector());
590 
591  theModifiedParticle->setTableMass();
592  theModifiedParticle->setMomentum(-newMomentum);
593  theModifiedParticle->adjustEnergyFromMomentum();
594  theModifiedParticle->boost(beta);
595 
596  theStore->addToOutgoing(theCreatedParticle);
597  }
598  else {
599  INCL_ERROR("Wrong number (!= 1) of created particles during the decay of a sigma zero");
600  }
601  delete fs;
602  delete decay;
603  }
604 
605  return true;
606  }
607 
610  ParticleList neutralkaon;
611  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
612  if((*i)->getType() == KZero || (*i)->getType() == KZeroBar) neutralkaon.push_back((*i));
613  }
614  if(neutralkaon.empty()) return false;
615 
616  for(ParticleIter i=neutralkaon.begin(), e=neutralkaon.end(); i!=e; ++i) {
617  INCL_DEBUG("Transform outgoing neutral kaon:" << '\n'
618  << (*i)->print() << '\n');
619 
620  // Use a DecayAvatar
621  IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
622  FinalState *fs = decay->getFinalState();
623 
624  delete fs;
625  delete decay;
626  }
627 
628  return true;
629  }
630 
633  ParticleList clusters;
634  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
635  if((*i)->isCluster()) clusters.push_back((*i));
636  }
637  if(clusters.empty()) return false;
638 
639  for(ParticleIter i=clusters.begin(), e=clusters.end(); i!=e; ++i) {
640  Cluster *cluster = dynamic_cast<Cluster*>(*i); // Can't avoid using a cast here
641 // assert(cluster);
642 #ifdef INCLXX_IN_GEANT4_MODE
643  if(!cluster)
644  continue;
645 #endif
646  cluster->deleteParticles(); // Don't need them
647  ParticleList decayProducts = ClusterDecay::decay(cluster);
648  for(ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j){
649  (*j)->setBiasCollisionVector(cluster->getBiasCollisionVector());
650  theStore->addToOutgoing(*j);
651  }
652  }
653  return true;
654  }
655 
657  // Do the phase-space decay only if Z=0 or N=0
658  if(theA<=1 || (theZ!=0 && (theA+theS)!=theZ))
659  return false;
660 
661  ParticleList decayProducts = ClusterDecay::decay(this);
662  for(ParticleIter j=decayProducts.begin(), e=decayProducts.end(); j!=e; ++j){
663  (*j)->setBiasCollisionVector(this->getBiasCollisionVector());
664  theStore->addToOutgoing(*j);
665  }
666 
667  return true;
668  }
669 
671  /* Forcing emissions of all pions in the nucleus. This probably violates
672  * energy conservation (although the computation of the recoil kinematics
673  * might sweep this under the carpet).
674  */
675  INCL_WARN("Forcing emissions of all pions in the nucleus." << '\n');
676 
677  // Emit the pions with this kinetic energy
678  const G4double tinyPionEnergy = 0.1; // MeV
679 
680  // Push out the emitted pions
681  ParticleList const &inside = theStore->getParticles();
682  ParticleList toEject;
683  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
684  if((*i)->isPion()) {
685  Particle * const thePion = *i;
686  INCL_DEBUG("Forcing emission of the following particle: "
687  << thePion->print() << '\n');
689  // Correction for real masses
690  const G4double theQValueCorrection = thePion->getEmissionQValueCorrection(theA,theZ,theS);
691  const G4double kineticEnergyOutside = thePion->getKineticEnergy() - thePion->getPotentialEnergy() + theQValueCorrection;
692  thePion->setTableMass();
693  if(kineticEnergyOutside > 0.0)
694  thePion->setEnergy(thePion->getMass()+kineticEnergyOutside);
695  else
696  thePion->setEnergy(thePion->getMass()+tinyPionEnergy);
697  thePion->adjustMomentumFromEnergy();
698  thePion->setPotentialEnergy(0.);
699  theZ -= thePion->getZ();
700  toEject.push_back(thePion);
701  }
702  }
703  for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
705  theStore->addToOutgoing(*i);
706  (*i)->setParticleBias(Particle::getTotalBias());
707  }
708  }
709 
711  /* Forcing emissions of Sigmas and antiKaons.
712  * This probably violates energy conservation
713  * (although the computation of the recoil kinematics
714  * might sweep this under the carpet).
715  */
716  INCL_DEBUG("Forcing emissions of all strange particles in the nucleus." << '\n');
717 
718  // Emit the strange particles with this kinetic energy
719  const G4double tinyEnergy = 0.1; // MeV
720 
721  // Push out the emitted strange particles
722  ParticleList const &inside = theStore->getParticles();
723  ParticleList toEject;
724  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
725  if((*i)->isSigma() || (*i)->isAntiKaon()) {
726  Particle * const theParticle = *i;
727  INCL_DEBUG("Forcing emission of the following particle: "
728  << theParticle->print() << '\n');
729  theParticle->setEmissionTime(theStore->getBook().getCurrentTime());
730  // Correction for real masses
731  const G4double theQValueCorrection = theParticle->getEmissionQValueCorrection(theA,theZ,theS); // Does it work for strange particles? should be check
732  const G4double kineticEnergyOutside = theParticle->getKineticEnergy() - theParticle->getPotentialEnergy() + theQValueCorrection;
733  theParticle->setTableMass();
734  if(kineticEnergyOutside > 0.0)
735  theParticle->setEnergy(theParticle->getMass()+kineticEnergyOutside);
736  else
737  theParticle->setEnergy(theParticle->getMass()+tinyEnergy);
738  theParticle->adjustMomentumFromEnergy();
739  theParticle->setPotentialEnergy(0.);
740  theA -= theParticle->getA();
741  theZ -= theParticle->getZ();
742  theS -= theParticle->getS();
743  toEject.push_back(theParticle);
744  }
745  }
746  for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
748  theStore->addToOutgoing(*i);
749  (*i)->setParticleBias(Particle::getTotalBias());
750  }
751  }
752 
754  /* Forcing emissions of all Lambda in the nucleus.
755  * This probably violates energy conservation
756  * (although the computation of the recoil kinematics
757  * might sweep this under the carpet).
758  */
759  INCL_DEBUG("Forcing emissions of all Lambda in the nucleus." << '\n');
760 
761  // Emit the Lambda with this kinetic energy
762  const G4double tinyEnergy = 0.1; // MeV
763 
764  // Push out the emitted Lambda
765  ParticleList const &inside = theStore->getParticles();
766  ParticleList toEject;
767  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
768  if((*i)->isLambda()) {
769  Particle * const theLambda = *i;
770  INCL_DEBUG("Forcing emission of the following particle: "
771  << theLambda->print() << '\n');
773  // Correction for real masses
774  const G4double theQValueCorrection = theLambda->getEmissionQValueCorrection(theA,theZ,theS); // Does it work for strange particles? Should be check
775  const G4double kineticEnergyOutside = theLambda->getKineticEnergy() - theLambda->getPotentialEnergy() + theQValueCorrection;
776  theLambda->setTableMass();
777  if(kineticEnergyOutside > 0.0)
778  theLambda->setEnergy(theLambda->getMass()+kineticEnergyOutside);
779  else
780  theLambda->setEnergy(theLambda->getMass()+tinyEnergy);
781  theLambda->adjustMomentumFromEnergy();
782  theLambda->setPotentialEnergy(0.);
783  theA -= theLambda->getA();
784  theS -= theLambda->getS();
785  toEject.push_back(theLambda);
786  }
787  }
788  for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
790  theStore->addToOutgoing(*i);
791  (*i)->setParticleBias(Particle::getTotalBias());
792  }
793  return toEject.size();
794  }
795 
797  /* Forcing emissions of all Kaon (not antiKaons) in the nucleus.
798  * This probably violates energy conservation
799  * (although the computation of the recoil kinematics
800  * might sweep this under the carpet).
801  */
802  INCL_DEBUG("Forcing emissions of all Kaon in the nucleus." << '\n');
803 
804  // Emit the Kaon with this kinetic energy (not supposed to append
805  const G4double tinyEnergy = 0.1; // MeV
806 
807  // Push out the emitted kaon
808  ParticleList const &inside = theStore->getParticles();
809  ParticleList toEject;
810  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
811  if((*i)->isKaon()) {
812  Particle * const theKaon = *i;
813  INCL_DEBUG("Forcing emission of the following particle: "
814  << theKaon->print() << '\n');
816  // Correction for real masses
817  const G4double theQValueCorrection = theKaon->getEmissionQValueCorrection(theA,theZ,theS);
818  const G4double kineticEnergyOutside = theKaon->getKineticEnergy() - theKaon->getPotentialEnergy() + theQValueCorrection;
819  theKaon->setTableMass();
820  if(kineticEnergyOutside > 0.0)
821  theKaon->setEnergy(theKaon->getMass()+kineticEnergyOutside);
822  else
823  theKaon->setEnergy(theKaon->getMass()+tinyEnergy);
824  theKaon->adjustMomentumFromEnergy();
825  theKaon->setPotentialEnergy(0.);
826  theZ -= theKaon->getZ();
827  theS -= theKaon->getS();
828  toEject.push_back(theKaon);
829  }
830  }
831  for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
833  theStore->addToOutgoing(*i);
834  (*i)->setParticleBias(Particle::getTotalBias());
835  }
836  theNKaon -= 1;
837  return toEject.size() != 0;
838  }
839 
841 
842  Book const &theBook = theStore->getBook();
843  const G4int nEventCollisions = theBook.getAcceptedCollisions();
844  const G4int nEventDecays = theBook.getAcceptedDecays();
845  const G4int nEventClusters = theBook.getEmittedClusters();
846  if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
847  return true;
848 
849  return false;
850 
851  }
852 
854  // We should be here only if the nucleus contains only one nucleon
855 // assert(theStore->getParticles().size()==1);
856 
857  // No excitation energy!
858  theExcitationEnergy = 0.0;
859 
860  // Move the nucleon to the outgoing list
861  Particle *remN = theStore->getParticles().front();
862  theA -= remN->getA();
863  theZ -= remN->getZ();
864  theS -= remN->getS();
866  theStore->addToOutgoing(remN);
868 
869  // Treat the special case of a remaining delta
870  if(remN->isDelta()) {
871  IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
872  FinalState *fs = decay->getFinalState();
873  // Eject the pion
874  ParticleList const &created = fs->getCreatedParticles();
875  for(ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
876  theStore->addToOutgoing(*j);
877  delete fs;
878  delete decay;
879  }
880 
881  // Do different things depending on how many outgoing particles we have
882  ParticleList const &outgoing = theStore->getOutgoingParticles();
883  if(outgoing.size() == 2) {
884 
885  INCL_DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << '\n');
886 
887  // Can apply exact 2-body kinematics here. Keep the CM emission angle of
888  // the first particle.
889  Particle *p1 = outgoing.front(), *p2 = outgoing.back();
890  const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
891  // Boost to the initial CM
892  p1->boost(aBoostVector);
893  const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
894  const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
895  const G4double scale = pcm/(p1->getMomentum().mag());
896  // Reset the momenta
897  p1->setMomentum(p1->getMomentum()*scale);
898  p2->setMomentum(-p1->getMomentum());
900  p2->adjustEnergyFromMomentum();
901  // Unboost
902  p1->boost(-aBoostVector);
903  p2->boost(-aBoostVector);
904 
905  } else {
906 
907  INCL_DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << '\n');
908 
909  const G4int maxIterations=8;
910  G4double totalEnergy, energyScale;
911  G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
912  ThreeVector totalMomentum, deltaP;
913  std::vector<ThreeVector> minMomenta; // use it to store the particle momenta that minimize the merit function
914 
915  // Reserve the vector size
916  minMomenta.reserve(outgoing.size());
917 
918  // Compute the initial total momentum
919  totalMomentum.setX(0.0);
920  totalMomentum.setY(0.0);
921  totalMomentum.setZ(0.0);
922  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
923  totalMomentum += (*i)->getMomentum();
924 
925  // Compute the initial total energy
926  totalEnergy = 0.0;
927  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
928  totalEnergy += (*i)->getEnergy();
929 
930  // Iterative algorithm starts here:
931  for(G4int iterations=0; iterations < maxIterations; ++iterations) {
932 
933  // Save the old merit-function values
934  oldOldOldVal = oldOldVal;
935  oldOldVal = oldVal;
936  oldVal = val;
937 
938  if(iterations%2 == 0) {
939  INCL_DEBUG("Momentum step" << '\n');
940  // Momentum step: modify all the particle momenta
941  deltaP = incomingMomentum - totalMomentum;
942  G4double pOldTot = 0.0;
943  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
944  pOldTot += (*i)->getMomentum().mag();
945  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
946  const ThreeVector mom = (*i)->getMomentum();
947  (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
948  (*i)->adjustEnergyFromMomentum();
949  }
950  } else {
951  INCL_DEBUG("Energy step" << '\n');
952  // Energy step: modify all the particle momenta
953  energyScale = initialEnergy/totalEnergy;
954  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
955  const ThreeVector mom = (*i)->getMomentum();
956  G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
957  if(pScale>0) {
958  (*i)->setEnergy((*i)->getEnergy()*energyScale);
959  (*i)->adjustMomentumFromEnergy();
960  }
961  }
962  }
963 
964  // Compute the current total momentum and energy
965  totalMomentum.setX(0.0);
966  totalMomentum.setY(0.0);
967  totalMomentum.setZ(0.0);
968  totalEnergy = 0.0;
969  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
970  totalMomentum += (*i)->getMomentum();
971  totalEnergy += (*i)->getEnergy();
972  }
973 
974  // Merit factor
975  val = std::pow(totalEnergy - initialEnergy,2) +
976  0.25*(totalMomentum - incomingMomentum).mag2();
977  INCL_DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
978 
979  // Store the minimum
980  if(val < oldVal) {
981  INCL_DEBUG("New minimum found, storing the particle momenta" << '\n');
982  minMomenta.clear();
983  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
984  minMomenta.push_back((*i)->getMomentum());
985  }
986 
987  // Stop the algorithm if the search diverges
988  if(val > oldOldVal && oldVal > oldOldOldVal) {
989  INCL_DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
990  break;
991  }
992  }
993 
994  // We should have made at least one successful iteration here
995 // assert(minMomenta.size()==outgoing.size());
996 
997  // Apply the optimal momenta
998  INCL_DEBUG("Applying the solution" << '\n');
999  std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
1000  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
1001  (*i)->setMomentum(*v);
1002  (*i)->adjustEnergyFromMomentum();
1003  INCL_DATABLOCK((*i)->print());
1004  }
1005 
1006  }
1007 
1008  }
1009 
1011  eventInfo->nParticles = 0;
1012  G4bool isNucleonAbsorption = false;
1013 
1014  G4bool isPionAbsorption = false;
1015  // It is possible to have pion absorption event only if the
1016  // projectile is pion.
1017  if(eventInfo->projectileType == PiPlus ||
1018  eventInfo->projectileType == PiMinus ||
1019  eventInfo->projectileType == PiZero) {
1020  isPionAbsorption = true;
1021  }
1022 
1023  // Forced CN
1024  eventInfo->forcedCompoundNucleus = tryCN;
1025 
1026  // Outgoing particles
1027  ParticleList const &outgoingParticles = getStore()->getOutgoingParticles();
1028 
1029  // Check if we have a nucleon absorption event: nucleon projectile
1030  // and no ejected particles.
1031  if(outgoingParticles.size() == 0 &&
1032  (eventInfo->projectileType == Proton ||
1033  eventInfo->projectileType == Neutron)) {
1034  isNucleonAbsorption = true;
1035  }
1036 
1037  // Reset the remnant counter
1038  eventInfo->nRemnants = 0;
1039  eventInfo->history.clear();
1040 
1041  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1042  // We have a pion absorption event only if the projectile is
1043  // pion and there are no ejected pions.
1044  if(isPionAbsorption) {
1045  if((*i)->isPion()) {
1046  isPionAbsorption = false;
1047  }
1048  }
1049 
1050  eventInfo->ParticleBias[eventInfo->nParticles] = (*i)->getParticleBias();
1051 
1052  eventInfo->A[eventInfo->nParticles] = (*i)->getA();
1053  eventInfo->Z[eventInfo->nParticles] = (*i)->getZ();
1054  eventInfo->S[eventInfo->nParticles] = (*i)->getS();
1055  eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
1056  eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
1057  ThreeVector mom = (*i)->getMomentum();
1058  eventInfo->px[eventInfo->nParticles] = mom.getX();
1059  eventInfo->py[eventInfo->nParticles] = mom.getY();
1060  eventInfo->pz[eventInfo->nParticles] = mom.getZ();
1061  eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
1062  eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
1063  eventInfo->origin[eventInfo->nParticles] = -1;
1064  eventInfo->history.push_back("");
1065  if ((*i)->getType() != Composite) {
1066  ParticleSpecies pt((*i)->getType());
1067  eventInfo->PDGCode[eventInfo->nParticles] = pt.getPDGCode();
1068  }
1069  else {
1070  ParticleSpecies pt((*i)->getA(), (*i)->getZ(), (*i)->getS());
1071  eventInfo->PDGCode[eventInfo->nParticles] = pt.getPDGCode();
1072  }
1073  eventInfo->nParticles++;
1074  }
1075  eventInfo->nucleonAbsorption = isNucleonAbsorption;
1076  eventInfo->pionAbsorption = isPionAbsorption;
1077  eventInfo->nCascadeParticles = eventInfo->nParticles;
1078 
1079  // Projectile-like remnant characteristics
1081  eventInfo->ARem[eventInfo->nRemnants] = theProjectileRemnant->getA();
1082  eventInfo->ZRem[eventInfo->nRemnants] = theProjectileRemnant->getZ();
1083  eventInfo->SRem[eventInfo->nRemnants] = theProjectileRemnant->getS();
1085  if(std::abs(eStar)<1E-10)
1086  eStar = 0.0; // blame rounding and set the excitation energy to zero
1087  eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
1088  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1089  INCL_WARN("Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << '\n');
1090  }
1091  const ThreeVector &spin = theProjectileRemnant->getSpin();
1092  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1093  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
1094  } else { // odd-A nucleus
1095  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
1096  }
1097  eventInfo->EKinRem[eventInfo->nRemnants] = theProjectileRemnant->getKineticEnergy();
1099  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1100  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1101  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1102  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
1103  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
1104  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
1105  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1106  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1107  eventInfo->nRemnants++;
1108  }
1109 
1110  // Target-like remnant characteristics
1111  if(hasRemnant()) {
1112  eventInfo->ARem[eventInfo->nRemnants] = getA();
1113  eventInfo->ZRem[eventInfo->nRemnants] = getZ();
1114  eventInfo->SRem[eventInfo->nRemnants] = getS();
1115  eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
1116  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1117  INCL_WARN("Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << " eventNumber=" << eventInfo->eventNumber << '\n');
1118  }
1119  const ThreeVector &spin = getSpin();
1120  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1121  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
1122  } else { // odd-A nucleus
1123  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
1124  }
1125  eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
1126  const ThreeVector &mom = getMomentum();
1127  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1128  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1129  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1130  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
1131  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
1132  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
1133  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1134  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1135  eventInfo->nRemnants++;
1136  }
1137 
1138  // Global counters, flags, etc.
1139  Book const &theBook = theStore->getBook();
1140  eventInfo->nCollisions = theBook.getAcceptedCollisions();
1141  eventInfo->nBlockedCollisions = theBook.getBlockedCollisions();
1142  eventInfo->nDecays = theBook.getAcceptedDecays();
1143  eventInfo->nBlockedDecays = theBook.getBlockedDecays();
1144  eventInfo->firstCollisionTime = theBook.getFirstCollisionTime();
1145  eventInfo->firstCollisionXSec = theBook.getFirstCollisionXSec();
1148  eventInfo->firstCollisionIsElastic = theBook.getFirstCollisionIsElastic();
1149  eventInfo->nReflectionAvatars = theBook.getAvatars(SurfaceAvatarType);
1150  eventInfo->nCollisionAvatars = theBook.getAvatars(CollisionAvatarType);
1151  eventInfo->nDecayAvatars = theBook.getAvatars(DecayAvatarType);
1153  }
1154 
1155  Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
1156  ConservationBalance theBalance;
1157  // Initialise balance variables with the incoming values
1158  theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
1159  theBalance.A = theEventInfo.Ap + theEventInfo.At;
1160  theBalance.S = theEventInfo.Sp + theEventInfo.St;
1161 
1162  theBalance.energy = getInitialEnergy();
1163  theBalance.momentum = getIncomingMomentum();
1164 
1165  // Process outgoing particles
1166  ParticleList const &outgoingParticles = theStore->getOutgoingParticles();
1167  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1168  theBalance.Z -= (*i)->getZ();
1169  theBalance.A -= (*i)->getA();
1170  theBalance.S -= (*i)->getS();
1171  // For outgoing clusters, the total energy automatically includes the
1172  // excitation energy
1173  theBalance.energy -= (*i)->getEnergy(); // Note that outgoing particles should have the real mass
1174  theBalance.momentum -= (*i)->getMomentum();
1175  }
1176 
1177  // Projectile-like remnant contribution, if present
1179  theBalance.Z -= theProjectileRemnant->getZ();
1180  theBalance.A -= theProjectileRemnant->getA();
1181  theBalance.S -= theProjectileRemnant->getS();
1184  theBalance.energy -= theProjectileRemnant->getKineticEnergy();
1185  theBalance.momentum -= theProjectileRemnant->getMomentum();
1186  }
1187 
1188  // Target-like remnant contribution, if present
1189  if(hasRemnant()) {
1190  theBalance.Z -= getZ();
1191  theBalance.A -= getA();
1192  theBalance.S -= getS();
1193  theBalance.energy -= ParticleTable::getTableMass(getA(),getZ(),getS()) +
1195  if(afterRecoil)
1196  theBalance.energy -= getKineticEnergy();
1197  theBalance.momentum -= getMomentum();
1198  }
1199 
1200  return theBalance;
1201  }
1202 
1209  }
1210 
1211  void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
1212  // Deal with the projectile remnant
1213  const G4int prA = theProjectileRemnant->getA();
1214  if(prA>=1) {
1215  // Set the mass
1217  theProjectileRemnant->setMass(aMass);
1218 
1219  // Compute the excitation energy from the invariant mass
1220  const G4double anExcitationEnergy = aMass
1222 
1223  // Set the excitation energy
1224  theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
1225 
1226  // No spin!
1228 
1229  // Set the emission time
1230  theProjectileRemnant->setEmissionTime(anEmissionTime);
1231  }
1232  }
1233 
1234 }
1235 
1236 #endif