| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 H
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, G
71 AnnihilationType AType) //D
72 : Cluster(charge,mass,strangess,true),
73 theInitialZ(charge), theInitialA(mass), t
74 theNpInitial(0), theNnInitial(0),
75 theNpionplusInitial(0), theNpionminusInit
76 theNkaonplusInitial(0), theNkaonminusInit
77 theNantiprotonInitial(0),
78 initialInternalEnergy(0.),
79 incomingAngularMomentum(0.,0.,0.), incomi
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
97 // potential. This is convenient for som
98 potentialType = IsospinPotential;
99 pionPotential = true;
100 }
101
102 thePotential = NuclearPotential::createPot
103
104 ParticleTable::setProtonSeparationEnergy(t
105 ParticleTable::setNeutronSeparationEnergy(
106
107 if (theAType==PType) theDensity = NuclearD
108 else if (theAType==NType) theDensity = Nuc
109 else
110 theDensity = NuclearDensityFactory::create
111
112 theParticleSampler->setPotential(thePotent
113 theParticleSampler->setDensity(theDensity)
114
115 if(theUniverseRadius<0)
116 theUniverseRadius = theDensity->getMaxim
117 theStore = new Store(conf);
118 }
119
120 Nucleus::~Nucleus() {
121 delete theStore;
122 deleteProjectileRemnant();
123 /* We don't delete the potential and the d
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
139 delete theProjectileRemnant;
140 theProjectileRemnant = NULL;
141 Cluster::initializeParticles();
142
143 for(ParticleIter i=particles.begin(), e=pa
144 updatePotentialEnergy(*i);
145 }
146 theStore->add(particles);
147 particles.clear();
148 initialInternalEnergy = computeTotalEnergy
149 initialCenterOfMass = thePosition;
150 }
151
152 void Nucleus::applyFinalState(FinalState *fi
153 if(!finalstate) // do nothing if no final
154 return;
155
156 G4double totalEnergy = 0.0;
157
158 FinalStateValidity const validity = finals
159 if(validity == ValidFS) {
160
161 ParticleList const &created = finalstate
162 for(ParticleIter iter=created.begin(), e
163 theStore->add((*iter));
164 if(!(*iter)->isOutOfWell()) {
165 totalEnergy += (*iter)->getEnergy()
166 }
167 }
168
169 ParticleList const &deleted = finalstate
170 for(ParticleIter iter=deleted.begin(), e
171 theStore->particleHasBeenDestroyed(*it
172 }
173
174 ParticleList const &modified = finalstat
175 for(ParticleIter iter=modified.begin(),
176 theStore->particleHasBeenUpdated(*iter
177 totalEnergy += (*iter)->getEnergy() -
178 }
179
180 ParticleList const &out = finalstate->ge
181 for(ParticleIter iter=out.begin(), e=out
182 if((*iter)->isCluster()) {
183 Cluster *clusterOut = dynamic_cast<Clust
184 // assert(clusterOut);
185 #ifdef INCLXX_IN_GEANT4_MODE
186 if(!clusterOut)
187 continue;
188 #endif
189 ParticleList const &components = clu
190 for(ParticleIter in=components.begin
191 theStore->particleHasBeenEjected(*
192 } else {
193 theStore->particleHasBeenEjected(*it
194 }
195 totalEnergy += (*iter)->getEnergy();
196 theA -= (*iter)->getA();
197 theZ -= (*iter)->getZ();
198 theS -= (*iter)->getS();
199 theStore->addToOutgoing(*iter);
200 (*iter)->setEmissionTime(theStore->get
201 }
202
203 ParticleList const &entering = finalstat
204 for(ParticleIter iter=entering.begin(),
205 insertParticle(*iter);
206 totalEnergy += (*iter)->getEnergy() -
207 }
208
209 // actually perform the removal of the s
210 theStore->removeScheduledAvatars();
211 } else if(validity == ParticleBelowFermiFS
212 INCL_DEBUG("A Particle is entering below
213 tryCN = true;
214 ParticleList const &entering = finalstat
215 for(ParticleIter iter=entering.begin(),
216 insertParticle(*iter);
217 }
218 }
219
220 if(validity==ValidFS &&
221 std::abs(totalEnergy - finalstate->get
222 INCL_ERROR("Energy nonconservation! Ener
223 << finalstate->getTotalEnergyBeforeI
224 <<" and after interaction = "
225 << totalEnergy << '\n'
226 << finalstate->print());
227 }
228 }
229
230 void Nucleus::propagateParticles(G4double /*
231 INCL_WARN("Useless Nucleus::propagateParti
232 }
233
234 G4double Nucleus::computeTotalEnergy() const
235 G4double totalEnergy = 0.0;
236 ParticleList const &inside = theStore->get
237 for(ParticleIter p=inside.begin(), e=insid
238 if((*p)->isNucleon()) // Ugly: we should
239 totalEnergy += (*p)->getKineticEnergy(
240 else if((*p)->isResonance())
241 totalEnergy += (*p)->getEnergy() - (*p
242 else if((*p)->isHyperon())
243 totalEnergy += (*p)->getEnergy() - (*p
244 else if((*p)->isAntiNucleon())
245 totalEnergy += (*p)->getEnergy() - (*p
246 else if((*p)->isAntiLambda())
247 totalEnergy += (*p)->getEnergy() - (*p
248 //std::cout << ParticleTable::getRealM
249 else
250 totalEnergy += (*p)->getEnergy() - (*p
251 }
252
253 return totalEnergy;
254 }
255
256 void Nucleus::computeRecoilKinematics() {
257 // If the remnant consists of only one nuc
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
267 theMomentum = incomingMomentum;
268 theSpin = incomingAngularMomentum;
269
270 ParticleList const &outgoing = theStore->g
271 for(ParticleIter p=outgoing.begin(), e=out
272 theMomentum -= (*p)->getMomentum();
273 theSpin -= (*p)->getAngularMomentum();
274 }
275 if(theProjectileRemnant) {
276 theMomentum -= theProjectileRemnant->get
277 theSpin -= theProjectileRemnant->getAngu
278 }
279
280 // Subtract orbital angular momentum
281 thePosition = computeCenterOfMass();
282 theSpin -= (thePosition-initialCenterOfMas
283
284 setMass(ParticleTable::getTableMass(theA,t
285 adjustEnergyFromMomentum();
286 remnant=true;
287 }
288
289 ThreeVector Nucleus::computeCenterOfMass() c
290 ThreeVector cm(0.,0.,0.);
291 G4double totalMass = 0.0;
292 ParticleList const &inside = theStore->get
293 for(ParticleIter p=inside.begin(), e=insid
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()
303 const G4double totalEnergy = computeTotalE
304 const G4double separationEnergies = comput
305
306 G4double eSep = 0;
307 if (getAType() == AnnihilationType::Def) {
308 } else if (getAType() == AnnihilationType::PT
309 } else if (getAType() == AnnihilationType::NT
310 } else if (getAType() == AnnihilationType::PT
311 eSep = ParticleTable::getProtonSeparationEne
312 } else if (getAType() == AnnihilationType::NT
313 eSep = ParticleTable::getNeutronSeparationEn
314 } else if (getAType() == AnnihilationType::Nb
315 eSep = ParticleTable::getProtonSeparationEne
316 } else if (getAType() == AnnihilationType::Nb
317 eSep = ParticleTable::getNeutronSeparationEn
318 }
319
320 if (eSep > 0. && (totalEnergy - initialInter
321 INCL_DEBUG("Negative Excitation Energy du
322 }
323
324 return totalEnergy - initialInternalEnergy -
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->get
337 for(ParticleIter p=inside.begin(), e=insid
338 ss << "index = " << counter << '\n'
339 << (*p)->print();
340 counter++;
341 }
342 ss <<"Outgoing:" << '\n';
343 ParticleList const &outgoing = theStore->g
344 for(ParticleIter p=outgoing.begin(), e=out
345 ss << (*p)->print();
346
347 return ss.str();
348 }
349
350 G4bool Nucleus::decayOutgoingDeltas() {
351 ParticleList const &out = theStore->getOut
352 ParticleList deltas;
353 for(ParticleIter i=out.begin(), e=out.end(
354 if((*i)->isDelta()) deltas.push_back((*i
355 }
356 if(deltas.empty()) return false;
357
358 for(ParticleIter i=deltas.begin(), e=delta
359 INCL_DEBUG("Decay outgoing delta particl
360 << (*i)->print() << '\n');
361 const ThreeVector beta = -(*i)->boostVec
362 const G4double deltaMass = (*i)->getMass
363
364 // Set the delta momentum to zero and sa
365 // This makes life simpler if we are usi
366 (*i)->setMomentum(ThreeVector());
367 (*i)->setEnergy((*i)->getMass());
368
369 // Use a DecayAvatar
370 IAvatar *decay = new DecayAvatar((*i), 0
371 FinalState *fs = decay->getFinalState();
372 Particle * const pion = fs->getCreatedPa
373 Particle * const nucleon = fs->getModifi
374
375 // Adjust the decay momentum if we are u
376 const G4double decayMomentum = Kinematic
377 nucleon->getTableMass(),
378 pion->getTableMass());
379 ThreeVector newMomentum = pion->getMomen
380 newMomentum *= decayMomentum / newMoment
381
382 pion->setTableMass();
383 pion->setMomentum(newMomentum);
384 pion->adjustEnergyFromMomentum();
385 pion->setEmissionTime(nucleon->getEmissi
386 pion->boost(beta);
387 pion->setBiasCollisionVector(nucleon->ge
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 nothin
405 * excitation energy).
406 * If, however, the remnant is unphysical
407 * decay and get rid of all the pions. In
408 * end up with Z<0 or Z>A if the remnant c
409 * more pi+ than neutrons, respectively.
410 */
411 const G4bool unphysicalRemnant = (theZ<0 |
412 if(thePotential->hasPionPotential() && !un
413 return false;
414
415 // Build a list of deltas (avoid modifying
416 ParticleList const &inside = theStore->get
417 ParticleList deltas;
418 for(ParticleIter i=inside.begin(), e=insid
419 if((*i)->isDelta()) deltas.push_back((*i
420
421 // Loop over the deltas, make them decay
422 for(ParticleIter i=deltas.begin(), e=delta
423 INCL_DEBUG("Decay inside delta particle:
424 << (*i)->print() << '\n');
425 // Create a forced-decay avatar. Note th
426 // also that if the remnant is unphysica
427 // up energy conservation and CDPP by pa
428 // nucleus.
429 IAvatar *decay;
430 if(unphysicalRemnant) {
431 INCL_WARN("Forcing delta decay inside
432 << ", Z=" << theZ << "). Might le
433 << '\n');
434 decay = new DecayAvatar((*i), 0.0, NUL
435 } else
436 decay = new DecayAvatar((*i), 0.0, thi
437 FinalState *fs = decay->getFinalState();
438
439 // The pion can be ejected only if we ma
440 // conservation and if pion emission doe
441 // energies.
442 if(fs->getValidity()==ValidFS) {
443 // Apply the final state to the nucleu
444 applyFinalState(fs);
445 }
446 delete fs;
447 delete decay;
448 }
449
450 // If the remnant is unphysical, emit all
451 if(unphysicalRemnant) {
452 INCL_DEBUG("Remnant is unphysical: Z=" <
453 emitInsidePions();
454 }
455
456 return true;
457 }
458
459 G4bool Nucleus::decayInsideStrangeParticles(
460
461 /* Transform each strange particles into a
462 * Every Kaon (KPlus and KZero) are emited
463 */
464 const G4bool unphysicalRemnant = (theZ<0 |
465 if(unphysicalRemnant){
466 emitInsideStrangeParticles();
467 INCL_WARN("Remnant is unphysical: Z="
468 return false;
469 }
470
471 /* Build a list of particles with a strang
472 * and two other one for proton and neutro
473 */
474 ParticleList const &inside = theStore->get
475 ParticleList stranges;
476 ParticleList protons;
477 ParticleList neutrons;
478 for(ParticleIter i=inside.begin(), e=insid
479 if((*i)->isSigma() || (*i)->isAntiKaon
480 else if((*i)->isNucleon() && (*i)->get
481 else if((*i)->isNucleon() && (*i)->get
482 }
483
484 if((stranges.size() > protons.size()) || (
485 INCL_WARN("Remnant is unphysical: Npro
486 emitInsideStrangeParticles();
487 return false;
488 }
489
490 // Loop over the strange particles, make t
491 ParticleIter protonIter = protons.begin();
492 ParticleIter neutronIter = neutrons.begin(
493 for(ParticleIter i=stranges.begin(), e=str
494 INCL_DEBUG("Absorbe inside strange parti
495 << (*i)->print() << '\n');
496 IAvatar *decay;
497 if((*i)->getType() == SigmaMinus){
498 decay = new DecayAvatar((*protonIter
499 ++protonIter;
500 }
501 else if((*i)->getType() == SigmaPlus){
502 decay = new DecayAvatar((*neutronIte
503 ++neutronIter;
504 }
505 else if(Random::shoot()*(protons.size()
506 decay = new DecayAvatar((*protonIter
507 ++protonIter;
508 }
509 else {
510 decay = new DecayAvatar((*neutronIte
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(
523 ParticleList const &out = theStore->ge
524 ParticleList pionResonances;
525 for(ParticleIter i=out.begin(), e=out.
526 // if((*i)->isEta() || (*i)->isOmeg
527 if(((*i)->isEta() && timeThreshold
528 }
529 if(pionResonances.empty()) return fals
530
531 for(ParticleIter i=pionResonances.begi
532 INCL_DEBUG("Decay outgoing pionRes
533 << (*i)->print() << '\n
534 const ThreeVector beta = -(*i)->bo
535 const G4double pionResonanceMass =
536
537 // Set the pionResonance momentum
538 // This makes life simpler if we a
539 (*i)->setMomentum(ThreeVector());
540 (*i)->setEnergy((*i)->getMass());
541
542 // Use a DecayAvatar
543 IAvatar *decay = new DecayAvatar((
544 FinalState *fs = decay->getFinalSt
545
546 Particle * const theModifiedPartic
547 ParticleList const &created = fs->
548 Particle * const theCreatedParticl
549
550 if (created.size() == 1) {
551
552 // Adjust the decay momentum i
553 const G4double decayMomentum =
554 ThreeVector newMomentum = theC
555 newMomentum *= decayMomentum /
556
557 theCreatedParticle1->setTableM
558 theCreatedParticle1->setMoment
559 theCreatedParticle1->adjustEne
560 //theCreatedParticle1->setEmis
561 theCreatedParticle1->boost(bet
562 theCreatedParticle1->setBiasCo
563
564 theModifiedParticle->setTableM
565 theModifiedParticle->setMoment
566 theModifiedParticle->adjustEne
567 theModifiedParticle->boost(bet
568
569 theStore->addToOutgoing(theCre
570 }
571 else if (created.size() == 2) {
572 Particle * const theCreatedPar
573
574 theCreatedParticle1->boost(bet
575 theCreatedParticle1->setBiasCo
576 theCreatedParticle2->boost(bet
577 theCreatedParticle2->setBiasCo
578 theModifiedParticle->boost(bet
579
580 theStore->addToOutgoing(theCre
581 theStore->addToOutgoing(theCre
582 }
583 else {
584 INCL_ERROR("Wrong number (< 2)
585 }
586 delete fs;
587 delete decay;
588 }
589
590 return true;
591 }
592
593 G4bool Nucleus::decayOutgoingSigmaZero(G4dou
594 ParticleList const &out = theStore->ge
595 ParticleList neutralsigma;
596 for(ParticleIter i=out.begin(), e=out.
597 if((*i)->getType() == SigmaZero &&
598 }
599 if(neutralsigma.empty()) return false;
600
601 for(ParticleIter i=neutralsigma.begin(
602 INCL_DEBUG("Decay outgoing neutral
603 << (*i)->print() << '\n
604 const ThreeVector beta = -(*i)->bo
605 const G4double neutralsigmaMass =
606
607 // Set the neutral sigma momentum
608 // This makes life simpler if we a
609 (*i)->setMomentum(ThreeVector());
610 (*i)->setEnergy((*i)->getMass());
611
612 // Use a DecayAvatar
613 IAvatar *decay = new DecayAvatar((
614 FinalState *fs = decay->getFinalSt
615
616 Particle * const theModifiedPartic
617 ParticleList const &created = fs->
618 Particle * const theCreatedParticl
619
620 if (created.size() == 1) {
621
622 // Adjust the decay momentum i
623 const G4double decayMomentum =
624 ThreeVector newMomentum = theC
625 newMomentum *= decayMomentum /
626
627 theCreatedParticle->setTableMa
628 theCreatedParticle->setMomentu
629 theCreatedParticle->adjustEner
630 theCreatedParticle->boost(beta
631 theCreatedParticle->setBiasCol
632
633 theModifiedParticle->setTableM
634 theModifiedParticle->setMoment
635 theModifiedParticle->adjustEne
636 theModifiedParticle->boost(bet
637
638 theStore->addToOutgoing(theCre
639 }
640 else {
641 INCL_ERROR("Wrong number (!= 1
642 }
643 delete fs;
644 delete decay;
645 }
646
647 return true;
648 }
649
650 G4bool Nucleus::decayOutgoingNeutralKaon() {
651 ParticleList const &out = theStore->ge
652 ParticleList neutralkaon;
653 for(ParticleIter i=out.begin(), e=out.
654 if((*i)->getType() == KZero || (*
655 }
656 if(neutralkaon.empty()) return false;
657
658 for(ParticleIter i=neutralkaon.begin()
659 INCL_DEBUG("Transform outgoing neu
660 << (*i)->print() << '\n
661
662 // Use a DecayAvatar
663 IAvatar *decay = new DecayAvatar((
664 FinalState *fs = decay->getFinalSt
665
666 delete fs;
667 delete decay;
668 }
669
670 return true;
671 }
672
673 G4bool Nucleus::decayOutgoingClusters() {
674 ParticleList const &out = theStore->getOut
675 ParticleList clusters;
676 for(ParticleIter i=out.begin(), e=out.end(
677 if((*i)->isCluster()) clusters.push_back
678 }
679 if(clusters.empty()) return false;
680
681 for(ParticleIter i=clusters.begin(), e=clu
682 Cluster *cluster = dynamic_cast<Cluster*
683 // assert(cluster);
684 #ifdef INCLXX_IN_GEANT4_MODE
685 if(!cluster)
686 continue;
687 #endif
688 cluster->deleteParticles(); // Don't nee
689 ParticleList decayProducts = ClusterDeca
690 for(ParticleIter j=decayProducts.begin()
691 (*j)->setBiasCollisionVector(cluster->
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
700 if(theA<=1 || (theZ!=0 && (theA+theS)!=the
701 return false;
702
703 ParticleList decayProducts = ClusterDecay:
704 for(ParticleIter j=decayProducts.begin(),
705 (*j)->setBiasCollisionVector(this->getBi
706 theStore->addToOutgoing(*j);
707 }
708
709 return true;
710 }
711
712 void Nucleus::emitInsidePions() {
713 /* Forcing emissions of all pions in the n
714 * energy conservation (although the compu
715 * might sweep this under the carpet).
716 */
717 INCL_WARN("Forcing emissions of all pions
718
719 // Emit the pions with this kinetic energy
720 const G4double tinyPionEnergy = 0.1; // Me
721
722 // Push out the emitted pions
723 ParticleList const &inside = theStore->get
724 ParticleList toEject;
725 for(ParticleIter i=inside.begin(), e=insid
726 if((*i)->isPion()) {
727 Particle * const thePion = *i;
728 INCL_DEBUG("Forcing emission of the fo
729 << thePion->print() << '\n'
730 thePion->setEmissionTime(theStore->get
731 // Correction for real masses
732 const G4double theQValueCorrection = t
733 const G4double kineticEnergyOutside =
734 thePion->setTableMass();
735 if(kineticEnergyOutside > 0.0)
736 thePion->setEnergy(thePion->getMass(
737 else
738 thePion->setEnergy(thePion->getMass(
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=toEj
746 theStore->particleHasBeenEjected(*i);
747 theStore->addToOutgoing(*i);
748 (*i)->setParticleBias(Particle::getTotal
749 }
750 }
751
752 void Nucleus::emitInsideStrangeParticles() {
753 /* Forcing emissions of Sigmas and antiKao
754 * This probably violates energy conservat
755 * (although the computation of the recoil
756 * might sweep this under the carpet).
757 */
758 INCL_DEBUG("Forcing emissions of all stran
759
760 // Emit the strange particles with this ki
761 const G4double tinyEnergy = 0.1; // MeV
762
763 // Push out the emitted strange particles
764 ParticleList const &inside = theStore->get
765 ParticleList toEject;
766 for(ParticleIter i=inside.begin(), e=insid
767 if((*i)->isSigma() || (*i)->isAntiKaon()
768 Particle * const theParticle = *i;
769 INCL_DEBUG("Forcing emission of the fo
770 << theParticle->print() <<
771 theParticle->setEmissionTime(theStore-
772 // Correction for real masses
773 const G4double theQValueCorrection = t
774 const G4double kineticEnergyOutside =
775 theParticle->setTableMass();
776 if(kineticEnergyOutside > 0.0)
777 theParticle->setEnergy(theParticle->
778 else
779 theParticle->setEnergy(theParticle->
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=toEj
789 theStore->particleHasBeenEjected(*i);
790 theStore->addToOutgoing(*i);
791 (*i)->setParticleBias(Particle::getTotal
792 }
793 }
794
795 G4int Nucleus::emitInsideLambda() {
796 /* Forcing emissions of all Lambda in the
797 * This probably violates energy conservat
798 * (although the computation of the recoil
799 * might sweep this under the carpet).
800 */
801 INCL_DEBUG("Forcing emissions of all Lambd
802
803 // Emit the Lambda with this kinetic energ
804 const G4double tinyEnergy = 0.1; // MeV
805
806 // Push out the emitted Lambda
807 ParticleList const &inside = theStore->get
808 ParticleList toEject;
809 for(ParticleIter i=inside.begin(), e=insid
810 if((*i)->isLambda()) {
811 Particle * const theLambda = *i;
812 INCL_DEBUG("Forcing emission of the fo
813 << theLambda->print() << '\
814 theLambda->setEmissionTime(theStore->g
815 // Correction for real masses
816 const G4double theQValueCorrection = t
817 const G4double kineticEnergyOutside =
818 theLambda->setTableMass();
819 if(kineticEnergyOutside > 0.0)
820 theLambda->setEnergy(theLambda->getM
821 else
822 theLambda->setEnergy(theLambda->getM
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=toEj
831 theStore->particleHasBeenEjected(*i);
832 theStore->addToOutgoing(*i);
833 (*i)->setParticleBias(Particle::getTotal
834 }
835 return (G4int)toEject.size();
836 }
837
838 G4bool Nucleus::emitInsideKaon() {
839 /* Forcing emissions of all Kaon (not anti
840 * This probably violates energy conservat
841 * (although the computation of the recoil
842 * might sweep this under the carpet).
843 */
844 INCL_DEBUG("Forcing emissions of all Kaon
845
846 // Emit the Kaon with this kinetic energy
847 const G4double tinyEnergy = 0.1; // MeV
848
849 // Push out the emitted kaon
850 ParticleList const &inside = theStore->get
851 ParticleList toEject;
852 for(ParticleIter i=inside.begin(), e=insid
853 if((*i)->isKaon()) {
854 Particle * const theKaon = *i;
855 INCL_DEBUG("Forcing emission of the fo
856 << theKaon->print() << '\n'
857 theKaon->setEmissionTime(theStore->get
858 // Correction for real masses
859 const G4double theQValueCorrection = t
860 const G4double kineticEnergyOutside =
861 theKaon->setTableMass();
862 if(kineticEnergyOutside > 0.0)
863 theKaon->setEnergy(theKaon->getMass(
864 else
865 theKaon->setEnergy(theKaon->getMass(
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=toEj
874 theStore->particleHasBeenEjected(*i);
875 theStore->addToOutgoing(*i);
876 (*i)->setParticleBias(Particle::getTotal
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.get
886 const G4int nEventDecays = theBook.getAcce
887 const G4int nEventClusters = theBook.getEm
888 if(nEventCollisions==0 && nEventDecays==0
889 return true;
890
891 return false;
892
893 }
894
895 void Nucleus::computeOneNucleonRecoilKinemat
896 // We should be here only if the nucleus c
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().
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().
910
911 // Treat the special case of a remaining d
912 if(remN->isDelta()) {
913 IAvatar *decay = new DecayAvatar(remN, 0
914 FinalState *fs = decay->getFinalState();
915 // Eject the pion
916 ParticleList const &created = fs->getCre
917 for(ParticleIter j=created.begin(), e=cr
918 theStore->addToOutgoing(*j);
919 delete fs;
920 delete decay;
921 }
922
923 // Do different things depending on how ma
924 ParticleList const &outgoing = theStore->g
925 if(outgoing.size() == 2) {
926
927 INCL_DEBUG("Two particles in the outgoin
928
929 // Can apply exact 2-body kinematics her
930 // the first particle.
931 Particle *p1 = outgoing.front(), *p2 = o
932 const ThreeVector aBoostVector = incomin
933 // Boost to the initial CM
934 p1->boost(aBoostVector);
935 const G4double sqrts = std::sqrt(initial
936 const G4double pcm = KinematicsUtils::mo
937 const G4double scale = pcm/(p1->getMomen
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
950
951 const G4int maxIterations=8;
952 G4double totalEnergy, energyScale;
953 G4double val=1.E+100, oldVal=1.E+100, ol
954 ThreeVector totalMomentum, deltaP;
955 std::vector<ThreeVector> minMomenta; //
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=o
965 totalMomentum += (*i)->getMomentum();
966
967 // Compute the initial total energy
968 totalEnergy = 0.0;
969 for(ParticleIter i=outgoing.begin(), e=o
970 totalEnergy += (*i)->getEnergy();
971
972 // Iterative algorithm starts here:
973 for(G4int iterations=0; iterations < max
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 par
983 deltaP = incomingMomentum - totalMom
984 G4double pOldTot = 0.0;
985 for(ParticleIter i=outgoing.begin(),
986 pOldTot += (*i)->getMomentum().mag
987 for(ParticleIter i=outgoing.begin(),
988 const ThreeVector mom = (*i)->getM
989 (*i)->setMomentum(mom + deltaP*mom
990 (*i)->adjustEnergyFromMomentum();
991 }
992 } else {
993 INCL_DEBUG("Energy step" << '\n');
994 // Energy step: modify all the parti
995 energyScale = initialEnergy/totalEne
996 for(ParticleIter i=outgoing.begin(),
997 const ThreeVector mom = (*i)->getM
998 G4double pScale = ((*i)->getEnergy
999 if(pScale>0) {
1000 (*i)->setEnergy((*i)->getEnergy
1001 (*i)->adjustMomentumFromEnergy(
1002 }
1003 }
1004 }
1005
1006 // Compute the current total momentum
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(),
1012 totalMomentum += (*i)->getMomentum(
1013 totalEnergy += (*i)->getEnergy();
1014 }
1015
1016 // Merit factor
1017 val = std::pow(totalEnergy - initialE
1018 0.25*(totalMomentum - incomingMomen
1019 INCL_DEBUG("Merit function: val=" <<
1020
1021 // Store the minimum
1022 if(val < oldVal) {
1023 INCL_DEBUG("New minimum found, stor
1024 minMomenta.clear();
1025 for(ParticleIter i=outgoing.begin()
1026 minMomenta.push_back((*i)->getMom
1027 }
1028
1029 // Stop the algorithm if the search d
1030 if(val > oldOldVal && oldVal > oldOld
1031 INCL_DEBUG("Search is diverging, br
1032 break;
1033 }
1034 }
1035
1036 // We should have made at least one suc
1037 // assert(minMomenta.size()==outgoing.size())
1038
1039 // Apply the optimal momenta
1040 INCL_DEBUG("Applying the solution" << '
1041 std::vector<ThreeVector>::const_iterato
1042 for(ParticleIter i=outgoing.begin(), e=
1043 (*i)->setMomentum(*v);
1044 (*i)->adjustEnergyFromMomentum();
1045 INCL_DATABLOCK((*i)->print());
1046 }
1047
1048 }
1049
1050 }
1051
1052 void Nucleus::fillEventInfo(EventInfo *even
1053 eventInfo->nParticles = 0;
1054 G4bool isNucleonAbsorption = false;
1055 G4bool isPionAbsorption = false;
1056 // It is possible to have pion absorption
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 = g
1069
1070 // Check if we have a nucleon absorption
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.begi
1083 // We have a pion absorption event only
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->nPar
1092
1093 #ifdef INCLXX_IN_GEANT4_MODE
1094 eventInfo->A[eventInfo->nParticles] = (
1095 eventInfo->Z[eventInfo->nParticles] = (
1096 eventInfo->S[eventInfo->nParticles] = (
1097 #else
1098 eventInfo->A[eventInfo->nParticles] = (
1099 eventInfo->Z[eventInfo->nParticles] = (
1100 eventInfo->S[eventInfo->nParticles] = (
1101 #endif
1102 eventInfo->emissionTime[eventInfo->nPar
1103 eventInfo->EKin[eventInfo->nParticles]
1104 ThreeVector mom = (*i)->getMomentum();
1105 eventInfo->px[eventInfo->nParticles] =
1106 eventInfo->py[eventInfo->nParticles] =
1107 eventInfo->pz[eventInfo->nParticles] =
1108 eventInfo->theta[eventInfo->nParticles]
1109 eventInfo->phi[eventInfo->nParticles] =
1110 eventInfo->origin[eventInfo->nParticles
1111 #ifdef INCLXX_IN_GEANT4_MODE
1112 eventInfo->parentResonancePDGCode[event
1113 eventInfo->parentResonanceID[eventInfo-
1114 #endif
1115 eventInfo->history.push_back("");
1116 if ((*i)->getType() != Composite) {
1117 ParticleSpecies pt((*i)->getType());
1118 eventInfo->PDGCode[eventInfo->nPartic
1119 }
1120 else {
1121 ParticleSpecies pt((*i)->getA(), (*i)
1122 eventInfo->PDGCode[eventInfo->nPartic
1123 }
1124 eventInfo->nParticles++;
1125 }
1126 eventInfo->nucleonAbsorption = isNucleonA
1127 eventInfo->pionAbsorption = isPionAbsorpt
1128 eventInfo->nCascadeParticles = eventInfo-
1129
1130 // Projectile-like remnant characteristic
1131 if(theProjectileRemnant && theProjectileR
1132 #ifdef INCLXX_IN_GEANT4_MODE
1133 eventInfo->ARem[eventInfo->nRemnants] =
1134 eventInfo->ZRem[eventInfo->nRemnants] =
1135 eventInfo->SRem[eventInfo->nRemnants] =
1136 #else
1137 eventInfo->ARem[eventInfo->nRemnants] =
1138 eventInfo->ZRem[eventInfo->nRemnants] =
1139 eventInfo->SRem[eventInfo->nRemnants] =
1140 #endif
1141 G4double eStar = theProjectileRemnant->
1142 if(std::abs(eStar)<1E-10)
1143 eStar = 0.0; // blame rounding and se
1144 eventInfo->EStarRem[eventInfo->nRemnant
1145 if(eventInfo->EStarRem[eventInfo->nRemn
1146 INCL_WARN("Negative excitation energy in
1147 }
1148 const ThreeVector &spin = theProjectile
1149 if(eventInfo->ARem[eventInfo->nRemnants
1150 eventInfo->JRem[eventInfo->nRemnants] = (
1151 } else { // odd-A nucleus
1152 eventInfo->JRem[eventInfo->nRemnants] = (
1153 }
1154 eventInfo->EKinRem[eventInfo->nRemnants
1155 const ThreeVector &mom = theProjectileR
1156 eventInfo->pxRem[eventInfo->nRemnants]
1157 eventInfo->pyRem[eventInfo->nRemnants]
1158 eventInfo->pzRem[eventInfo->nRemnants]
1159 eventInfo->jxRem[eventInfo->nRemnants]
1160 eventInfo->jyRem[eventInfo->nRemnants]
1161 eventInfo->jzRem[eventInfo->nRemnants]
1162 eventInfo->thetaRem[eventInfo->nRemnant
1163 eventInfo->phiRem[eventInfo->nRemnants]
1164 eventInfo->nRemnants++;
1165 }
1166
1167 // Target-like remnant characteristics
1168 if(hasRemnant()) {
1169 #ifdef INCLXX_IN_GEANT4_MODE
1170 eventInfo->ARem[eventInfo->nRemnants] =
1171 eventInfo->ZRem[eventInfo->nRemnants] =
1172 eventInfo->SRem[eventInfo->nRemnants] =
1173 #else
1174 eventInfo->ARem[eventInfo->nRemnants] =
1175 eventInfo->ZRem[eventInfo->nRemnants] =
1176 eventInfo->SRem[eventInfo->nRemnants] =
1177 #endif
1178 eventInfo->EStarRem[eventInfo->nRemnant
1179 if(eventInfo->EStarRem[eventInfo->nRemn
1180 INCL_WARN("Negative excitation energy in
1181 }
1182 const ThreeVector &spin = getSpin();
1183 if(eventInfo->ARem[eventInfo->nRemnants
1184 eventInfo->JRem[eventInfo->nRemnants] = (
1185 } else { // odd-A nucleus
1186 eventInfo->JRem[eventInfo->nRemnants] = (
1187 }
1188 eventInfo->EKinRem[eventInfo->nRemnants
1189 const ThreeVector &mom = getMomentum();
1190 eventInfo->pxRem[eventInfo->nRemnants]
1191 eventInfo->pyRem[eventInfo->nRemnants]
1192 eventInfo->pzRem[eventInfo->nRemnants]
1193 eventInfo->jxRem[eventInfo->nRemnants]
1194 eventInfo->jyRem[eventInfo->nRemnants]
1195 eventInfo->jzRem[eventInfo->nRemnants]
1196 eventInfo->thetaRem[eventInfo->nRemnant
1197 eventInfo->phiRem[eventInfo->nRemnants]
1198 eventInfo->nRemnants++;
1199 }
1200
1201 // Global counters, flags, etc.
1202 Book const &theBook = theStore->getBook()
1203 eventInfo->nCollisions = theBook.getAccep
1204 eventInfo->nBlockedCollisions = theBook.g
1205 eventInfo->nDecays = theBook.getAcceptedD
1206 eventInfo->nBlockedDecays = theBook.getBl
1207 eventInfo->firstCollisionTime = theBook.g
1208 eventInfo->firstCollisionXSec = theBook.g
1209 eventInfo->firstCollisionSpectatorPositio
1210 eventInfo->firstCollisionSpectatorMomentu
1211 eventInfo->firstCollisionIsElastic = theB
1212 eventInfo->nReflectionAvatars = theBook.g
1213 eventInfo->nCollisionAvatars = theBook.ge
1214 eventInfo->nDecayAvatars = theBook.getAva
1215 eventInfo->nEnergyViolationInteraction =
1216 }
1217
1218
1219
1220 Nucleus::ConservationBalance Nucleus::getCo
1221 ConservationBalance theBalance;
1222 // Initialise balance variables with the
1223 INCL_DEBUG("theEventInfo " << theEventInf
1224 theBalance.Z = theEventInfo.Zp + theEvent
1225 theBalance.A = theEventInfo.Ap + theEvent
1226 theBalance.S = theEventInfo.Sp + theEvent
1227 INCL_DEBUG("theBalance Z and A " << theBa
1228 theBalance.energy = getInitialEnergy();
1229 theBalance.momentum = getIncomingMomentum
1230
1231 // Process outgoing particles
1232 ParticleList const &outgoingParticles = t
1233 for(ParticleIter i=outgoingParticles.begi
1234 theBalance.Z -= (*i)->getZ();
1235 theBalance.A -= (*i)->getA();
1236 theBalance.S -= (*i)->getS();
1237 // For outgoing clusters, the total ene
1238 // excitation energy
1239 theBalance.energy -= (*i)->getEnergy();
1240 theBalance.momentum -= (*i)->getMomentu
1241 }
1242
1243 // Projectile-like remnant contribution,
1244 if(theProjectileRemnant && theProjectileR
1245 theBalance.Z -= theProjectileRemnant->g
1246 theBalance.A -= theProjectileRemnant->g
1247 theBalance.S -= theProjectileRemnant->g
1248 theBalance.energy -= ParticleTable::get
1249 theProjectileRemnant->getExcitationEn
1250 theBalance.energy -= theProjectileRemna
1251 theBalance.momentum -= theProjectileRem
1252 }
1253
1254 // Target-like remnant contribution, if p
1255 if(hasRemnant()) {
1256 theBalance.Z -= getZ();
1257 theBalance.A -= getA();
1258 theBalance.S -= getS();
1259 theBalance.energy -= ParticleTable::get
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
1274 setMass(getTableMass() + theExcitationEne
1275 }
1276
1277 void Nucleus::finalizeProjectileRemnant(con
1278 // Deal with the projectile remnant
1279 const G4int prA = theProjectileRemnant->g
1280 if(prA>=1) {
1281 // Set the mass
1282 const G4double aMass = theProjectileRem
1283 theProjectileRemnant->setMass(aMass);
1284
1285 // Compute the excitation energy from t
1286 const G4double anExcitationEnergy = aMa
1287 - ParticleTable::getTableMass(prA, th
1288
1289 // Set the excitation energy
1290 theProjectileRemnant->setExcitationEner
1291
1292 // No spin!
1293 theProjectileRemnant->setSpin(ThreeVect
1294
1295 // Set the emission time
1296 theProjectileRemnant->setEmissionTime(a
1297 }
1298 }
1299
1300 }
1301
1302 #endif
1303