Geant4  v4-10.4-release
 모두 클래스 네임스페이스들 파일들 함수 변수 타입정의 열거형 타입 열거형 멤버 Friends 매크로 그룹들 페이지들
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)->isLambda())
227  totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
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 Z=A
658  if(theA<=1 || (theZ!=0 && theA!=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);
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 all strange particles exept Lambda particles in the nucleus.
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); // 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); // 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 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);
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  return toEject.size() != 0;
837  }
838 
840 
841  Book const &theBook = theStore->getBook();
842  const G4int nEventCollisions = theBook.getAcceptedCollisions();
843  const G4int nEventDecays = theBook.getAcceptedDecays();
844  const G4int nEventClusters = theBook.getEmittedClusters();
845  if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
846  return true;
847 
848  return false;
849 
850  }
851 
853  // We should be here only if the nucleus contains only one nucleon
854 // assert(theStore->getParticles().size()==1);
855 
856  // No excitation energy!
857  theExcitationEnergy = 0.0;
858 
859  // Move the nucleon to the outgoing list
860  Particle *remN = theStore->getParticles().front();
861  theA -= remN->getA();
862  theZ -= remN->getZ();
863  theS -= remN->getS();
865  theStore->addToOutgoing(remN);
867 
868  // Treat the special case of a remaining delta
869  if(remN->isDelta()) {
870  IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
871  FinalState *fs = decay->getFinalState();
872  // Eject the pion
873  ParticleList const &created = fs->getCreatedParticles();
874  for(ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
875  theStore->addToOutgoing(*j);
876  delete fs;
877  delete decay;
878  }
879 
880  // Do different things depending on how many outgoing particles we have
881  ParticleList const &outgoing = theStore->getOutgoingParticles();
882  if(outgoing.size() == 2) {
883 
884  INCL_DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << '\n');
885 
886  // Can apply exact 2-body kinematics here. Keep the CM emission angle of
887  // the first particle.
888  Particle *p1 = outgoing.front(), *p2 = outgoing.back();
889  const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
890  // Boost to the initial CM
891  p1->boost(aBoostVector);
892  const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
893  const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
894  const G4double scale = pcm/(p1->getMomentum().mag());
895  // Reset the momenta
896  p1->setMomentum(p1->getMomentum()*scale);
897  p2->setMomentum(-p1->getMomentum());
899  p2->adjustEnergyFromMomentum();
900  // Unboost
901  p1->boost(-aBoostVector);
902  p2->boost(-aBoostVector);
903 
904  } else {
905 
906  INCL_DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << '\n');
907 
908  const G4int maxIterations=8;
909  G4double totalEnergy, energyScale;
910  G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
911  ThreeVector totalMomentum, deltaP;
912  std::vector<ThreeVector> minMomenta; // use it to store the particle momenta that minimize the merit function
913 
914  // Reserve the vector size
915  minMomenta.reserve(outgoing.size());
916 
917  // Compute the initial total momentum
918  totalMomentum.setX(0.0);
919  totalMomentum.setY(0.0);
920  totalMomentum.setZ(0.0);
921  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
922  totalMomentum += (*i)->getMomentum();
923 
924  // Compute the initial total energy
925  totalEnergy = 0.0;
926  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
927  totalEnergy += (*i)->getEnergy();
928 
929  // Iterative algorithm starts here:
930  for(G4int iterations=0; iterations < maxIterations; ++iterations) {
931 
932  // Save the old merit-function values
933  oldOldOldVal = oldOldVal;
934  oldOldVal = oldVal;
935  oldVal = val;
936 
937  if(iterations%2 == 0) {
938  INCL_DEBUG("Momentum step" << '\n');
939  // Momentum step: modify all the particle momenta
940  deltaP = incomingMomentum - totalMomentum;
941  G4double pOldTot = 0.0;
942  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
943  pOldTot += (*i)->getMomentum().mag();
944  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
945  const ThreeVector mom = (*i)->getMomentum();
946  (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
947  (*i)->adjustEnergyFromMomentum();
948  }
949  } else {
950  INCL_DEBUG("Energy step" << '\n');
951  // Energy step: modify all the particle momenta
952  energyScale = initialEnergy/totalEnergy;
953  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
954  const ThreeVector mom = (*i)->getMomentum();
955  G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
956  if(pScale>0) {
957  (*i)->setEnergy((*i)->getEnergy()*energyScale);
958  (*i)->adjustMomentumFromEnergy();
959  }
960  }
961  }
962 
963  // Compute the current total momentum and energy
964  totalMomentum.setX(0.0);
965  totalMomentum.setY(0.0);
966  totalMomentum.setZ(0.0);
967  totalEnergy = 0.0;
968  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
969  totalMomentum += (*i)->getMomentum();
970  totalEnergy += (*i)->getEnergy();
971  }
972 
973  // Merit factor
974  val = std::pow(totalEnergy - initialEnergy,2) +
975  0.25*(totalMomentum - incomingMomentum).mag2();
976  INCL_DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
977 
978  // Store the minimum
979  if(val < oldVal) {
980  INCL_DEBUG("New minimum found, storing the particle momenta" << '\n');
981  minMomenta.clear();
982  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
983  minMomenta.push_back((*i)->getMomentum());
984  }
985 
986  // Stop the algorithm if the search diverges
987  if(val > oldOldVal && oldVal > oldOldOldVal) {
988  INCL_DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
989  break;
990  }
991  }
992 
993  // We should have made at least one successful iteration here
994 // assert(minMomenta.size()==outgoing.size());
995 
996  // Apply the optimal momenta
997  INCL_DEBUG("Applying the solution" << '\n');
998  std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
999  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
1000  (*i)->setMomentum(*v);
1001  (*i)->adjustEnergyFromMomentum();
1002  INCL_DATABLOCK((*i)->print());
1003  }
1004 
1005  }
1006 
1007  }
1008 
1010  eventInfo->nParticles = 0;
1011  G4bool isNucleonAbsorption = false;
1012 
1013  G4bool isPionAbsorption = false;
1014  // It is possible to have pion absorption event only if the
1015  // projectile is pion.
1016  if(eventInfo->projectileType == PiPlus ||
1017  eventInfo->projectileType == PiMinus ||
1018  eventInfo->projectileType == PiZero) {
1019  isPionAbsorption = true;
1020  }
1021 
1022  // Forced CN
1023  eventInfo->forcedCompoundNucleus = tryCN;
1024 
1025  // Outgoing particles
1026  ParticleList const &outgoingParticles = getStore()->getOutgoingParticles();
1027 
1028  // Check if we have a nucleon absorption event: nucleon projectile
1029  // and no ejected particles.
1030  if(outgoingParticles.size() == 0 &&
1031  (eventInfo->projectileType == Proton ||
1032  eventInfo->projectileType == Neutron)) {
1033  isNucleonAbsorption = true;
1034  }
1035 
1036  // Reset the remnant counter
1037  eventInfo->nRemnants = 0;
1038  eventInfo->history.clear();
1039 
1040  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1041  // We have a pion absorption event only if the projectile is
1042  // pion and there are no ejected pions.
1043  if(isPionAbsorption) {
1044  if((*i)->isPion()) {
1045  isPionAbsorption = false;
1046  }
1047  }
1048 
1049  eventInfo->ParticleBias[eventInfo->nParticles] = (*i)->getParticleBias();
1050 
1051  eventInfo->A[eventInfo->nParticles] = (*i)->getA();
1052  eventInfo->Z[eventInfo->nParticles] = (*i)->getZ();
1053  eventInfo->S[eventInfo->nParticles] = (*i)->getS();
1054  eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
1055  eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
1056  ThreeVector mom = (*i)->getMomentum();
1057  eventInfo->px[eventInfo->nParticles] = mom.getX();
1058  eventInfo->py[eventInfo->nParticles] = mom.getY();
1059  eventInfo->pz[eventInfo->nParticles] = mom.getZ();
1060  eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
1061  eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
1062  eventInfo->origin[eventInfo->nParticles] = -1;
1063  eventInfo->history.push_back("");
1064  if ((*i)->getType() != Composite) {
1065  ParticleSpecies pt((*i)->getType());
1066  eventInfo->PDGCode[eventInfo->nParticles] = pt.getPDGCode();
1067  }
1068  else {
1069  ParticleSpecies pt((*i)->getA(), (*i)->getZ(), (*i)->getS());
1070  eventInfo->PDGCode[eventInfo->nParticles] = pt.getPDGCode();
1071  }
1072  eventInfo->nParticles++;
1073  }
1074  eventInfo->nucleonAbsorption = isNucleonAbsorption;
1075  eventInfo->pionAbsorption = isPionAbsorption;
1076  eventInfo->nCascadeParticles = eventInfo->nParticles;
1077 
1078  // Projectile-like remnant characteristics
1080  eventInfo->ARem[eventInfo->nRemnants] = theProjectileRemnant->getA();
1081  eventInfo->ZRem[eventInfo->nRemnants] = theProjectileRemnant->getZ();
1082  eventInfo->SRem[eventInfo->nRemnants] = theProjectileRemnant->getS();
1084  if(std::abs(eStar)<1E-10)
1085  eStar = 0.0; // blame rounding and set the excitation energy to zero
1086  eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
1087  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1088  INCL_WARN("Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << '\n');
1089  }
1090  const ThreeVector &spin = theProjectileRemnant->getSpin();
1091  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1092  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
1093  } else { // odd-A nucleus
1094  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
1095  }
1096  eventInfo->EKinRem[eventInfo->nRemnants] = theProjectileRemnant->getKineticEnergy();
1098  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1099  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1100  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1101  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
1102  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
1103  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
1104  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1105  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1106  eventInfo->nRemnants++;
1107  }
1108 
1109  // Target-like remnant characteristics
1110  if(hasRemnant()) {
1111  eventInfo->ARem[eventInfo->nRemnants] = getA();
1112  eventInfo->ZRem[eventInfo->nRemnants] = getZ();
1113  eventInfo->SRem[eventInfo->nRemnants] = getS();
1114  eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
1115  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1116  INCL_WARN("Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << " eventNumber=" << eventInfo->eventNumber << '\n');
1117  }
1118  const ThreeVector &spin = getSpin();
1119  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1120  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
1121  } else { // odd-A nucleus
1122  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
1123  }
1124  eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
1125  const ThreeVector &mom = getMomentum();
1126  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1127  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1128  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1129  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
1130  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
1131  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
1132  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1133  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1134  eventInfo->nRemnants++;
1135  }
1136 
1137  // Global counters, flags, etc.
1138  Book const &theBook = theStore->getBook();
1139  eventInfo->nCollisions = theBook.getAcceptedCollisions();
1140  eventInfo->nBlockedCollisions = theBook.getBlockedCollisions();
1141  eventInfo->nDecays = theBook.getAcceptedDecays();
1142  eventInfo->nBlockedDecays = theBook.getBlockedDecays();
1143  eventInfo->firstCollisionTime = theBook.getFirstCollisionTime();
1144  eventInfo->firstCollisionXSec = theBook.getFirstCollisionXSec();
1147  eventInfo->firstCollisionIsElastic = theBook.getFirstCollisionIsElastic();
1148  eventInfo->nReflectionAvatars = theBook.getAvatars(SurfaceAvatarType);
1149  eventInfo->nCollisionAvatars = theBook.getAvatars(CollisionAvatarType);
1150  eventInfo->nDecayAvatars = theBook.getAvatars(DecayAvatarType);
1152  }
1153 
1154  Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
1155  ConservationBalance theBalance;
1156  // Initialise balance variables with the incoming values
1157  theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
1158  theBalance.A = theEventInfo.Ap + theEventInfo.At;
1159  theBalance.S = theEventInfo.Sp + theEventInfo.St;
1160 
1161  theBalance.energy = getInitialEnergy();
1162  theBalance.momentum = getIncomingMomentum();
1163 
1164  // Process outgoing particles
1165  ParticleList const &outgoingParticles = theStore->getOutgoingParticles();
1166  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
1167  theBalance.Z -= (*i)->getZ();
1168  theBalance.A -= (*i)->getA();
1169  theBalance.S += (*i)->getS();
1170  // For outgoing clusters, the total energy automatically includes the
1171  // excitation energy
1172  theBalance.energy -= (*i)->getEnergy(); // Note that outgoing particles should have the real mass
1173  theBalance.momentum -= (*i)->getMomentum();
1174  }
1175 
1176  // Projectile-like remnant contribution, if present
1178  theBalance.Z -= theProjectileRemnant->getZ();
1179  theBalance.A -= theProjectileRemnant->getA();
1180  theBalance.S += theProjectileRemnant->getS();
1183  if(theProjectileRemnant->getS() != 0) theBalance.energy -= theProjectileRemnant->getS()*(ParticleTable::effectiveNucleonMass - ParticleTable::effectiveLambdaMass); // Trick because hypernuclus mass unkown
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()) +
1195  if(getS() != 0) theBalance.energy -= getS()*(ParticleTable::effectiveNucleonMass - ParticleTable::effectiveLambdaMass); // Trick because hypernuclus mass unkown
1196  if(afterRecoil)
1197  theBalance.energy -= getKineticEnergy();
1198  theBalance.momentum -= getMomentum();
1199  }
1200 
1201  return theBalance;
1202  }
1203 
1210  }
1211 
1212  void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
1213  // Deal with the projectile remnant
1214  const G4int prA = theProjectileRemnant->getA();
1215  if(prA>=1) {
1216  // Set the mass
1218  theProjectileRemnant->setMass(aMass);
1219 
1220  // Compute the excitation energy from the invariant mass
1221  const G4double anExcitationEnergy = aMass
1223 
1224  // Set the excitation energy
1225  theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
1226 
1227  // No spin!
1229 
1230  // Set the emission time
1231  theProjectileRemnant->setEmissionTime(anEmissionTime);
1232  }
1233  }
1234 
1235 }
1236 
1237 #endif
G4double getFirstCollisionTime() const
Definition: G4INCLBook.hh:83
std::vector< G4int > getBiasCollisionVector() const
Get the vector list of biased vertices on the particle path.
ParticleList::const_iterator ParticleIter
G4bool pion(G4int ityp)
void setZ(G4double az)
Set the z coordinate.
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
virtual G4double getTableMass() const
Get the tabulated particle mass.
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
Int_t nDecayAvatars
Number of decay avatars.
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
Float_t JRem[maxSizeRemnants]
Remnant spin [ ].
void setProtonSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
Int_t nBlockedDecays
Number of decays blocked by Pauli or CDPP.
G4double getWidth(const ParticleType t)
Get particle width (in s)
ProjectileRemnant * theProjectileRemnant
Pointer to the quasi-projectile.
G4int getEmittedClusters() const
Definition: G4INCLBook.hh:106
G4double computeExcitationEnergy() const
Compute the current excitation energy.
G4int emitInsideLambda()
Force emission of all Lambda (desexitation code with strangeness not implanted yet) ...
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
ParticleList const & getParticles() const
Definition: G4INCLStore.hh:253
void setNeutronSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
void setTableMass()
Set the mass of the Particle to its table mass.
ParticleSampler * theParticleSampler
const G4double effectiveNucleonMass
const G4double effectiveLambdaMass
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
Float_t firstCollisionSpectatorPosition
Position of the spectator on the first collision (fm)
Short_t St
Strangeness number of the target nucleus.
G4double getKineticEnergy() const
Get the particle kinetic energy.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
G4double getFirstCollisionXSec() const
Definition: G4INCLBook.hh:86
Float_t emissionTime[maxSizeParticles]
Emission time [fm/c].
void setMass(G4double mass)
const G4INCL::ThreeVector & getMomentum() const
G4int getBlockedDecays() const
Definition: G4INCLBook.hh:103
G4int getS() const
Returns the strangeness number.
G4INCL::ThreeVector theMomentum
Int_t nReflectionAvatars
Number of reflection avatars.
G4bool decayInsideStrangeParticles()
Force the transformation of strange particles into a Lambda;.
G4double initialEnergy
NuclearPotential::INuclearPotential const * thePotential
Pointer to the NuclearPotential object.
const char * p
Definition: xmltok.h:285
FinalStateValidity getValidity() const
G4bool getFirstCollisionIsElastic() const
Definition: G4INCLBook.hh:95
G4int getEnergyViolationInteraction() const
Definition: G4INCLBook.hh:107
G4bool getPionPotential() const
Do we want the pion potential?
void emitInsidePions()
Force emission of all pions inside the nucleus.
void particleHasBeenUpdated(Particle *const)
Notify the Store about a particle update.
Definition: G4INCLStore.cc:127
std::vector< std::string > history
History of the particle.
G4double initialInternalEnergy
Book & getBook()
Definition: G4INCLStore.hh:259
Short_t nParticles
Number of particles in the final state.
Short_t Z[maxSizeParticles]
Particle charge number.
G4bool nucleon(G4int ityp)
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
G4double getExcitationEnergy() const
Get the excitation energy of the cluster.
void propagateParticles(G4double step)
Short_t nRemnants
Number of remnants.
G4double getMass() const
Get the cached particle mass.
void applyFinalState(FinalState *)
virtual G4double getTableMass() const
Get the real particle mass.
NuclearDensity const * theDensity
Pointer to the NuclearDensity object.
NuclearDensity const * createDensity(const G4int A, const G4int Z)
void deleteProjectileRemnant()
Delete the projectile remnant.
void setX(G4double ax)
Set the x coordinate.
Double_t beta
ParticleList particles
void setBiasCollisionVector(std::vector< G4int > BiasCollisionVector)
Set the vector list of biased vertices on the particle path.
G4double theExcitationEnergy
void computeOneNucleonRecoilKinematics()
Compute the recoil kinematics for a 1-nucleon remnant.
FinalState * getFinalState()
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
Float_t phi[maxSizeParticles]
Particle momentum azimuthal angle [radians].
void updatePotentialEnergy(Particle *p) const
Update the particle potential energy.
Short_t origin[maxSizeParticles]
Origin of the particle.
G4double getFirstCollisionSpectatorMomentum() const
Definition: G4INCLBook.hh:92
void add(Particle *p)
Definition: G4INCLStore.cc:58
G4double computeTotalEnergy() const
Compute the current total energy.
#define INCL_DEBUG(x)
void particleHasBeenEjected(Particle *const)
Definition: G4INCLStore.cc:175
void initializeParticles()
G4double getZ() const
std::string print() const
static G4double getTotalBias()
General bias vector function.
G4double phi() const
void fillEventInfo(EventInfo *eventInfo)
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
INuclearPotential const * createPotential(const PotentialType type, const G4int theA, const G4int theZ, const G4bool pionPotential)
Create an INuclearPotential object.
Bool_t forcedCompoundNucleus
True if the event is a forced CN.
Int_t nCollisionAvatars
Number of collision avatars.
Float_t phiRem[maxSizeRemnants]
Remnant momentum azimuthal angle [radians].
ParticleList const & getParticles() const
double G4double
Definition: G4Types.hh:76
G4bool hasPionPotential() const
Do we have a pion potential?
bool G4bool
Definition: G4Types.hh:79
Float_t firstCollisionTime
Time of the first collision [fm/c].
Double_t scale
G4double getEmissionTime()
Bool_t pionAbsorption
True if the event is a pion absorption.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
G4double shoot()
Definition: G4INCLRandom.cc:93
G4double getPotentialEnergy() const
Get the particle potential energy.
G4bool isDelta() const
Is it a Delta?
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
void setPotential(NuclearPotential::INuclearPotential const *const p)
Setter for thePotential.
Short_t Sp
Strangeness number of the projectile nucleus.
G4double getInvariantMass() const
Get the the particle invariant mass.
G4double getSeparationEnergy(const Particle *const p) const
Return the separation energy for a particle.
void removeScheduledAvatars()
Remove avatars that have been scheduled.
Definition: G4INCLStore.cc:134
#define INCL_ERROR(x)
Short_t ZRem[maxSizeRemnants]
Remnant charge number.
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
Float_t firstCollisionXSec
Cross section of the first collision (mb)
Short_t Zt
Charge number of the target nucleus.
Short_t ARem[maxSizeRemnants]
Remnant mass number.
Float_t thetaRem[maxSizeRemnants]
Remnant momentum polar angle [radians].
G4INCL::ThreeVector thePosition
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
TMarker * pt
Definition: egs.C:25
void setPotentialEnergy(G4double v)
Set the particle potential energy.
ParticleList const & getEnteringParticles() const
std::string print() const
ParticleList const & getModifiedParticles() const
G4int getZ() const
Returns the charge number.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
PotentialType getPotentialType() const
Get the type of the potential for nucleons.
G4double getX() const
Short_t S[maxSizeParticles]
Particle strangeness number.
void particleHasBeenDestroyed(Particle *const)
Definition: G4INCLStore.cc:181
G4int getA() const
Returns the baryon number.
G4bool decayMe()
Force the phase-space decay of the Nucleus.
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
void emitInsideStrangeParticles()
Force emission of all strange particles inside the nucleus.
Simple class implementing De Jong&#39;s spin model for nucleus-nucleus collisions.
G4bool isEventTransparent() const
Is the event transparent?
ThreeVector initialCenterOfMass
G4double getFirstCollisionSpectatorPosition() const
Definition: G4INCLBook.hh:89
Nucleus(G4int mass, G4int charge, G4int strangess, Config const *const conf, const G4double universeRadius=-1.)
G4int getAcceptedDecays() const
Definition: G4INCLBook.hh:102
G4bool emitInsideKaon()
Force emission of all Kaon inside the nucleus.
ParticleList const & getDestroyedParticles() const
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
Bool_t firstCollisionIsElastic
True if the first collision was elastic.
Struct for conservation laws.
ParticleList const & getOutgoingParticles() const
G4double getTotalEnergyBeforeInteraction() const
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
static G4ThreadLocal Int_t eventNumber
Number of the event.
G4double theUniverseRadius
The radius of the universe.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4double toDegrees(G4double radians)
G4int getAvatars(AvatarType type) const
Definition: G4INCLBook.hh:104
virtual ~Nucleus()
ThreeVector incomingAngularMomentum
ThreeVector const & getSpin() const
Get the spin of the nucleus.
G4int getAcceptedCollisions() const
Definition: G4INCLBook.hh:100
#define INCL_WARN(x)
int G4int
Definition: G4Types.hh:78
G4bool decayOutgoingSigmaZero(G4double timeThreshold)
Force the decay of outgoing Neutral Sigma.
ifstream in
Definition: comparison.C:7
Int_t PDGCode[maxSizeParticles]
PDG numbering of the particles.
ThreeVector incomingMomentum
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
Short_t nCascadeParticles
Number of cascade particles.
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
G4int getBlockedCollisions() const
Definition: G4INCLBook.hh:101
#define INCL_DATABLOCK(x)
void boost(const ThreeVector &aBoostVector)
Int_t nCollisions
Number of accepted two-body collisions.
G4double mag2() const
Short_t Ap
Mass number of the projectile nucleus.
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
void setEnergy(G4double energy)
G4double getInitialEnergy() const
Get the initial energy.
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
Definition: G4INCLStore.hh:190
static constexpr double cm
Definition: G4SIunits.hh:119
G4double mag() const
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
virtual void initializeParticles()
Initialise the NuclearDensity pointer and sample the particles.
G4double theta() const
G4double computeSeparationEnergyBalance() const
Outgoing - incoming separation energies.
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
Float_t ParticleBias[maxSizeParticles]
Particle weight due to the bias.
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
ThreeVector theSpin
Short_t A[maxSizeParticles]
Particle mass number.
Short_t Zp
Charge number of the projectile nucleus.
G4bool decayOutgoingNeutralKaon()
Force the transformation of outgoing Neutral Kaon into propation eigenstate.
Static class for carrying out cluster decays.
ThreeVector computeCenterOfMass() const
Compute the current center-of-mass position.
Int_t projectileType
Projectile particle type.
std::string print()
void setEmissionTime(G4double t)
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
ParticleList const & getOutgoingParticles() const
Definition: G4INCLStore.hh:223
Int_t nBlockedCollisions
Number of two-body collisions blocked by Pauli or CDPP.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
void insertParticle(Particle *p)
Insert a new particle (e.g. a projectile) in the nucleus.
Short_t SRem[maxSizeRemnants]
Remnant strangeness number.
ParticleList const & getCreatedParticles() const
G4double getCurrentTime() const
Definition: G4INCLBook.hh:98
Short_t At
Mass number of the target nucleus.
G4double getY() const
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
Store * getStore() const
Float_t theta[maxSizeParticles]
Particle momentum polar angle [radians].
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
G4bool decayOutgoingPionResonances(G4double timeThreshold)
Force the decay of outgoing PionResonances (eta/omega).
const G4double hc
[MeV*fm]
void setY(G4double ay)
Set the y coordinate.
Int_t nDecays
Number of accepted Delta decays.
Float_t firstCollisionSpectatorMomentum
Momentum of the spectator on the first collision (fm)