Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLNucleus.cc
Go to the documentation of this file.
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"
55 #include "G4INCLCluster.hh"
56 #include "G4INCLClusterDecay.hh"
57 #include "G4INCLDeJongSpin.hh"
58 #include "G4INCLParticleSpecies.hh"
59 #include <iterator>
60 #include <cstdlib>
61 #include <sstream>
62 // #include <cassert>
63 
64 namespace G4INCL {
65 
66  Nucleus::Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius)
67  : Cluster(charge,mass,true),
68  theInitialZ(charge), theInitialA(mass),
69  theNpInitial(0), theNnInitial(0),
70  initialInternalEnergy(0.),
71  incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
72  initialCenterOfMass(0.,0.,0.),
73  remnant(true),
74  initialEnergy(0.),
75  tryCN(false),
76  theUniverseRadius(universeRadius),
77  isNucleusNucleus(false),
78  theProjectileRemnant(NULL),
79  theDensity(NULL),
80  thePotential(NULL)
81  {
82  PotentialType potentialType;
83  G4bool pionPotential;
84  if(conf) {
85  potentialType = conf->getPotentialType();
86  pionPotential = conf->getPionPotential();
87  } else { // By default we don't use energy dependent
88  // potential. This is convenient for some tests.
89  potentialType = IsospinPotential;
90  pionPotential = true;
91  }
92 
93  thePotential = NuclearPotential::createPotential(potentialType, theA, theZ, pionPotential);
94 
97 
99 
100  theParticleSampler->setPotential(thePotential);
101  theParticleSampler->setDensity(theDensity);
102 
103  if(theUniverseRadius<0)
104  theUniverseRadius = theDensity->getMaximumRadius();
105  theStore = new Store(conf);
106  }
107 
109  delete theStore;
111  /* We don't delete the potential and the density here any more -- Factories
112  * are caching them
113  delete thePotential;
114  delete theDensity;*/
115  }
116 
118  // Reset the variables connected with the projectile remnant
119  delete theProjectileRemnant;
120  theProjectileRemnant = NULL;
121 
123  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
125  }
126  theStore->add(particles);
127  particles.clear();
128  initialInternalEnergy = computeTotalEnergy();
129  initialCenterOfMass = thePosition;
130  }
131 
133  if(!finalstate) // do nothing if no final state was returned
134  return;
135 
136  G4double totalEnergy = 0.0;
137 
138  FinalStateValidity const validity = finalstate->getValidity();
139  if(validity == ValidFS) {
140 
141  ParticleList const &created = finalstate->getCreatedParticles();
142  for(ParticleIter iter=created.begin(), e=created.end(); iter!=e; ++iter) {
143  theStore->add((*iter));
144  if(!(*iter)->isOutOfWell()) {
145  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
146  }
147  }
148 
149  ParticleList const &deleted = finalstate->getDestroyedParticles();
150  for(ParticleIter iter=deleted.begin(), e=deleted.end(); iter!=e; ++iter) {
151  theStore->particleHasBeenDestroyed(*iter);
152  }
153 
154  ParticleList const &modified = finalstate->getModifiedParticles();
155  for(ParticleIter iter=modified.begin(), e=modified.end(); iter!=e; ++iter) {
156  theStore->particleHasBeenUpdated(*iter);
157  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
158  }
159 
160  ParticleList const &out = finalstate->getOutgoingParticles();
161  for(ParticleIter iter=out.begin(), e=out.end(); iter!=e; ++iter) {
162  if((*iter)->isCluster()) {
163  Cluster *clusterOut = dynamic_cast<Cluster*>((*iter));
164 // assert(clusterOut);
165 #ifdef INCLXX_IN_GEANT4_MODE
166  if(!clusterOut)
167  continue;
168 #endif
169  ParticleList const &components = clusterOut->getParticles();
170  for(ParticleIter in=components.begin(), end=components.end(); in!=end; ++in)
171  theStore->particleHasBeenEjected(*in);
172  } else {
173  theStore->particleHasBeenEjected(*iter);
174  }
175  totalEnergy += (*iter)->getEnergy(); // No potential here because the particle is gone
176  theA -= (*iter)->getA();
177  theZ -= (*iter)->getZ();
178  theStore->addToOutgoing(*iter);
179  (*iter)->setEmissionTime(theStore->getBook().getCurrentTime());
180  }
181 
182  ParticleList const &entering = finalstate->getEnteringParticles();
183  for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
184  insertParticle(*iter);
185  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
186  }
187 
188  // actually perform the removal of the scheduled avatars
189  theStore->removeScheduledAvatars();
190  } else if(validity == ParticleBelowFermiFS || validity == ParticleBelowZeroFS) {
191  INCL_DEBUG("A Particle is entering below the Fermi sea:" << '\n' << finalstate->print() << '\n');
192  tryCN = true;
193  ParticleList const &entering = finalstate->getEnteringParticles();
194  for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
195  insertParticle(*iter);
196  }
197  }
198 
199  if(validity==ValidFS &&
200  std::abs(totalEnergy - finalstate->getTotalEnergyBeforeInteraction()) > 0.1) {
201  INCL_ERROR("Energy nonconservation! Energy at the beginning of the event = "
202  << finalstate->getTotalEnergyBeforeInteraction()
203  <<" and after interaction = "
204  << totalEnergy << '\n'
205  << finalstate->print());
206  }
207  }
208 
210  INCL_WARN("Useless Nucleus::propagateParticles -method called." << '\n');
211  }
212 
214  G4double totalEnergy = 0.0;
215  ParticleList const &inside = theStore->getParticles();
216  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
217  if((*p)->isNucleon()) // Ugly: we should calculate everything using total energies!
218  totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
219  else if((*p)->isResonance())
220  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::effectiveNucleonMass;
221  else
222  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
223  }
224  return totalEnergy;
225  }
226 
228  // If the remnant consists of only one nucleon, we need to apply a special
229  // procedure to put it on mass shell.
230  if(theA==1) {
231  emitInsidePions();
232  computeOneNucleonRecoilKinematics();
233  remnant=false;
234  return;
235  }
236 
237  // Compute the recoil momentum and angular momentum
238  theMomentum = incomingMomentum;
239  theSpin = incomingAngularMomentum;
240 
241  ParticleList const &outgoing = theStore->getOutgoingParticles();
242  for(ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p) {
243  theMomentum -= (*p)->getMomentum();
244  theSpin -= (*p)->getAngularMomentum();
245  }
246  if(theProjectileRemnant) {
247  theMomentum -= theProjectileRemnant->getMomentum();
248  theSpin -= theProjectileRemnant->getAngularMomentum();
249  }
250 
251  // Subtract orbital angular momentum
253  theSpin -= (thePosition-initialCenterOfMass).vector(theMomentum);
254 
257  remnant=true;
258  }
259 
261  ThreeVector cm(0.,0.,0.);
262  G4double totalMass = 0.0;
263  ParticleList const &inside = theStore->getParticles();
264  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
265  const G4double mass = (*p)->getMass();
266  cm += (*p)->getPosition() * mass;
267  totalMass += mass;
268  }
269  cm /= totalMass;
270  return cm;
271  }
272 
274  const G4double totalEnergy = computeTotalEnergy();
275  const G4double separationEnergies = computeSeparationEnergyBalance();
276 
277  return totalEnergy - initialInternalEnergy - separationEnergies;
278  }
279 
280  std::string Nucleus::print()
281  {
282  std::stringstream ss;
283  ss << "Particles in the nucleus:" << '\n'
284  << "Inside:" << '\n';
285  G4int counter = 1;
286  ParticleList const &inside = theStore->getParticles();
287  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
288  ss << "index = " << counter << '\n'
289  << (*p)->print();
290  counter++;
291  }
292  ss <<"Outgoing:" << '\n';
293  ParticleList const &outgoing = theStore->getOutgoingParticles();
294  for(ParticleIter p=outgoing.begin(), e=outgoing.end(); p!=e; ++p)
295  ss << (*p)->print();
296 
297  return ss.str();
298  }
299 
301  ParticleList const &out = theStore->getOutgoingParticles();
302  ParticleList deltas;
303  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
304  if((*i)->isDelta()) deltas.push_back((*i));
305  }
306  if(deltas.empty()) return false;
307 
308  for(ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
309  INCL_DEBUG("Decay outgoing delta particle:" << '\n'
310  << (*i)->print() << '\n');
311  const ThreeVector beta = -(*i)->boostVector();
312  const G4double deltaMass = (*i)->getMass();
313 
314  // Set the delta momentum to zero and sample the decay in the CM frame.
315  // This makes life simpler if we are using real particle masses.
316  (*i)->setMomentum(ThreeVector());
317  (*i)->setEnergy((*i)->getMass());
318 
319  // Use a DecayAvatar
320  IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
321  FinalState *fs = decay->getFinalState();
322  Particle * const pion = fs->getCreatedParticles().front();
323  Particle * const nucleon = fs->getModifiedParticles().front();
324 
325  // Adjust the decay momentum if we are using the real masses
326  const G4double decayMomentum = KinematicsUtils::momentumInCM(deltaMass,
327  nucleon->getTableMass(),
328  pion->getTableMass());
329  ThreeVector newMomentum = pion->getMomentum();
330  newMomentum *= decayMomentum / newMomentum.mag();
331 
332  pion->setTableMass();
333  pion->setMomentum(newMomentum);
334  pion->adjustEnergyFromMomentum();
335  pion->setEmissionTime(nucleon->getEmissionTime());
336  pion->boost(beta);
337 
338  nucleon->setTableMass();
339  nucleon->setMomentum(-newMomentum);
340  nucleon->adjustEnergyFromMomentum();
341  nucleon->boost(beta);
342 
343  theStore->addToOutgoing(pion);
344 
345  delete fs;
346  delete decay;
347  }
348 
349  return true;
350  }
351 
353  /* If there is a pion potential, do nothing (deltas will be counted as
354  * excitation energy).
355  * If, however, the remnant is unphysical (Z<0 or Z>A), force the deltas to
356  * decay and get rid of all the pions. In case you're wondering, you can
357  * end up with Z<0 or Z>A if the remnant contains more pi- than protons or
358  * more pi+ than neutrons, respectively.
359  */
360  const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
361  if(thePotential->hasPionPotential() && !unphysicalRemnant)
362  return false;
363 
364  // Build a list of deltas (avoid modifying the list you are iterating on).
365  ParticleList const &inside = theStore->getParticles();
366  ParticleList deltas;
367  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i)
368  if((*i)->isDelta()) deltas.push_back((*i));
369 
370  // Loop over the deltas, make them decay
371  for(ParticleIter i=deltas.begin(), e=deltas.end(); i!=e; ++i) {
372  INCL_DEBUG("Decay inside delta particle:" << '\n'
373  << (*i)->print() << '\n');
374  // Create a forced-decay avatar. Note the last boolean parameter. Note
375  // also that if the remnant is unphysical we more or less explicitly give
376  // up energy conservation and CDPP by passing a NULL pointer for the
377  // nucleus.
378  IAvatar *decay;
379  if(unphysicalRemnant) {
380  INCL_WARN("Forcing delta decay inside an unphysical remnant (A=" << theA
381  << ", Z=" << theZ << "). Might lead to energy-violation warnings."
382  << '\n');
383  decay = new DecayAvatar((*i), 0.0, NULL, true);
384  } else
385  decay = new DecayAvatar((*i), 0.0, this, true);
386  FinalState *fs = decay->getFinalState();
387 
388  // The pion can be ejected only if we managed to satisfy energy
389  // conservation and if pion emission does not lead to negative excitation
390  // energies.
391  if(fs->getValidity()==ValidFS) {
392  // Apply the final state to the nucleus
393  applyFinalState(fs);
394  }
395  delete fs;
396  delete decay;
397  }
398 
399  // If the remnant is unphysical, emit all the pions
400  if(unphysicalRemnant) {
401  INCL_DEBUG("Remnant is unphysical: Z=" << theZ << ", A=" << theA << ", emitting all the pions" << '\n');
402  emitInsidePions();
403  }
404 
405  return true;
406  }
407 
409  ParticleList const &out = theStore->getOutgoingParticles();
410  ParticleList pionResonances;
411  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
412 // if((*i)->isEta() || (*i)->isOmega()) pionResonances.push_back((*i));
413  if(((*i)->isEta() && timeThreshold > ParticleTable::getWidth(Eta)) || ((*i)->isOmega() && timeThreshold > ParticleTable::getWidth(Omega))) pionResonances.push_back((*i));
414  }
415  if(pionResonances.empty()) return false;
416 
417  for(ParticleIter i=pionResonances.begin(), e=pionResonances.end(); i!=e; ++i) {
418  INCL_DEBUG("Decay outgoing pionResonances particle:" << '\n'
419  << (*i)->print() << '\n');
420  const ThreeVector beta = -(*i)->boostVector();
421  const G4double pionResonanceMass = (*i)->getMass();
422 
423  // Set the pionResonance momentum to zero and sample the decay in the CM frame.
424  // This makes life simpler if we are using real particle masses.
425  (*i)->setMomentum(ThreeVector());
426  (*i)->setEnergy((*i)->getMass());
427 
428  // Use a DecayAvatar
429  IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
430  FinalState *fs = decay->getFinalState();
431 
432  Particle * const theModifiedParticle = fs->getModifiedParticles().front();
433  ParticleList const &created = fs->getCreatedParticles();
434  Particle * const theCreatedParticle1 = created.front();
435 
436  if (created.size() == 1) {
437 
438  // Adjust the decay momentum if we are using the real masses
439  const G4double decayMomentum = KinematicsUtils::momentumInCM(pionResonanceMass,theModifiedParticle->getTableMass(),theCreatedParticle1->getTableMass());
440  ThreeVector newMomentum = theCreatedParticle1->getMomentum();
441  newMomentum *= decayMomentum / newMomentum.mag();
442 
443  theCreatedParticle1->setTableMass();
444  theCreatedParticle1->setMomentum(newMomentum);
445  theCreatedParticle1->adjustEnergyFromMomentum();
446  //theCreatedParticle1->setEmissionTime(nucleon->getEmissionTime());
447  theCreatedParticle1->boost(beta);
448 
449  theModifiedParticle->setTableMass();
450  theModifiedParticle->setMomentum(-newMomentum);
451  theModifiedParticle->adjustEnergyFromMomentum();
452  theModifiedParticle->boost(beta);
453 
454  theStore->addToOutgoing(theCreatedParticle1);
455  }
456  else if (created.size() == 2) {
457  Particle * const theCreatedParticle2 = created.back();
458 
459  theCreatedParticle1->boost(beta);
460  theCreatedParticle2->boost(beta);
461  theModifiedParticle->boost(beta);
462 
463  theStore->addToOutgoing(theCreatedParticle1);
464  theStore->addToOutgoing(theCreatedParticle2);
465  }
466  else {
467  INCL_ERROR("Wrong number (< 2) of created particles during the decay of a pion resonance");
468  }
469  delete fs;
470  delete decay;
471  }
472 
473  return true;
474  }
475 
477  ParticleList const &out = theStore->getOutgoingParticles();
478  ParticleList clusters;
479  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
480  if((*i)->isCluster()) clusters.push_back((*i));
481  }
482  if(clusters.empty()) return false;
483 
484  for(ParticleIter i=clusters.begin(), e=clusters.end(); i!=e; ++i) {
485  Cluster *cluster = dynamic_cast<Cluster*>(*i); // Can't avoid using a cast here
486 // assert(cluster);
487 #ifdef INCLXX_IN_GEANT4_MODE
488  if(!cluster)
489  continue;
490 #endif
491  cluster->deleteParticles(); // Don't need them
492  ParticleList decayProducts = ClusterDecay::decay(cluster);
493  for(ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j)
494  theStore->addToOutgoing(*j);
495  }
496  return true;
497  }
498 
500  // Do the phase-space decay only if Z=0 or Z=A
501  if(theA<=1 || (theZ!=0 && theA!=theZ))
502  return false;
503 
504  ParticleList decayProducts = ClusterDecay::decay(this);
505  for(ParticleIter j=decayProducts.begin(), e=decayProducts.end(); j!=e; ++j)
506  theStore->addToOutgoing(*j);
507 
508  return true;
509  }
510 
512  /* Forcing emissions of all pions in the nucleus. This probably violates
513  * energy conservation (although the computation of the recoil kinematics
514  * might sweep this under the carpet).
515  */
516  INCL_WARN("Forcing emissions of all pions in the nucleus." << '\n');
517 
518  // Emit the pions with this kinetic energy
519  const G4double tinyPionEnergy = 0.1; // MeV
520 
521  // Push out the emitted pions
522  ParticleList const &inside = theStore->getParticles();
523  ParticleList toEject;
524  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
525  if((*i)->isPion()) {
526  Particle * const thePion = *i;
527  INCL_DEBUG("Forcing emission of the following particle: "
528  << thePion->print() << '\n');
529  thePion->setEmissionTime(theStore->getBook().getCurrentTime());
530  // Correction for real masses
531  const G4double theQValueCorrection = thePion->getEmissionQValueCorrection(theA,theZ);
532  const G4double kineticEnergyOutside = thePion->getKineticEnergy() - thePion->getPotentialEnergy() + theQValueCorrection;
533  thePion->setTableMass();
534  if(kineticEnergyOutside > 0.0)
535  thePion->setEnergy(thePion->getMass()+kineticEnergyOutside);
536  else
537  thePion->setEnergy(thePion->getMass()+tinyPionEnergy);
538  thePion->adjustMomentumFromEnergy();
539  thePion->setPotentialEnergy(0.);
540  theZ -= thePion->getZ();
541  toEject.push_back(thePion);
542  }
543  }
544  for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
545  theStore->particleHasBeenEjected(*i);
546  theStore->addToOutgoing(*i);
547  }
548  }
549 
551 
552  Book const &theBook = theStore->getBook();
553  const G4int nEventCollisions = theBook.getAcceptedCollisions();
554  const G4int nEventDecays = theBook.getAcceptedDecays();
555  const G4int nEventClusters = theBook.getEmittedClusters();
556  if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
557  return true;
558 
559  return false;
560 
561  }
562 
563  void Nucleus::computeOneNucleonRecoilKinematics() {
564  // We should be here only if the nucleus contains only one nucleon
565 // assert(theStore->getParticles().size()==1);
566 
567  // No excitation energy!
568  theExcitationEnergy = 0.0;
569 
570  // Move the nucleon to the outgoing list
571  Particle *remN = theStore->getParticles().front();
572  theA -= remN->getA();
573  theZ -= remN->getZ();
574  theStore->particleHasBeenEjected(remN);
575  theStore->addToOutgoing(remN);
576  remN->setEmissionTime(theStore->getBook().getCurrentTime());
577 
578  // Treat the special case of a remaining delta
579  if(remN->isDelta()) {
580  IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
581  FinalState *fs = decay->getFinalState();
582  // Eject the pion
583  ParticleList const &created = fs->getCreatedParticles();
584  for(ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
585  theStore->addToOutgoing(*j);
586  delete fs;
587  delete decay;
588  }
589 
590  // Do different things depending on how many outgoing particles we have
591  ParticleList const &outgoing = theStore->getOutgoingParticles();
592  if(outgoing.size() == 2) {
593 
594  INCL_DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << '\n');
595 
596  // Can apply exact 2-body kinematics here. Keep the CM emission angle of
597  // the first particle.
598  Particle *p1 = outgoing.front(), *p2 = outgoing.back();
599  const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
600  // Boost to the initial CM
601  p1->boost(aBoostVector);
602  const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
603  const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
604  const G4double scale = pcm/(p1->getMomentum().mag());
605  // Reset the momenta
606  p1->setMomentum(p1->getMomentum()*scale);
607  p2->setMomentum(-p1->getMomentum());
608  p1->adjustEnergyFromMomentum();
609  p2->adjustEnergyFromMomentum();
610  // Unboost
611  p1->boost(-aBoostVector);
612  p2->boost(-aBoostVector);
613 
614  } else {
615 
616  INCL_DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << '\n');
617 
618  const G4int maxIterations=8;
619  G4double totalEnergy, energyScale;
620  G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
621  ThreeVector totalMomentum, deltaP;
622  std::vector<ThreeVector> minMomenta; // use it to store the particle momenta that minimize the merit function
623 
624  // Reserve the vector size
625  minMomenta.reserve(outgoing.size());
626 
627  // Compute the initial total momentum
628  totalMomentum.setX(0.0);
629  totalMomentum.setY(0.0);
630  totalMomentum.setZ(0.0);
631  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
632  totalMomentum += (*i)->getMomentum();
633 
634  // Compute the initial total energy
635  totalEnergy = 0.0;
636  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
637  totalEnergy += (*i)->getEnergy();
638 
639  // Iterative algorithm starts here:
640  for(G4int iterations=0; iterations < maxIterations; ++iterations) {
641 
642  // Save the old merit-function values
643  oldOldOldVal = oldOldVal;
644  oldOldVal = oldVal;
645  oldVal = val;
646 
647  if(iterations%2 == 0) {
648  INCL_DEBUG("Momentum step" << '\n');
649  // Momentum step: modify all the particle momenta
650  deltaP = incomingMomentum - totalMomentum;
651  G4double pOldTot = 0.0;
652  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
653  pOldTot += (*i)->getMomentum().mag();
654  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
655  const ThreeVector mom = (*i)->getMomentum();
656  (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
657  (*i)->adjustEnergyFromMomentum();
658  }
659  } else {
660  INCL_DEBUG("Energy step" << '\n');
661  // Energy step: modify all the particle momenta
662  energyScale = initialEnergy/totalEnergy;
663  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
664  const ThreeVector mom = (*i)->getMomentum();
665  G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
666  if(pScale>0) {
667  (*i)->setEnergy((*i)->getEnergy()*energyScale);
668  (*i)->adjustMomentumFromEnergy();
669  }
670  }
671  }
672 
673  // Compute the current total momentum and energy
674  totalMomentum.setX(0.0);
675  totalMomentum.setY(0.0);
676  totalMomentum.setZ(0.0);
677  totalEnergy = 0.0;
678  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
679  totalMomentum += (*i)->getMomentum();
680  totalEnergy += (*i)->getEnergy();
681  }
682 
683  // Merit factor
684  val = std::pow(totalEnergy - initialEnergy,2) +
685  0.25*(totalMomentum - incomingMomentum).mag2();
686  INCL_DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
687 
688  // Store the minimum
689  if(val < oldVal) {
690  INCL_DEBUG("New minimum found, storing the particle momenta" << '\n');
691  minMomenta.clear();
692  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
693  minMomenta.push_back((*i)->getMomentum());
694  }
695 
696  // Stop the algorithm if the search diverges
697  if(val > oldOldVal && oldVal > oldOldOldVal) {
698  INCL_DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << '\n');
699  break;
700  }
701  }
702 
703  // We should have made at least one successful iteration here
704 // assert(minMomenta.size()==outgoing.size());
705 
706  // Apply the optimal momenta
707  INCL_DEBUG("Applying the solution" << '\n');
708  std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
709  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
710  (*i)->setMomentum(*v);
711  (*i)->adjustEnergyFromMomentum();
712  INCL_DATABLOCK((*i)->print());
713  }
714 
715  }
716 
717  }
718 
720  eventInfo->nParticles = 0;
721  G4bool isNucleonAbsorption = false;
722 
723  G4bool isPionAbsorption = false;
724  // It is possible to have pion absorption event only if the
725  // projectile is pion.
726  if(eventInfo->projectileType == PiPlus ||
727  eventInfo->projectileType == PiMinus ||
728  eventInfo->projectileType == PiZero) {
729  isPionAbsorption = true;
730  }
731 
732  // Forced CN
733  eventInfo->forcedCompoundNucleus = tryCN;
734 
735  // Outgoing particles
736  ParticleList const &outgoingParticles = getStore()->getOutgoingParticles();
737 
738  // Check if we have a nucleon absorption event: nucleon projectile
739  // and no ejected particles.
740  if(outgoingParticles.size() == 0 &&
741  (eventInfo->projectileType == Proton ||
742  eventInfo->projectileType == Neutron)) {
743  isNucleonAbsorption = true;
744  }
745 
746  // Reset the remnant counter
747  eventInfo->nRemnants = 0;
748  eventInfo->history.clear();
749 
750  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
751  // We have a pion absorption event only if the projectile is
752  // pion and there are no ejected pions.
753  if(isPionAbsorption) {
754  if((*i)->isPion()) {
755  isPionAbsorption = false;
756  }
757  }
758 
759  eventInfo->A[eventInfo->nParticles] = (*i)->getA();
760  eventInfo->Z[eventInfo->nParticles] = (*i)->getZ();
761  eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
762  eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
763  ThreeVector mom = (*i)->getMomentum();
764  eventInfo->px[eventInfo->nParticles] = mom.getX();
765  eventInfo->py[eventInfo->nParticles] = mom.getY();
766  eventInfo->pz[eventInfo->nParticles] = mom.getZ();
767  eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
768  eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
769  eventInfo->origin[eventInfo->nParticles] = -1;
770  eventInfo->history.push_back("");
771  if ((*i)->getType() != Composite) {
772  ParticleSpecies pt((*i)->getType());
773  eventInfo->PDGCode[eventInfo->nParticles] = pt.ParticleSpecies::getPDGCode();
774  }
775  else {
776  ParticleSpecies pt((*i)->getA(), (*i)->getZ());
777  eventInfo->PDGCode[eventInfo->nParticles] = pt.ParticleSpecies::getPDGCode();
778  }
779  eventInfo->nParticles++;
780  }
781  eventInfo->nucleonAbsorption = isNucleonAbsorption;
782  eventInfo->pionAbsorption = isPionAbsorption;
783  eventInfo->nCascadeParticles = eventInfo->nParticles;
784 
785  // Projectile-like remnant characteristics
786  if(theProjectileRemnant && theProjectileRemnant->getA()>0) {
787  eventInfo->ARem[eventInfo->nRemnants] = theProjectileRemnant->getA();
788  eventInfo->ZRem[eventInfo->nRemnants] = theProjectileRemnant->getZ();
789  G4double eStar = theProjectileRemnant->getExcitationEnergy();
790  if(std::abs(eStar)<1E-10)
791  eStar = 0.0; // blame rounding and set the excitation energy to zero
792  eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
793  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
794  INCL_WARN("Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << '\n');
795  }
796  const ThreeVector &spin = theProjectileRemnant->getSpin();
797  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
798  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
799  } else { // odd-A nucleus
800  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
801  }
802  eventInfo->EKinRem[eventInfo->nRemnants] = theProjectileRemnant->getKineticEnergy();
803  const ThreeVector &mom = theProjectileRemnant->getMomentum();
804  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
805  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
806  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
807  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
808  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
809  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
810  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
811  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
812  eventInfo->nRemnants++;
813  }
814 
815  // Target-like remnant characteristics
816  if(hasRemnant()) {
817  eventInfo->ARem[eventInfo->nRemnants] = getA();
818  eventInfo->ZRem[eventInfo->nRemnants] = getZ();
819  eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
820  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
821  INCL_WARN("Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << '\n');
822  }
823  const ThreeVector &spin = getSpin();
824  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
825  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
826  } else { // odd-A nucleus
827  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
828  }
829  eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
830  const ThreeVector &mom = getMomentum();
831  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
832  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
833  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
834  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
835  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
836  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
837  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
838  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
839  eventInfo->nRemnants++;
840  }
841 
842  // Global counters, flags, etc.
843  Book const &theBook = theStore->getBook();
844  eventInfo->nCollisions = theBook.getAcceptedCollisions();
845  eventInfo->nBlockedCollisions = theBook.getBlockedCollisions();
846  eventInfo->nDecays = theBook.getAcceptedDecays();
847  eventInfo->nBlockedDecays = theBook.getBlockedDecays();
848  eventInfo->firstCollisionTime = theBook.getFirstCollisionTime();
849  eventInfo->firstCollisionXSec = theBook.getFirstCollisionXSec();
853  eventInfo->nReflectionAvatars = theBook.getAvatars(SurfaceAvatarType);
854  eventInfo->nCollisionAvatars = theBook.getAvatars(CollisionAvatarType);
855  eventInfo->nDecayAvatars = theBook.getAvatars(DecayAvatarType);
857  }
858 
859  Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
860  ConservationBalance theBalance;
861  // Initialise balance variables with the incoming values
862  theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
863  theBalance.A = theEventInfo.Ap + theEventInfo.At;
864 
865  theBalance.energy = getInitialEnergy();
866  theBalance.momentum = getIncomingMomentum();
867 
868  // Process outgoing particles
869  ParticleList const &outgoingParticles = theStore->getOutgoingParticles();
870  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
871  theBalance.Z -= (*i)->getZ();
872  theBalance.A -= (*i)->getA();
873  // For outgoing clusters, the total energy automatically includes the
874  // excitation energy
875  theBalance.energy -= (*i)->getEnergy(); // Note that outgoing particles should have the real mass
876  theBalance.momentum -= (*i)->getMomentum();
877  }
878 
879  // Projectile-like remnant contribution, if present
880  if(theProjectileRemnant && theProjectileRemnant->getA()>0) {
881  theBalance.Z -= theProjectileRemnant->getZ();
882  theBalance.A -= theProjectileRemnant->getA();
883  theBalance.energy -= ParticleTable::getTableMass(theProjectileRemnant->getA(),theProjectileRemnant->getZ()) +
884  theProjectileRemnant->getExcitationEnergy();
885  theBalance.energy -= theProjectileRemnant->getKineticEnergy();
886  theBalance.momentum -= theProjectileRemnant->getMomentum();
887  }
888 
889  // Target-like remnant contribution, if present
890  if(hasRemnant()) {
891  theBalance.Z -= getZ();
892  theBalance.A -= getA();
893  theBalance.energy -= ParticleTable::getTableMass(getA(),getZ()) +
895  if(afterRecoil)
896  theBalance.energy -= getKineticEnergy();
897  theBalance.momentum -= getMomentum();
898  }
899 
900  return theBalance;
901  }
902 
904  setEnergy(initialEnergy);
905  setMomentum(incomingMomentum);
906  setSpin(incomingAngularMomentum);
909  }
910 
911  void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
912  // Deal with the projectile remnant
913  const G4int prA = theProjectileRemnant->getA();
914  if(prA>=1) {
915  // Set the mass
916  const G4double aMass = theProjectileRemnant->getInvariantMass();
917  theProjectileRemnant->setMass(aMass);
918 
919  // Compute the excitation energy from the invariant mass
920  const G4double anExcitationEnergy = aMass
921  - ParticleTable::getTableMass(prA, theProjectileRemnant->getZ());
922 
923  // Set the excitation energy
924  theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
925 
926  // No spin!
927  theProjectileRemnant->setSpin(ThreeVector());
928 
929  // Set the emission time
930  theProjectileRemnant->setEmissionTime(anEmissionTime);
931  }
932  }
933 
934 }
935 
936 #endif
Short_t Zp
Charge number of the projectile nucleus.
G4int getA() const
Returns the baryon number.
ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const
Compute charge, mass, energy and momentum balance.
void initializeParticles()
Short_t nRemnants
Number of remnants.
Bool_t nucleonAbsorption
True if the event is a nucleon absorption.
void setMass(G4double mass)
void updatePotentialEnergy(Particle *p) const
Update the particle potential energy.
G4double getFirstCollisionSpectatorMomentum() const
Definition: G4INCLBook.hh:92
FinalStateValidity getValidity() const
virtual G4double getTableMass() const
Get the real particle mass.
G4double getFirstCollisionSpectatorPosition() const
Definition: G4INCLBook.hh:89
Simple class implementing De Jong&#39;s spin model for nucleus-nucleus collisions.
G4double getMass() const
Get the cached particle mass.
Int_t PDGCode[maxSizeParticles]
PDG numbering of the particles.
Int_t nEnergyViolationInteraction
Number of attempted collisions/decays for which the energy-conservation algorithm failed to find a so...
ParticleList const & getParticles() const
Definition: G4INCLStore.hh:253
Float_t firstCollisionTime
Time of the first collision [fm/c].
std::string print()
G4int getEmittedClusters() const
Definition: G4INCLBook.hh:106
G4bool pion(G4int ityp)
G4int getAcceptedCollisions() const
Definition: G4INCLBook.hh:100
Float_t emissionTime[maxSizeParticles]
Emission time [fm/c].
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
virtual ~Nucleus()
ParticleList const & getParticles() const
G4double getEmissionTime()
Int_t nCollisions
Number of accepted two-body collisions.
Float_t py[maxSizeParticles]
Particle momentum, y component [MeV/c].
G4bool decayOutgoingClusters()
Force the decay of unstable outgoing clusters.
Float_t firstCollisionXSec
Cross section of the first collision (mb)
ParticleList const & getModifiedParticles() const
void setEmissionTime(G4double t)
Static class for carrying out cluster decays.
const char * p
Definition: xmltok.h:285
void applyFinalState(FinalState *)
#define INCL_ERROR(x)
const G4double hc
[MeV*fm]
void boost(const ThreeVector &aBoostVector)
const G4INCL::ThreeVector & getMomentum() const
G4double adjustEnergyFromMomentum()
Recompute the energy to match the momentum.
Bool_t forcedCompoundNucleus
True if the event is a forced CN.
Store * getStore() const
Float_t jzRem[maxSizeRemnants]
Remnant angular momentum, z component [ ].
G4bool isDelta() const
Is it a Delta?
std::string print() const
void add(Particle *p)
Definition: G4INCLStore.cc:58
Int_t nBlockedCollisions
Number of two-body collisions blocked by Pauli or CDPP.
Struct for conservation laws.
#define INCL_WARN(x)
Bool_t pionAbsorption
True if the event is a pion absorption.
G4double toDegrees(G4double radians)
Float_t EKin[maxSizeParticles]
Particle kinetic energy [MeV].
G4bool hasRemnant() const
Does the nucleus give a cascade remnant?
G4double momentumInCM(Particle const *const p1, Particle const *const p2)
gives the momentum in the CM frame of two particles.
Float_t JRem[maxSizeRemnants]
Remnant spin [ ].
Int_t nReflectionAvatars
Number of reflection avatars.
Float_t EKinRem[maxSizeRemnants]
Remnant kinetic energy [MeV].
int G4int
Definition: G4Types.hh:78
G4int getBlockedCollisions() const
Definition: G4INCLBook.hh:101
G4double mag2() const
void setProtonSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
G4double getFirstCollisionTime() const
Definition: G4INCLBook.hh:83
Float_t pyRem[maxSizeRemnants]
Remnant momentum, y component [MeV/c].
G4double getPotentialEnergy() const
Get the particle potential energy.
Int_t nCollisionAvatars
Number of collision avatars.
Float_t pz[maxSizeParticles]
Particle momentum, z component [MeV/c].
G4double getInvariantMass() const
Get the the particle invariant mass.
Short_t At
Mass number of the target nucleus.
void setEnergy(G4double energy)
Short_t ZRem[maxSizeRemnants]
Remnant charge number.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
ParticleList const & getOutgoingParticles() const
Short_t ARem[maxSizeRemnants]
Remnant mass number.
ParticleList const & getCreatedParticles() const
Short_t Zt
Charge number of the target nucleus.
G4int getAvatars(AvatarType type) const
Definition: G4INCLBook.hh:104
void removeScheduledAvatars()
Remove avatars that have been scheduled.
Definition: G4INCLStore.cc:134
Float_t theta[maxSizeParticles]
Particle momentum polar angle [radians].
void useFusionKinematics()
Adjust the kinematics for complete-fusion events.
G4double getInitialEnergy() const
Get the initial energy.
G4double getSeparationEnergy(const Particle *const p) const
Return the separation energy for a particle.
ParticleList decay(Cluster *const c)
Carries out a cluster decay.
G4bool decayOutgoingDeltas()
Force the decay of outgoing deltas.
ParticleList const & getOutgoingParticles() const
Definition: G4INCLStore.hh:223
Short_t Z[maxSizeParticles]
Particle charge number.
ThreeVector theSpin
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
Book & getBook()
Definition: G4INCLStore.hh:259
Short_t nParticles
Number of particles in the final state.
G4bool nucleon(G4int ityp)
G4int getBlockedDecays() const
Definition: G4INCLBook.hh:103
void particleHasBeenDestroyed(Particle *const)
Definition: G4INCLStore.cc:181
Float_t jxRem[maxSizeRemnants]
Remnant angular momentum, x component [ ].
bool G4bool
Definition: G4Types.hh:79
virtual void initializeParticles()
Initialise the NuclearDensity pointer and sample the particles.
G4double phi() const
void particleHasBeenEjected(Particle *const)
Definition: G4INCLStore.cc:175
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
static constexpr double cm
Definition: G4SIunits.hh:119
ParticleList particles
G4int getZ() const
Returns the charge number.
G4bool getPionPotential() const
Do we want the pion potential?
void setPotentialEnergy(G4double v)
Set the particle potential energy.
Float_t EStarRem[maxSizeRemnants]
Remnant excitation energy [MeV].
G4double computeSeparationEnergyBalance() const
Outgoing - incoming separation energies.
G4double getExcitationEnergy() const
Get the excitation energy of the nucleus.
NuclearDensity const * createDensity(const G4int A, const G4int Z)
G4double computeExcitationEnergy() const
Compute the current excitation energy.
Int_t nDecayAvatars
Number of decay avatars.
void deleteProjectileRemnant()
Delete the projectile remnant.
G4INCL::ThreeVector thePosition
G4bool decayMe()
Force the phase-space decay of the Nucleus.
G4int getAcceptedDecays() const
Definition: G4INCLBook.hh:102
tuple v
Definition: test.py:18
Float_t jyRem[maxSizeRemnants]
Remnant angular momentum, y component [ ].
ParticleSampler * theParticleSampler
Float_t firstCollisionSpectatorPosition
Position of the spectator on the first collision (fm)
Float_t phiRem[maxSizeRemnants]
Remnant momentum azimuthal angle [radians].
G4bool getFirstCollisionIsElastic() const
Definition: G4INCLBook.hh:95
Int_t nBlockedDecays
Number of decays blocked by Pauli or CDPP.
G4bool hasPionPotential() const
Do we have a pion potential?
ThreeVector computeCenterOfMass() const
Compute the current center-of-mass position.
G4double getTotalEnergyBeforeInteraction() const
G4double getFirstCollisionXSec() const
Definition: G4INCLBook.hh:86
Float_t pxRem[maxSizeRemnants]
Remnant momentum, x component [MeV/c].
Short_t A[maxSizeParticles]
Particle mass number.
Int_t nDecays
Number of accepted Delta decays.
Bool_t firstCollisionIsElastic
True if the first collision was elastic.
G4double theExcitationEnergy
void particleHasBeenUpdated(Particle *const)
Notify the Store about a particle update.
Definition: G4INCLStore.cc:127
G4double getWidth(const ParticleType t)
Get particle width (in s)
void propagateParticles(G4double step)
G4bool decayInsideDeltas()
Force the decay of deltas inside the nucleus.
void finalizeProjectileRemnant(const G4double emissionTime)
Finalise the projectile remnant.
PotentialType getPotentialType() const
Get the type of the potential for nucleons.
std::string print() const
void emitInsidePions()
Force emission of all pions inside the nucleus.
G4ThreadLocal NuclearMassFn getTableMass
Static pointer to the mass function for nuclei.
void setTableMass()
Set the mass of the Particle to its table mass.
ParticleList const & getEnteringParticles() const
ParticleList const & getDestroyedParticles() const
void setPotential(NuclearPotential::INuclearPotential const *const p)
Setter for thePotential.
Float_t thetaRem[maxSizeRemnants]
Remnant momentum polar angle [radians].
G4bool isEventTransparent() const
Is the event transparent?
FinalState * getFinalState()
void fillEventInfo(EventInfo *eventInfo)
G4double computeTotalEnergy() const
Compute the current total energy.
G4INCL::ThreeVector theMomentum
ThreeVector const & getSpin() const
Get the spin of the nucleus.
virtual G4double getTableMass() const
Get the tabulated particle mass.
G4double getX() const
G4double getExcitationEnergy() const
Get the excitation energy of the cluster.
G4double getKineticEnergy() const
Get the particle kinetic energy.
void insertParticle(Particle *p)
Insert a new particle (e.g. a projectile) in the nucleus.
G4double theta() const
Float_t phi[maxSizeParticles]
Particle momentum azimuthal angle [radians].
G4double mag() const
double G4double
Definition: G4Types.hh:76
G4double getCurrentTime() const
Definition: G4INCLBook.hh:98
Nucleus(G4int mass, G4int charge, Config const *const conf, const G4double universeRadius=-1.)
G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const
Computes correction on the emission Q-value.
#define INCL_DEBUG(x)
void computeRecoilKinematics()
Compute the recoil momentum and spin of the nucleus.
INuclearPotential const * createPotential(const PotentialType type, const G4int theA, const G4int theZ, const G4bool pionPotential)
Create an INuclearPotential object.
Float_t firstCollisionSpectatorMomentum
Momentum of the spectator on the first collision (fm)
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
void addToOutgoing(Particle *p)
add the particle to the outgoing particle list.
Definition: G4INCLStore.hh:190
G4int getEnergyViolationInteraction() const
Definition: G4INCLBook.hh:107
Float_t px[maxSizeParticles]
Particle momentum, x component [MeV/c].
Short_t Ap
Mass number of the projectile nucleus.
Short_t nCascadeParticles
Number of cascade particles.
G4bool decayOutgoingPionResonances(G4double timeThreshold)
Force the decay of outgoing PionResonances (eta/omega).
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
const G4double effectiveNucleonMass
Short_t origin[maxSizeParticles]
Origin of the particle.
#define INCL_DATABLOCK(x)
void setNeutronSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
ParticleList::const_iterator ParticleIter
G4double getZ() const
Int_t projectileType
Projectile particle type.
virtual void setMomentum(const G4INCL::ThreeVector &momentum)
G4double getY() const
std::vector< std::string > history
History of the particle.