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