Geant4_10
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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Davide Mancusi, CEA
29 // Alain Boudard, CEA
30 // Sylvie Leray, CEA
31 // Joseph Cugnon, University of Liege
32 //
33 #define INCLXX_IN_GEANT4_MODE 1
34 
35 #include "globals.hh"
36 
37 /*
38  * G4INCLNucleus.cc
39  *
40  * \date Jun 5, 2009
41  * \author Pekka Kaitaniemi
42  */
43 
44 #ifndef G4INCLNucleus_hh
45 #define G4INCLNucleus_hh 1
46 
47 #include "G4INCLGlobals.hh"
48 #include "G4INCLLogger.hh"
49 #include "G4INCLParticle.hh"
50 #include "G4INCLIAvatar.hh"
51 #include "G4INCLNucleus.hh"
52 #include "G4INCLKinematicsUtils.hh"
53 #include "G4INCLDecayAvatar.hh"
54 #include "G4INCLCluster.hh"
55 #include "G4INCLClusterDecay.hh"
56 #include "G4INCLDeJongSpin.hh"
57 #include <iterator>
58 #include <cstdlib>
59 #include <sstream>
60 // #include <cassert>
61 
62 namespace G4INCL {
63 
64  Nucleus::Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius)
65  : Cluster(charge,mass,true),
66  theInitialZ(charge), theInitialA(mass),
67  theNpInitial(0), theNnInitial(0),
68  initialInternalEnergy(0.),
69  incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
70  initialCenterOfMass(0.,0.,0.),
71  remnant(true),
72  blockedDelta(NULL),
73  initialEnergy(0.),
74  tryCN(false),
75  projectileZ(0),
76  projectileA(0),
77  theUniverseRadius(universeRadius),
78  isNucleusNucleus(false),
79  theProjectileRemnant(NULL),
80  theDensity(NULL),
81  thePotential(NULL)
82  {
83  PotentialType potentialType;
84  G4bool pionPotential;
85  if(conf) {
86  potentialType = conf->getPotentialType();
87  pionPotential = conf->getPionPotential();
88  } else { // By default we don't use energy dependent
89  // potential. This is convenient for some tests.
90  potentialType = IsospinPotential;
91  pionPotential = true;
92  }
93 
94  thePotential = NuclearPotential::createPotential(potentialType, theA, theZ, pionPotential);
95 
98 
100 
101  theParticleSampler->setPotential(thePotential);
102  theParticleSampler->setDensity(theDensity);
103 
104  if(theUniverseRadius<0)
105  theUniverseRadius = theDensity->getMaximumRadius();
106  theStore = new Store(conf);
107  toBeUpdated.clear();
108  }
109 
111  delete theStore;
113  /* We don't delete the potential and the density here any more -- Factories
114  * are caching them
115  delete thePotential;
116  delete theDensity;*/
117  }
118 
120  // Reset the variables connected with the projectile remnant
121  delete theProjectileRemnant;
122  theProjectileRemnant = NULL;
123 
125  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
127  theStore->add(*i);
128  }
129  particles.clear();
130  initialInternalEnergy = computeTotalEnergy();
131  initialCenterOfMass = thePosition;
132  }
133 
135  justCreated.clear();
136  toBeUpdated.clear(); // Clear the list of particles to be updated by the propagation model.
137  blockedDelta = NULL;
138 
139  if(!finalstate) // do nothing if no final state was returned
140  return;
141 
142  G4double totalEnergy = 0.0;
143 
144  FinalStateValidity const validity = finalstate->getValidity();
145  if(validity == ValidFS) {
146 
147  ParticleList const &created = finalstate->getCreatedParticles();
148  for(ParticleIter iter=created.begin(), e=created.end(); iter!=e; ++iter) {
149  theStore->add((*iter));
150  if(!(*iter)->isOutOfWell()) {
151  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
152  justCreated.push_back((*iter)); // New particle, so we must create avatars for it
153  }
154  }
155 
156  ParticleList const &deleted = finalstate->getDestroyedParticles();
157  for(ParticleIter iter=deleted.begin(), e=deleted.end(); iter!=e; ++iter) {
158  theStore->particleHasBeenDestroyed(*iter);
159  }
160 
161  ParticleList const &modified = finalstate->getModifiedParticles();
162  for(ParticleIter iter=modified.begin(), e=modified.end(); iter!=e; ++iter) {
163  theStore->particleHasBeenUpdated(*iter);
164  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
165  toBeUpdated.push_back((*iter)); // Particle is modified so we have to create new avatars for it.
166  }
167 
168  ParticleList const &out = finalstate->getOutgoingParticles();
169  for(ParticleIter iter=out.begin(), e=out.end(); iter!=e; ++iter) {
170  if((*iter)->isCluster()) {
171  Cluster *clusterOut = dynamic_cast<Cluster*>((*iter));
172 // assert(clusterOut);
173 #ifdef INCLXX_IN_GEANT4_MODE
174  if(!clusterOut)
175  continue;
176 #endif
177  ParticleList const &components = clusterOut->getParticles();
178  for(ParticleIter in=components.begin(), end=components.end(); in!=end; ++in)
179  theStore->particleHasBeenEjected(*in);
180  } else {
181  theStore->particleHasBeenEjected(*iter);
182  }
183  totalEnergy += (*iter)->getEnergy(); // No potential here because the particle is gone
184  theA -= (*iter)->getA();
185  theZ -= (*iter)->getZ();
186  theStore->addToOutgoing(*iter);
187  (*iter)->setEmissionTime(theStore->getBook().getCurrentTime());
188  }
189 
190  ParticleList const &entering = finalstate->getEnteringParticles();
191  for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
192  insertParticle(*iter);
193  totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
194  toBeUpdated.push_back((*iter)); // Particle is modified so we have to create new avatars for it.
195  }
196  } else if(validity == PauliBlockedFS) {
197  blockedDelta = finalstate->getBlockedDelta();
198  } else if(validity == ParticleBelowFermiFS || validity == ParticleBelowZeroFS) {
199  INCL_DEBUG("A Particle is entering below the Fermi sea:" << std::endl << finalstate->print() << std::endl);
200  tryCN = true;
201  ParticleList const &entering = finalstate->getEnteringParticles();
202  for(ParticleIter iter=entering.begin(), e=entering.end(); iter!=e; ++iter) {
203  insertParticle(*iter);
204  }
205  }
206 
207  if(validity==ValidFS &&
208  std::abs(totalEnergy - finalstate->getTotalEnergyBeforeInteraction()) > 0.1) {
209  INCL_ERROR("Energy nonconservation! Energy at the beginning of the event = "
210  << finalstate->getTotalEnergyBeforeInteraction()
211  <<" and after interaction = "
212  << totalEnergy << std::endl
213  << finalstate->print());
214  }
215  }
216 
218  INCL_WARN("Useless Nucleus::propagateParticles -method called." << std::endl);
219  }
220 
222  G4double totalEnergy = 0.0;
223  ParticleList const &inside = theStore->getParticles();
224  for(ParticleIter p=inside.begin(), e=inside.end(); p!=e; ++p) {
225  if((*p)->isNucleon()) // Ugly: we should calculate everything using total energies!
226  totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
227  else if((*p)->isResonance())
228  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::effectiveNucleonMass;
229  else
230  totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
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();
240  computeOneNucleonRecoilKinematics();
241  remnant=false;
242  return;
243  }
244 
245  // Compute the recoil momentum and angular momentum
246  theMomentum = incomingMomentum;
247  theSpin = incomingAngularMomentum;
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  }
254  if(theProjectileRemnant) {
255  theMomentum -= theProjectileRemnant->getMomentum();
256  theSpin -= theProjectileRemnant->getAngularMomentum();
257  }
258 
259  // Subtract orbital angular momentum
261  theSpin -= (thePosition-initialCenterOfMass).vector(theMomentum);
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:" << std::endl
292  << "Inside:" << std::endl;
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 << std::endl
297  << (*p)->print();
298  counter++;
299  }
300  ss <<"Outgoing:" << std::endl;
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 
309  ParticleList const &out = theStore->getOutgoingParticles();
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:" << std::endl
318  << (*i)->print() << std::endl);
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(theStore->getBook().getCurrentTime());
344  pion->boost(beta);
345 
346  nucleon->setTableMass();
347  nucleon->setMomentum(-newMomentum);
348  nucleon->adjustEnergyFromMomentum();
349  nucleon->setEmissionTime(theStore->getBook().getCurrentTime());
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:" << std::endl
382  << (*i)->print() << std::endl);
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  << std::endl);
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" << std::endl);
411  emitInsidePions();
412  }
413 
414  return true;
415  }
416 
418  ParticleList const &out = theStore->getOutgoingParticles();
419  ParticleList clusters;
420  for(ParticleIter i=out.begin(), e=out.end(); i!=e; ++i) {
421  if((*i)->isCluster()) clusters.push_back((*i));
422  }
423  if(clusters.empty()) return false;
424 
425  for(ParticleIter i=clusters.begin(), e=clusters.end(); i!=e; ++i) {
426  Cluster *cluster = dynamic_cast<Cluster*>(*i); // Can't avoid using a cast here
427 // assert(cluster);
428 #ifdef INCLXX_IN_GEANT4_MODE
429  if(!cluster)
430  continue;
431 #endif
432  cluster->deleteParticles(); // Don't need them
433  ParticleList decayProducts = ClusterDecay::decay(cluster);
434  for(ParticleIter j=decayProducts.begin(), end=decayProducts.end(); j!=end; ++j)
435  theStore->addToOutgoing(*j);
436  }
437  return true;
438  }
439 
441  // Do the phase-space decay only if Z=0 or Z=A
442  if(theA<=1 || (theZ!=0 && theA!=theZ))
443  return false;
444 
445  ParticleList decayProducts = ClusterDecay::decay(this);
446  for(ParticleIter j=decayProducts.begin(), e=decayProducts.end(); j!=e; ++j)
447  theStore->addToOutgoing(*j);
448 
449  return true;
450  }
451 
453  /* Forcing emissions of all pions in the nucleus. This probably violates
454  * energy conservation (although the computation of the recoil kinematics
455  * might sweep this under the carpet).
456  */
457  INCL_WARN("Forcing emissions of all pions in the nucleus." << std::endl);
458 
459  // Emit the pions with this kinetic energy
460  const G4double tinyPionEnergy = 0.1; // MeV
461 
462  // Push out the emitted pions
463  ParticleList const &inside = theStore->getParticles();
464  ParticleList toEject;
465  for(ParticleIter i=inside.begin(), e=inside.end(); i!=e; ++i) {
466  if((*i)->isPion()) {
467  Particle * const thePion = *i;
468  INCL_DEBUG("Forcing emission of the following particle: "
469  << thePion->print() << std::endl);
470  thePion->setEmissionTime(theStore->getBook().getCurrentTime());
471  // Correction for real masses
472  const G4double theQValueCorrection = thePion->getEmissionQValueCorrection(theA,theZ);
473  const G4double kineticEnergyOutside = thePion->getKineticEnergy() - thePion->getPotentialEnergy() + theQValueCorrection;
474  thePion->setTableMass();
475  if(kineticEnergyOutside > 0.0)
476  thePion->setEnergy(thePion->getMass()+kineticEnergyOutside);
477  else
478  thePion->setEnergy(thePion->getMass()+tinyPionEnergy);
479  thePion->adjustMomentumFromEnergy();
480  thePion->setPotentialEnergy(0.);
481  theZ -= thePion->getZ();
482  toEject.push_back(thePion);
483  }
484  }
485  for(ParticleIter i=toEject.begin(), e=toEject.end(); i!=e; ++i) {
486  theStore->particleHasBeenEjected(*i);
487  theStore->addToOutgoing(*i);
488  }
489  }
490 
492 
493  Book const &theBook = theStore->getBook();
494  const G4int nEventCollisions = theBook.getAcceptedCollisions();
495  const G4int nEventDecays = theBook.getAcceptedDecays();
496  const G4int nEventClusters = theBook.getEmittedClusters();
497  if(nEventCollisions==0 && nEventDecays==0 && nEventClusters==0)
498  return true;
499 
500  return false;
501 
502  }
503 
504  void Nucleus::computeOneNucleonRecoilKinematics() {
505  // We should be here only if the nucleus contains only one nucleon
506 // assert(theStore->getParticles().size()==1);
507 
508  INCL_ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << std::endl);
509 
510  // No excitation energy!
511  theExcitationEnergy = 0.0;
512 
513  // Move the nucleon to the outgoing list
514  Particle *remN = theStore->getParticles().front();
515  theA -= remN->getA();
516  theZ -= remN->getZ();
517  theStore->particleHasBeenEjected(remN);
518  theStore->addToOutgoing(remN);
519  remN->setEmissionTime(theStore->getBook().getCurrentTime());
520 
521  // Treat the special case of a remaining delta
522  if(remN->isDelta()) {
523  IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
524  FinalState *fs = decay->getFinalState();
525  // Eject the pion
526  ParticleList const &created = fs->getCreatedParticles();
527  for(ParticleIter j=created.begin(), e=created.end(); j!=e; ++j)
528  theStore->addToOutgoing(*j);
529  delete fs;
530  delete decay;
531  }
532 
533  // Do different things depending on how many outgoing particles we have
534  ParticleList const &outgoing = theStore->getOutgoingParticles();
535  if(outgoing.size() == 2) {
536 
537  INCL_DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << std::endl);
538 
539  // Can apply exact 2-body kinematics here. Keep the CM emission angle of
540  // the first particle.
541  Particle *p1 = outgoing.front(), *p2 = outgoing.back();
542  const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
543  // Boost to the initial CM
544  p1->boost(aBoostVector);
545  const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
546  const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
547  const G4double scale = pcm/(p1->getMomentum().mag());
548  // Reset the momenta
549  p1->setMomentum(p1->getMomentum()*scale);
550  p2->setMomentum(-p1->getMomentum());
551  p1->adjustEnergyFromMomentum();
552  p2->adjustEnergyFromMomentum();
553  // Unboost
554  p1->boost(-aBoostVector);
555  p2->boost(-aBoostVector);
556 
557  } else {
558 
559  INCL_DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << std::endl);
560 
561  const G4int maxIterations=8;
562  G4double totalEnergy, energyScale;
563  G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
564  ThreeVector totalMomentum, deltaP;
565  std::vector<ThreeVector> minMomenta; // use it to store the particle momenta that minimize the merit function
566 
567  // Reserve the vector size
568  minMomenta.reserve(outgoing.size());
569 
570  // Compute the initial total momentum
571  totalMomentum.setX(0.0);
572  totalMomentum.setY(0.0);
573  totalMomentum.setZ(0.0);
574  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
575  totalMomentum += (*i)->getMomentum();
576 
577  // Compute the initial total energy
578  totalEnergy = 0.0;
579  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
580  totalEnergy += (*i)->getEnergy();
581 
582  // Iterative algorithm starts here:
583  for(G4int iterations=0; iterations < maxIterations; ++iterations) {
584 
585  // Save the old merit-function values
586  oldOldOldVal = oldOldVal;
587  oldOldVal = oldVal;
588  oldVal = val;
589 
590  if(iterations%2 == 0) {
591  INCL_DEBUG("Momentum step" << std::endl);
592  // Momentum step: modify all the particle momenta
593  deltaP = incomingMomentum - totalMomentum;
594  G4double pOldTot = 0.0;
595  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
596  pOldTot += (*i)->getMomentum().mag();
597  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
598  const ThreeVector mom = (*i)->getMomentum();
599  (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
600  (*i)->adjustEnergyFromMomentum();
601  }
602  } else {
603  INCL_DEBUG("Energy step" << std::endl);
604  // Energy step: modify all the particle momenta
605  energyScale = initialEnergy/totalEnergy;
606  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
607  const ThreeVector mom = (*i)->getMomentum();
608  G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
609  if(pScale>0) {
610  (*i)->setEnergy((*i)->getEnergy()*energyScale);
611  (*i)->adjustMomentumFromEnergy();
612  }
613  }
614  }
615 
616  // Compute the current total momentum and energy
617  totalMomentum.setX(0.0);
618  totalMomentum.setY(0.0);
619  totalMomentum.setZ(0.0);
620  totalEnergy = 0.0;
621  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i) {
622  totalMomentum += (*i)->getMomentum();
623  totalEnergy += (*i)->getEnergy();
624  }
625 
626  // Merit factor
627  val = std::pow(totalEnergy - initialEnergy,2) +
628  0.25*(totalMomentum - incomingMomentum).mag2();
629  INCL_DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << std::endl);
630 
631  // Store the minimum
632  if(val < oldVal) {
633  INCL_DEBUG("New minimum found, storing the particle momenta" << std::endl);
634  minMomenta.clear();
635  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i)
636  minMomenta.push_back((*i)->getMomentum());
637  }
638 
639  // Stop the algorithm if the search diverges
640  if(val > oldOldVal && oldVal > oldOldOldVal) {
641  INCL_DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << std::endl);
642  break;
643  }
644  }
645 
646  // We should have made at least one successful iteration here
647 // assert(minMomenta.size()==outgoing.size());
648 
649  // Apply the optimal momenta
650  INCL_DEBUG("Applying the solution" << std::endl);
651  std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
652  for(ParticleIter i=outgoing.begin(), e=outgoing.end(); i!=e; ++i, ++v) {
653  (*i)->setMomentum(*v);
654  (*i)->adjustEnergyFromMomentum();
655  INCL_DATABLOCK((*i)->print());
656  }
657 
658  }
659 
660  }
661 
663  eventInfo->nParticles = 0;
664  G4bool isNucleonAbsorption = false;
665 
666  G4bool isPionAbsorption = false;
667  // It is possible to have pion absorption event only if the
668  // projectile is pion.
669  if(eventInfo->projectileType == PiPlus ||
670  eventInfo->projectileType == PiMinus ||
671  eventInfo->projectileType == PiZero) {
672  isPionAbsorption = true;
673  }
674 
675  // Forced CN
676  eventInfo->forcedCompoundNucleus = tryCN;
677 
678  // Outgoing particles
679  ParticleList const &outgoingParticles = getStore()->getOutgoingParticles();
680 
681  // Check if we have a nucleon absorption event: nucleon projectile
682  // and no ejected particles.
683  if(outgoingParticles.size() == 0 &&
684  (eventInfo->projectileType == Proton ||
685  eventInfo->projectileType == Neutron)) {
686  isNucleonAbsorption = true;
687  }
688 
689  // Reset the remnant counter
690  eventInfo->nRemnants = 0;
691  eventInfo->history.clear();
692 
693  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
694  // We have a pion absorption event only if the projectile is
695  // pion and there are no ejected pions.
696  if(isPionAbsorption) {
697  if((*i)->isPion()) {
698  isPionAbsorption = false;
699  }
700  }
701 
702  eventInfo->A[eventInfo->nParticles] = (*i)->getA();
703  eventInfo->Z[eventInfo->nParticles] = (*i)->getZ();
704  eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
705  eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
706  ThreeVector mom = (*i)->getMomentum();
707  eventInfo->px[eventInfo->nParticles] = mom.getX();
708  eventInfo->py[eventInfo->nParticles] = mom.getY();
709  eventInfo->pz[eventInfo->nParticles] = mom.getZ();
710  eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
711  eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
712  eventInfo->origin[eventInfo->nParticles] = -1;
713  eventInfo->history.push_back("");
714  eventInfo->nParticles++;
715  }
716  eventInfo->nucleonAbsorption = isNucleonAbsorption;
717  eventInfo->pionAbsorption = isPionAbsorption;
718  eventInfo->nCascadeParticles = eventInfo->nParticles;
719 
720  // Projectile-like remnant characteristics
721  if(theProjectileRemnant && theProjectileRemnant->getA()>0) {
722  eventInfo->ARem[eventInfo->nRemnants] = theProjectileRemnant->getA();
723  eventInfo->ZRem[eventInfo->nRemnants] = theProjectileRemnant->getZ();
724  G4double eStar = theProjectileRemnant->getExcitationEnergy();
725  if(std::abs(eStar)<1E-10)
726  eStar = 0.0; // blame rounding and set the excitation energy to zero
727  eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
728  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
729  INCL_WARN("Negative excitation energy in projectile-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << std::endl);
730  }
731  const ThreeVector &spin = theProjectileRemnant->getSpin();
732  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
733  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
734  } else { // odd-A nucleus
735  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
736  }
737  eventInfo->EKinRem[eventInfo->nRemnants] = theProjectileRemnant->getKineticEnergy();
738  const ThreeVector &mom = theProjectileRemnant->getMomentum();
739  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
740  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
741  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
742  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
743  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
744  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
745  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
746  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
747  eventInfo->nRemnants++;
748  }
749 
750  // Target-like remnant characteristics
751  if(hasRemnant()) {
752  eventInfo->ARem[eventInfo->nRemnants] = getA();
753  eventInfo->ZRem[eventInfo->nRemnants] = getZ();
754  eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
755  if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
756  INCL_WARN("Negative excitation energy in target-like remnant! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << std::endl);
757  }
758  const ThreeVector &spin = getSpin();
759  if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
760  eventInfo->JRem[eventInfo->nRemnants] = (G4int) (spin.mag()/PhysicalConstants::hc + 0.5);
761  } else { // odd-A nucleus
762  eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (spin.mag()/PhysicalConstants::hc)) + 0.5;
763  }
764  eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
765  const ThreeVector &mom = getMomentum();
766  eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
767  eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
768  eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
769  eventInfo->jxRem[eventInfo->nRemnants] = spin.getX() / PhysicalConstants::hc;
770  eventInfo->jyRem[eventInfo->nRemnants] = spin.getY() / PhysicalConstants::hc;
771  eventInfo->jzRem[eventInfo->nRemnants] = spin.getZ() / PhysicalConstants::hc;
772  eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
773  eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
774  eventInfo->nRemnants++;
775  }
776 
777  // Global counters, flags, etc.
778  Book const &theBook = theStore->getBook();
779  eventInfo->nCollisions = theBook.getAcceptedCollisions();
780  eventInfo->nBlockedCollisions = theBook.getBlockedCollisions();
781  eventInfo->nDecays = theBook.getAcceptedDecays();
782  eventInfo->nBlockedDecays = theBook.getBlockedDecays();
783  eventInfo->firstCollisionTime = theBook.getFirstCollisionTime();
784  eventInfo->firstCollisionXSec = theBook.getFirstCollisionXSec();
788  eventInfo->nReflectionAvatars = theBook.getAvatars(SurfaceAvatarType);
789  eventInfo->nCollisionAvatars = theBook.getAvatars(CollisionAvatarType);
790  eventInfo->nDecayAvatars = theBook.getAvatars(DecayAvatarType);
792  }
793 
794  Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
795  ConservationBalance theBalance;
796  // Initialise balance variables with the incoming values
797  theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
798  theBalance.A = theEventInfo.Ap + theEventInfo.At;
799 
800  theBalance.energy = getInitialEnergy();
801  theBalance.momentum = getIncomingMomentum();
802 
803  // Process outgoing particles
804  ParticleList const &outgoingParticles = theStore->getOutgoingParticles();
805  for(ParticleIter i=outgoingParticles.begin(), e=outgoingParticles.end(); i!=e; ++i ) {
806  theBalance.Z -= (*i)->getZ();
807  theBalance.A -= (*i)->getA();
808  // For outgoing clusters, the total energy automatically includes the
809  // excitation energy
810  theBalance.energy -= (*i)->getEnergy(); // Note that outgoing particles should have the real mass
811  theBalance.momentum -= (*i)->getMomentum();
812  }
813 
814  // Projectile-like remnant contribution, if present
815  if(theProjectileRemnant && theProjectileRemnant->getA()>0) {
816  theBalance.Z -= theProjectileRemnant->getZ();
817  theBalance.A -= theProjectileRemnant->getA();
818  theBalance.energy -= ParticleTable::getTableMass(theProjectileRemnant->getA(),theProjectileRemnant->getZ()) +
819  theProjectileRemnant->getExcitationEnergy();
820  theBalance.energy -= theProjectileRemnant->getKineticEnergy();
821  theBalance.momentum -= theProjectileRemnant->getMomentum();
822  }
823 
824  // Target-like remnant contribution, if present
825  if(hasRemnant()) {
826  theBalance.Z -= getZ();
827  theBalance.A -= getA();
828  theBalance.energy -= ParticleTable::getTableMass(getA(),getZ()) +
830  if(afterRecoil)
831  theBalance.energy -= getKineticEnergy();
832  theBalance.momentum -= getMomentum();
833  }
834 
835  return theBalance;
836  }
837 
839  setEnergy(initialEnergy);
840  setMomentum(incomingMomentum);
841  setSpin(incomingAngularMomentum);
844  }
845 
846  void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
847  // Deal with the projectile remnant
848  const G4int prA = theProjectileRemnant->getA();
849  if(prA>=1) {
850  // Set the mass
851  const G4double aMass = theProjectileRemnant->getInvariantMass();
852  theProjectileRemnant->setMass(aMass);
853 
854  // Compute the excitation energy from the invariant mass
855  const G4double anExcitationEnergy = aMass
856  - ParticleTable::getTableMass(prA, theProjectileRemnant->getZ());
857 
858  // Set the excitation energy
859  theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
860 
861  // No spin!
862  theProjectileRemnant->setSpin(ThreeVector());
863 
864  // Set the emission time
865  theProjectileRemnant->setEmissionTime(anEmissionTime);
866  }
867  }
868 
869 }
870 
871 #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()
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:91
FinalStateValidity getValidity() const
virtual G4double getTableMass() const
Get the real particle mass.
G4double getFirstCollisionSpectatorPosition() const
Definition: G4INCLBook.hh:88
Simple class implementing De Jong&#39;s spin model for nucleus-nucleus collisions.
G4double getMass() const
Get the cached particle mass.
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:231
Float_t firstCollisionTime
Time of the first collision [fm/c].
std::string print()
G4int getEmittedClusters() const
Definition: G4INCLBook.hh:105
G4bool pion(G4int ityp)
G4int getAcceptedCollisions() const
Definition: G4INCLBook.hh:99
Float_t emissionTime[maxSizeParticles]
Emission time [fm/c].
const ThreeVector & getIncomingMomentum() const
Get the incoming momentum vector.
virtual ~Nucleus()
ParticleList const & getParticles() const
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
ifstream in
Definition: comparison.C:7
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:57
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:100
Particle * getBlockedDelta()
G4double mag2() const
void setProtonSeparationEnergy(const G4double s)
Setter for protonSeparationEnergy.
G4double getFirstCollisionTime() const
Definition: G4INCLBook.hh:82
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:103
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.
UnorderedVector< Particle * > ParticleList
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:201
Short_t Z[maxSizeParticles]
Particle charge number.
ThreeVector theSpin
void setDensity(NuclearDensity const *const d)
Setter for theDensity.
Book & getBook()
Definition: G4INCLStore.hh:237
Short_t nParticles
Number of particles in the final state.
G4bool nucleon(G4int ityp)
G4int getBlockedDecays() const
Definition: G4INCLBook.hh:102
void particleHasBeenDestroyed(Particle *const)
Definition: G4INCLStore.cc:240
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 getCurrentTime()
Definition: G4INCLBook.hh:97
G4double phi() const
void particleHasBeenEjected(Particle *const)
Definition: G4INCLStore.cc:231
Float_t pzRem[maxSizeRemnants]
Remnant momentum, z component [MeV/c].
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.
Double_t scale
Definition: plot.C:11
G4int getAcceptedDecays() const
Definition: G4INCLBook.hh:101
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:94
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:85
Int_t nRemnants
Number of remnants.
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)
Definition: G4INCLStore.cc:144
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?
G4INCL::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
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:168
G4int getEnergyViolationInteraction() const
Definition: G4INCLBook.hh:106
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.
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.