| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA
29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA
30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA
31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege
32 // Davide Mancusi, CEA-Saclay, France <<
33 // 32 //
34 #define INCLXX_IN_GEANT4_MODE 1 33 #define INCLXX_IN_GEANT4_MODE 1
35 34
36 #include "globals.hh" 35 #include "globals.hh"
37 36
38 /* 37 /*
39 * G4INCLNucleus.cc 38 * G4INCLNucleus.cc
40 * 39 *
41 * \date Jun 5, 2009 40 * \date Jun 5, 2009
42 * \author Pekka Kaitaniemi 41 * \author Pekka Kaitaniemi
43 */ 42 */
44 43
45 #ifndef G4INCLNucleus_hh 44 #ifndef G4INCLNucleus_hh
46 #define G4INCLNucleus_hh 1 45 #define G4INCLNucleus_hh 1
47 46
48 #include "G4INCLGlobals.hh" 47 #include "G4INCLGlobals.hh"
49 #include "G4INCLLogger.hh" 48 #include "G4INCLLogger.hh"
50 #include "G4INCLParticle.hh" 49 #include "G4INCLParticle.hh"
51 #include "G4INCLIAvatar.hh" 50 #include "G4INCLIAvatar.hh"
52 #include "G4INCLNucleus.hh" 51 #include "G4INCLNucleus.hh"
53 #include "G4INCLKinematicsUtils.hh" 52 #include "G4INCLKinematicsUtils.hh"
54 #include "G4INCLDecayAvatar.hh" 53 #include "G4INCLDecayAvatar.hh"
55 #include "G4INCLStrangeAbsorbtionChannel.hh" <<
56 #include "G4INCLCluster.hh" 54 #include "G4INCLCluster.hh"
57 #include "G4INCLClusterDecay.hh" 55 #include "G4INCLClusterDecay.hh"
58 #include "G4INCLDeJongSpin.hh" 56 #include "G4INCLDeJongSpin.hh"
59 #include "G4INCLParticleSpecies.hh" << 57 #include "G4INCLNuclearPotentialEnergyIsospinSmooth.hh"
60 #include "G4INCLParticleTable.hh" << 58 #include "G4INCLNuclearPotentialEnergyIsospin.hh"
>> 59 #include "G4INCLNuclearPotentialIsospin.hh"
>> 60 #include "G4INCLNuclearPotentialConstant.hh"
61 #include <iterator> 61 #include <iterator>
62 #include <cstdlib> 62 #include <cstdlib>
63 #include <sstream> 63 #include <sstream>
64 // #include <cassert> 64 // #include <cassert>
65 #include "G4INCLPbarAtrestEntryChannel.hh" <<
66 #include "G4INCLBinaryCollisionAvatar.hh" <<
67 65
68 namespace G4INCL { 66 namespace G4INCL {
69 << 67
70 Nucleus::Nucleus(G4int mass, G4int charge, G << 68 Nucleus::Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius)
71 AnnihilationType AType) //D << 69 : Cluster(charge,mass),
72 : Cluster(charge,mass,strangess,true), << 70 theInitialZ(charge), theInitialA(mass),
73 theInitialZ(charge), theInitialA(mass), t << 71 theNpInitial(0), theNnInitial(0),
74 theNpInitial(0), theNnInitial(0), <<
75 theNpionplusInitial(0), theNpionminusInit <<
76 theNkaonplusInitial(0), theNkaonminusInit <<
77 theNantiprotonInitial(0), <<
78 initialInternalEnergy(0.), 72 initialInternalEnergy(0.),
79 incomingAngularMomentum(0.,0.,0.), incomi 73 incomingAngularMomentum(0.,0.,0.), incomingMomentum(0.,0.,0.),
80 initialCenterOfMass(0.,0.,0.), 74 initialCenterOfMass(0.,0.,0.),
81 remnant(true), 75 remnant(true),
>> 76 blockedDelta(NULL),
82 initialEnergy(0.), 77 initialEnergy(0.),
83 tryCN(false), 78 tryCN(false),
>> 79 forceTransparent(false),
>> 80 projectileZ(0),
>> 81 projectileA(0),
84 theUniverseRadius(universeRadius), 82 theUniverseRadius(universeRadius),
85 isNucleusNucleus(false), 83 isNucleusNucleus(false),
86 theProjectileRemnant(NULL), 84 theProjectileRemnant(NULL),
87 theDensity(NULL), 85 theDensity(NULL),
88 thePotential(NULL), << 86 thePotential(NULL)
89 theAType(AType) <<
90 { 87 {
91 PotentialType potentialType; 88 PotentialType potentialType;
92 G4bool pionPotential; 89 G4bool pionPotential;
93 if(conf) { 90 if(conf) {
94 potentialType = conf->getPotentialType() 91 potentialType = conf->getPotentialType();
95 pionPotential = conf->getPionPotential() 92 pionPotential = conf->getPionPotential();
96 } else { // By default we don't use energy 93 } else { // By default we don't use energy dependent
97 // potential. This is convenient for som 94 // potential. This is convenient for some tests.
98 potentialType = IsospinPotential; 95 potentialType = IsospinPotential;
99 pionPotential = true; 96 pionPotential = true;
100 } 97 }
101 << 98 switch(potentialType) {
102 thePotential = NuclearPotential::createPot << 99 case IsospinEnergySmoothPotential:
>> 100 thePotential = new NuclearPotential::NuclearPotentialEnergyIsospinSmooth(theA, theZ, pionPotential);
>> 101 break;
>> 102 case IsospinEnergyPotential:
>> 103 thePotential = new NuclearPotential::NuclearPotentialEnergyIsospin(theA, theZ, pionPotential);
>> 104 break;
>> 105 case IsospinPotential:
>> 106 thePotential = new NuclearPotential::NuclearPotentialIsospin(theA, theZ, pionPotential);
>> 107 break;
>> 108 case ConstantPotential:
>> 109 thePotential = new NuclearPotential::NuclearPotentialConstant(theA, theZ, pionPotential);
>> 110 break;
>> 111 default:
>> 112 FATAL("Unrecognized potential type at Nucleus creation." << std::endl);
>> 113 break;
>> 114 }
103 115
104 ParticleTable::setProtonSeparationEnergy(t 116 ParticleTable::setProtonSeparationEnergy(thePotential->getSeparationEnergy(Proton));
105 ParticleTable::setNeutronSeparationEnergy( 117 ParticleTable::setNeutronSeparationEnergy(thePotential->getSeparationEnergy(Neutron));
106 118
107 if (theAType==PType) theDensity = NuclearD << 119 theDensity = NuclearDensityFactory::createDensity(theA, theZ);
108 else if (theAType==NType) theDensity = Nuc <<
109 else <<
110 theDensity = NuclearDensityFactory::create <<
111 120
112 theParticleSampler->setPotential(thePotent 121 theParticleSampler->setPotential(thePotential);
113 theParticleSampler->setDensity(theDensity) 122 theParticleSampler->setDensity(theDensity);
114 123
115 if(theUniverseRadius<0) 124 if(theUniverseRadius<0)
116 theUniverseRadius = theDensity->getMaxim 125 theUniverseRadius = theDensity->getMaximumRadius();
117 theStore = new Store(conf); 126 theStore = new Store(conf);
>> 127 toBeUpdated.clear();
118 } 128 }
119 129
120 Nucleus::~Nucleus() { 130 Nucleus::~Nucleus() {
121 delete theStore; 131 delete theStore;
122 deleteProjectileRemnant(); <<
123 /* We don't delete the potential and the d <<
124 * are caching them <<
125 delete thePotential; 132 delete thePotential;
>> 133 /* We don't delete the density here any more -- the Factory is caching them
126 delete theDensity;*/ 134 delete theDensity;*/
127 } 135 }
128 136
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() { 137 void Nucleus::initializeParticles() {
138 // Reset the variables connected with the 138 // Reset the variables connected with the projectile remnant
139 delete theProjectileRemnant; 139 delete theProjectileRemnant;
140 theProjectileRemnant = NULL; 140 theProjectileRemnant = NULL;
141 Cluster::initializeParticles(); <<
142 141
143 for(ParticleIter i=particles.begin(), e=pa << 142 Cluster::initializeParticles();
>> 143 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
144 updatePotentialEnergy(*i); 144 updatePotentialEnergy(*i);
>> 145 theStore->add(*i);
145 } 146 }
146 theStore->add(particles); <<
147 particles.clear(); 147 particles.clear();
148 initialInternalEnergy = computeTotalEnergy 148 initialInternalEnergy = computeTotalEnergy();
149 initialCenterOfMass = thePosition; 149 initialCenterOfMass = thePosition;
150 } 150 }
151 151
152 void Nucleus::applyFinalState(FinalState *fi << 152 std::string Nucleus::dump() {
153 if(!finalstate) // do nothing if no final << 153 std::stringstream ss;
154 return; << 154 ss <<"(list ;; List of participants " << std::endl;
>> 155 ParticleList participants = theStore->getParticipants();
>> 156 for(ParticleIter i = participants.begin(); i != participants.end(); ++i) {
>> 157 ss <<"(make-particle-avatar-map " << std::endl
>> 158 << (*i)->dump()
>> 159 << "(list ;; List of avatars in this particle" << std::endl
>> 160 << ")) ;; Close the list of avatars and the particle-avatar-map" << std::endl;
>> 161 }
>> 162 ss << ")" << std::endl;
>> 163 return ss.str();
>> 164 }
155 165
>> 166 void Nucleus::applyFinalState(FinalState *finalstate) {
>> 167 justCreated.clear();
>> 168 toBeUpdated.clear(); // Clear the list of particles to be updated by the propagation model.
>> 169 blockedDelta = NULL;
156 G4double totalEnergy = 0.0; 170 G4double totalEnergy = 0.0;
157 171
158 FinalStateValidity const validity = finals 172 FinalStateValidity const validity = finalstate->getValidity();
159 if(validity == ValidFS) { 173 if(validity == ValidFS) {
160 174
161 ParticleList const &created = finalstate 175 ParticleList const &created = finalstate->getCreatedParticles();
162 for(ParticleIter iter=created.begin(), e << 176 for(ParticleIter iter = created.begin(); iter != created.end(); ++iter) {
163 theStore->add((*iter)); 177 theStore->add((*iter));
164 if(!(*iter)->isOutOfWell()) { 178 if(!(*iter)->isOutOfWell()) {
165 totalEnergy += (*iter)->getEnergy() 179 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
>> 180 justCreated.push_back((*iter)); // New particle, so we must create avatars for it
166 } 181 }
167 } 182 }
168 183
169 ParticleList const &deleted = finalstate 184 ParticleList const &deleted = finalstate->getDestroyedParticles();
170 for(ParticleIter iter=deleted.begin(), e << 185 for(ParticleIter iter = deleted.begin(); iter != deleted.end(); ++iter) {
171 theStore->particleHasBeenDestroyed(*it << 186 theStore->particleHasBeenDestroyed((*iter)->getID());
172 } 187 }
173 188
174 ParticleList const &modified = finalstat 189 ParticleList const &modified = finalstate->getModifiedParticles();
175 for(ParticleIter iter=modified.begin(), << 190 for(ParticleIter iter = modified.begin(); iter != modified.end(); ++iter) {
176 theStore->particleHasBeenUpdated(*iter << 191 theStore->particleHasBeenUpdated((*iter)->getID());
177 totalEnergy += (*iter)->getEnergy() - 192 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
>> 193 toBeUpdated.push_back((*iter)); // Particle is modified so we have to create new avatars for it.
178 } 194 }
179 195
180 ParticleList const &out = finalstate->ge 196 ParticleList const &out = finalstate->getOutgoingParticles();
181 for(ParticleIter iter=out.begin(), e=out << 197 for(ParticleIter iter = out.begin(); iter != out.end(); ++iter) {
182 if((*iter)->isCluster()) { 198 if((*iter)->isCluster()) {
183 Cluster *clusterOut = dynamic_cast<Clust << 199 Cluster *clusterOut = dynamic_cast<Cluster*>((*iter));
184 // assert(clusterOut); << 200 ParticleList const components = clusterOut->getParticles();
185 #ifdef INCLXX_IN_GEANT4_MODE << 201 for(ParticleIter in = components.begin(); in != components.end(); ++in)
186 if(!clusterOut) << 202 theStore->particleHasBeenEjected((*in)->getID());
187 continue; <<
188 #endif <<
189 ParticleList const &components = clu <<
190 for(ParticleIter in=components.begin <<
191 theStore->particleHasBeenEjected(* <<
192 } else { 203 } else {
193 theStore->particleHasBeenEjected(*it << 204 theStore->particleHasBeenEjected((*iter)->getID());
194 } 205 }
195 totalEnergy += (*iter)->getEnergy(); 206 totalEnergy += (*iter)->getEnergy(); // No potential here because the particle is gone
196 theA -= (*iter)->getA(); 207 theA -= (*iter)->getA();
197 theZ -= (*iter)->getZ(); 208 theZ -= (*iter)->getZ();
198 theS -= (*iter)->getS(); <<
199 theStore->addToOutgoing(*iter); 209 theStore->addToOutgoing(*iter);
200 (*iter)->setEmissionTime(theStore->get << 210 (*iter)->setEmissionTime(theStore->getBook()->getCurrentTime());
201 } 211 }
202 212
203 ParticleList const &entering = finalstat 213 ParticleList const &entering = finalstate->getEnteringParticles();
204 for(ParticleIter iter=entering.begin(), << 214 for(ParticleIter iter = entering.begin(); iter != entering.end(); ++iter) {
205 insertParticle(*iter); 215 insertParticle(*iter);
206 totalEnergy += (*iter)->getEnergy() - 216 totalEnergy += (*iter)->getEnergy() - (*iter)->getPotentialEnergy();
>> 217 toBeUpdated.push_back((*iter)); // Particle is modified so we have to create new avatars for it.
207 } 218 }
208 << 219 } else if(validity == PauliBlockedFS) {
209 // actually perform the removal of the s << 220 blockedDelta = finalstate->getBlockedDelta();
210 theStore->removeScheduledAvatars(); << 221 } else if(validity == ParticleBelowFermiFS) {
211 } else if(validity == ParticleBelowFermiFS << 222 DEBUG("A Particle is entering below the Fermi sea:" << std::endl << finalstate->print() << std::endl);
212 INCL_DEBUG("A Particle is entering below <<
213 tryCN = true; 223 tryCN = true;
214 ParticleList const &entering = finalstat 224 ParticleList const &entering = finalstate->getEnteringParticles();
215 for(ParticleIter iter=entering.begin(), << 225 for(ParticleIter iter = entering.begin(); iter != entering.end(); ++iter) {
>> 226 insertParticle(*iter);
>> 227 }
>> 228 } else if(validity == ParticleBelowZeroFS) {
>> 229 DEBUG("A Particle is entering below zero energy:" << std::endl << finalstate->print() << std::endl);
>> 230 forceTransparent = true;
>> 231 ParticleList const &entering = finalstate->getEnteringParticles();
>> 232 for(ParticleIter iter = entering.begin(); iter != entering.end(); ++iter) {
216 insertParticle(*iter); 233 insertParticle(*iter);
217 } 234 }
218 } 235 }
219 236
220 if(validity==ValidFS && 237 if(validity==ValidFS &&
221 std::abs(totalEnergy - finalstate->get 238 std::abs(totalEnergy - finalstate->getTotalEnergyBeforeInteraction()) > 0.1) {
222 INCL_ERROR("Energy nonconservation! Ener << 239 ERROR("Energy nonconservation! Energy at the beginning of the event = "
223 << finalstate->getTotalEnergyBeforeI 240 << finalstate->getTotalEnergyBeforeInteraction()
224 <<" and after interaction = " 241 <<" and after interaction = "
225 << totalEnergy << '\n' << 242 << totalEnergy << std::endl
226 << finalstate->print()); 243 << finalstate->print());
227 } 244 }
228 } 245 }
229 246
230 void Nucleus::propagateParticles(G4double /* 247 void Nucleus::propagateParticles(G4double /*step*/) {
231 INCL_WARN("Useless Nucleus::propagateParti << 248 WARN("Useless Nucleus::propagateParticles -method called." << std::endl);
232 } 249 }
233 250
234 G4double Nucleus::computeTotalEnergy() const 251 G4double Nucleus::computeTotalEnergy() const {
235 G4double totalEnergy = 0.0; 252 G4double totalEnergy = 0.0;
236 ParticleList const &inside = theStore->get << 253 ParticleList inside = theStore->getParticles();
237 for(ParticleIter p=inside.begin(), e=insid << 254 for(ParticleIter p=inside.begin(); p!=inside.end(); ++p) {
238 if((*p)->isNucleon()) // Ugly: we should 255 if((*p)->isNucleon()) // Ugly: we should calculate everything using total energies!
239 totalEnergy += (*p)->getKineticEnergy( 256 totalEnergy += (*p)->getKineticEnergy() - (*p)->getPotentialEnergy();
240 else if((*p)->isResonance()) 257 else if((*p)->isResonance())
241 totalEnergy += (*p)->getEnergy() - (*p 258 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy() - ParticleTable::effectiveNucleonMass;
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 259 else
250 totalEnergy += (*p)->getEnergy() - (*p 260 totalEnergy += (*p)->getEnergy() - (*p)->getPotentialEnergy();
251 } 261 }
252 <<
253 return totalEnergy; 262 return totalEnergy;
254 } 263 }
255 264
256 void Nucleus::computeRecoilKinematics() { 265 void Nucleus::computeRecoilKinematics() {
257 // If the remnant consists of only one nuc 266 // If the remnant consists of only one nucleon, we need to apply a special
258 // procedure to put it on mass shell. 267 // procedure to put it on mass shell.
259 if(theA==1) { 268 if(theA==1) {
260 emitInsidePions(); 269 emitInsidePions();
261 computeOneNucleonRecoilKinematics(); 270 computeOneNucleonRecoilKinematics();
262 remnant=false; 271 remnant=false;
263 return; 272 return;
264 } 273 }
265 274
266 // Compute the recoil momentum and angular 275 // Compute the recoil momentum and angular momentum
267 theMomentum = incomingMomentum; 276 theMomentum = incomingMomentum;
268 theSpin = incomingAngularMomentum; 277 theSpin = incomingAngularMomentum;
269 278
270 ParticleList const &outgoing = theStore->g << 279 ParticleList outgoing = theStore->getOutgoingParticles();
271 for(ParticleIter p=outgoing.begin(), e=out << 280 for(ParticleIter p=outgoing.begin(); p!=outgoing.end(); ++p)
>> 281 {
272 theMomentum -= (*p)->getMomentum(); 282 theMomentum -= (*p)->getMomentum();
273 theSpin -= (*p)->getAngularMomentum(); 283 theSpin -= (*p)->getAngularMomentum();
274 } 284 }
275 if(theProjectileRemnant) { <<
276 theMomentum -= theProjectileRemnant->get <<
277 theSpin -= theProjectileRemnant->getAngu <<
278 } <<
279 285
280 // Subtract orbital angular momentum 286 // Subtract orbital angular momentum
281 thePosition = computeCenterOfMass(); 287 thePosition = computeCenterOfMass();
282 theSpin -= (thePosition-initialCenterOfMas 288 theSpin -= (thePosition-initialCenterOfMass).vector(theMomentum);
283 289
284 setMass(ParticleTable::getTableMass(theA,t << 290 setMass(ParticleTable::getTableMass(theA,theZ) + theExcitationEnergy);
285 adjustEnergyFromMomentum(); 291 adjustEnergyFromMomentum();
286 remnant=true; 292 remnant=true;
287 } 293 }
288 294
289 ThreeVector Nucleus::computeCenterOfMass() c 295 ThreeVector Nucleus::computeCenterOfMass() const {
290 ThreeVector cm(0.,0.,0.); 296 ThreeVector cm(0.,0.,0.);
291 G4double totalMass = 0.0; 297 G4double totalMass = 0.0;
292 ParticleList const &inside = theStore->get << 298 ParticleList inside = theStore->getParticles();
293 for(ParticleIter p=inside.begin(), e=insid << 299 for(ParticleIter p=inside.begin(); p!=inside.end(); ++p) {
294 const G4double mass = (*p)->getMass(); 300 const G4double mass = (*p)->getMass();
295 cm += (*p)->getPosition() * mass; 301 cm += (*p)->getPosition() * mass;
296 totalMass += mass; 302 totalMass += mass;
297 } 303 }
298 cm /= totalMass; 304 cm /= totalMass;
299 return cm; 305 return cm;
300 } 306 }
301 307
302 G4double Nucleus::computeExcitationEnergy() 308 G4double Nucleus::computeExcitationEnergy() const {
303 const G4double totalEnergy = computeTotalE 309 const G4double totalEnergy = computeTotalEnergy();
304 const G4double separationEnergies = comput 310 const G4double separationEnergies = computeSeparationEnergyBalance();
305 311
306 G4double eSep = 0; << 312 return totalEnergy - initialInternalEnergy - separationEnergies;
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 } 313 }
327 314
328 //thePotential->getSeparationEnergy(Proton) <<
329 <<
330 std::string Nucleus::print() 315 std::string Nucleus::print()
331 { 316 {
332 std::stringstream ss; 317 std::stringstream ss;
333 ss << "Particles in the nucleus:" << '\n' << 318 ss << "Particles in the nucleus:" << std::endl
334 << "Inside:" << '\n'; << 319 << "Participants:" << std::endl;
335 G4int counter = 1; 320 G4int counter = 1;
336 ParticleList const &inside = theStore->get << 321 ParticleList participants = theStore->getParticipants();
337 for(ParticleIter p=inside.begin(), e=insid << 322 for(ParticleIter p = participants.begin(); p != participants.end(); ++p) {
338 ss << "index = " << counter << '\n' << 323 ss << "index = " << counter << std::endl
339 << (*p)->print(); 324 << (*p)->print();
340 counter++; 325 counter++;
341 } 326 }
342 ss <<"Outgoing:" << '\n'; << 327 ss <<"Spectators:" << std::endl;
343 ParticleList const &outgoing = theStore->g << 328 ParticleList spectators = theStore->getSpectators();
344 for(ParticleIter p=outgoing.begin(), e=out << 329 for(ParticleIter p = spectators.begin(); p != spectators.end(); ++p)
>> 330 ss << (*p)->print();
>> 331 ss <<"Outgoing:" << std::endl;
>> 332 ParticleList outgoing = theStore->getOutgoingParticles();
>> 333 for(ParticleIter p = outgoing.begin(); p != outgoing.end(); ++p)
345 ss << (*p)->print(); 334 ss << (*p)->print();
346 335
347 return ss.str(); 336 return ss.str();
348 } 337 }
349 338
350 G4bool Nucleus::decayOutgoingDeltas() { 339 G4bool Nucleus::decayOutgoingDeltas() {
351 ParticleList const &out = theStore->getOut << 340 ParticleList out = theStore->getOutgoingParticles();
352 ParticleList deltas; 341 ParticleList deltas;
353 for(ParticleIter i=out.begin(), e=out.end( << 342 for(ParticleIter i = out.begin(); i != out.end(); ++i) {
354 if((*i)->isDelta()) deltas.push_back((*i 343 if((*i)->isDelta()) deltas.push_back((*i));
355 } 344 }
356 if(deltas.empty()) return false; 345 if(deltas.empty()) return false;
357 346
358 for(ParticleIter i=deltas.begin(), e=delta << 347 for(ParticleIter i = deltas.begin(); i != deltas.end(); ++i) {
359 INCL_DEBUG("Decay outgoing delta particl << 348 DEBUG("Decay outgoing delta particle:" << std::endl
360 << (*i)->print() << '\n'); << 349 << (*i)->print() << std::endl);
361 const ThreeVector beta = -(*i)->boostVec 350 const ThreeVector beta = -(*i)->boostVector();
362 const G4double deltaMass = (*i)->getMass 351 const G4double deltaMass = (*i)->getMass();
363 352
364 // Set the delta momentum to zero and sa 353 // Set the delta momentum to zero and sample the decay in the CM frame.
365 // This makes life simpler if we are usi 354 // This makes life simpler if we are using real particle masses.
366 (*i)->setMomentum(ThreeVector()); 355 (*i)->setMomentum(ThreeVector());
367 (*i)->setEnergy((*i)->getMass()); 356 (*i)->setEnergy((*i)->getMass());
368 357
369 // Use a DecayAvatar 358 // Use a DecayAvatar
370 IAvatar *decay = new DecayAvatar((*i), 0 359 IAvatar *decay = new DecayAvatar((*i), 0.0, NULL);
371 FinalState *fs = decay->getFinalState(); 360 FinalState *fs = decay->getFinalState();
372 Particle * const pion = fs->getCreatedPa 361 Particle * const pion = fs->getCreatedParticles().front();
373 Particle * const nucleon = fs->getModifi 362 Particle * const nucleon = fs->getModifiedParticles().front();
374 363
375 // Adjust the decay momentum if we are u 364 // Adjust the decay momentum if we are using the real masses
376 const G4double decayMomentum = Kinematic 365 const G4double decayMomentum = KinematicsUtils::momentumInCM(deltaMass,
377 nucleon->getTableMass(), 366 nucleon->getTableMass(),
378 pion->getTableMass()); 367 pion->getTableMass());
379 ThreeVector newMomentum = pion->getMomen 368 ThreeVector newMomentum = pion->getMomentum();
380 newMomentum *= decayMomentum / newMoment 369 newMomentum *= decayMomentum / newMomentum.mag();
381 370
382 pion->setTableMass(); 371 pion->setTableMass();
383 pion->setMomentum(newMomentum); 372 pion->setMomentum(newMomentum);
384 pion->adjustEnergyFromMomentum(); 373 pion->adjustEnergyFromMomentum();
385 pion->setEmissionTime(nucleon->getEmissi << 374 pion->setEmissionTime(theStore->getBook()->getCurrentTime());
386 pion->boost(beta); 375 pion->boost(beta);
387 pion->setBiasCollisionVector(nucleon->ge <<
388 376
389 nucleon->setTableMass(); 377 nucleon->setTableMass();
390 nucleon->setMomentum(-newMomentum); 378 nucleon->setMomentum(-newMomentum);
391 nucleon->adjustEnergyFromMomentum(); 379 nucleon->adjustEnergyFromMomentum();
>> 380 nucleon->setEmissionTime(theStore->getBook()->getCurrentTime());
392 nucleon->boost(beta); 381 nucleon->boost(beta);
393 382
394 theStore->addToOutgoing(pion); 383 theStore->addToOutgoing(pion);
395 384
396 delete fs; 385 delete fs;
397 delete decay; 386 delete decay;
398 } 387 }
399 388
400 return true; 389 return true;
401 } 390 }
402 391
403 G4bool Nucleus::decayInsideDeltas() { 392 G4bool Nucleus::decayInsideDeltas() {
404 /* If there is a pion potential, do nothin 393 /* If there is a pion potential, do nothing (deltas will be counted as
405 * excitation energy). 394 * excitation energy).
406 * If, however, the remnant is unphysical 395 * 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 396 * 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 c 397 * end up with Z<0 or Z>A if the remnant contains more pi- than protons or
409 * more pi+ than neutrons, respectively. 398 * more pi+ than neutrons, respectively.
410 */ 399 */
411 const G4bool unphysicalRemnant = (theZ<0 | 400 const G4bool unphysicalRemnant = (theZ<0 || theZ>theA);
412 if(thePotential->hasPionPotential() && !un 401 if(thePotential->hasPionPotential() && !unphysicalRemnant)
413 return false; 402 return false;
414 403
415 // Build a list of deltas (avoid modifying 404 // Build a list of deltas (avoid modifying the list you are iterating on).
416 ParticleList const &inside = theStore->get << 405 ParticleList inside = theStore->getParticles();
417 ParticleList deltas; 406 ParticleList deltas;
418 for(ParticleIter i=inside.begin(), e=insid << 407 for(ParticleIter i = inside.begin(); i != inside.end(); ++i)
419 if((*i)->isDelta()) deltas.push_back((*i 408 if((*i)->isDelta()) deltas.push_back((*i));
420 409
421 // Loop over the deltas, make them decay 410 // Loop over the deltas, make them decay
422 for(ParticleIter i=deltas.begin(), e=delta << 411 for(ParticleIter i = deltas.begin(); i != deltas.end(); ++i) {
423 INCL_DEBUG("Decay inside delta particle: << 412 DEBUG("Decay inside delta particle:" << std::endl
424 << (*i)->print() << '\n'); << 413 << (*i)->print() << std::endl);
425 // Create a forced-decay avatar. Note th 414 // Create a forced-decay avatar. Note the last boolean parameter. Note
426 // also that if the remnant is unphysica 415 // also that if the remnant is unphysical we more or less explicitly give
427 // up energy conservation and CDPP by pa 416 // up energy conservation and CDPP by passing a NULL pointer for the
428 // nucleus. 417 // nucleus.
429 IAvatar *decay; 418 IAvatar *decay;
430 if(unphysicalRemnant) { << 419 if(unphysicalRemnant)
431 INCL_WARN("Forcing delta decay inside <<
432 << ", Z=" << theZ << "). Might le <<
433 << '\n'); <<
434 decay = new DecayAvatar((*i), 0.0, NUL 420 decay = new DecayAvatar((*i), 0.0, NULL, true);
435 } else << 421 else
436 decay = new DecayAvatar((*i), 0.0, thi 422 decay = new DecayAvatar((*i), 0.0, this, true);
437 FinalState *fs = decay->getFinalState(); 423 FinalState *fs = decay->getFinalState();
438 424
439 // The pion can be ejected only if we ma 425 // The pion can be ejected only if we managed to satisfy energy
440 // conservation and if pion emission doe 426 // conservation and if pion emission does not lead to negative excitation
441 // energies. 427 // energies.
442 if(fs->getValidity()==ValidFS) { 428 if(fs->getValidity()==ValidFS) {
443 // Apply the final state to the nucleu 429 // Apply the final state to the nucleus
444 applyFinalState(fs); 430 applyFinalState(fs);
445 } 431 }
446 delete fs; 432 delete fs;
447 delete decay; 433 delete decay;
448 } 434 }
449 435
450 // If the remnant is unphysical, emit all 436 // If the remnant is unphysical, emit all the pions
451 if(unphysicalRemnant) { 437 if(unphysicalRemnant) {
452 INCL_DEBUG("Remnant is unphysical: Z=" < << 438 DEBUG("Remnant is unphysical: Z=" << theZ << ", A=" << theA << std::endl);
453 emitInsidePions(); 439 emitInsidePions();
454 } 440 }
455 441
456 return true; 442 return true;
457 } 443 }
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 444
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() { 445 G4bool Nucleus::decayOutgoingClusters() {
674 ParticleList const &out = theStore->getOut << 446 ParticleList out = theStore->getOutgoingParticles();
675 ParticleList clusters; 447 ParticleList clusters;
676 for(ParticleIter i=out.begin(), e=out.end( << 448 for(ParticleIter i = out.begin(); i != out.end(); ++i) {
677 if((*i)->isCluster()) clusters.push_back 449 if((*i)->isCluster()) clusters.push_back((*i));
678 } 450 }
679 if(clusters.empty()) return false; 451 if(clusters.empty()) return false;
680 452
681 for(ParticleIter i=clusters.begin(), e=clu << 453 for(ParticleIter i = clusters.begin(); i != clusters.end(); ++i) {
682 Cluster *cluster = dynamic_cast<Cluster* 454 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 nee 455 cluster->deleteParticles(); // Don't need them
689 ParticleList decayProducts = ClusterDeca 456 ParticleList decayProducts = ClusterDecay::decay(cluster);
690 for(ParticleIter j=decayProducts.begin() << 457 for(ParticleIter j = decayProducts.begin(); j!=decayProducts.end(); ++j)
691 (*j)->setBiasCollisionVector(cluster-> <<
692 theStore->addToOutgoing(*j); 458 theStore->addToOutgoing(*j);
693 } <<
694 } 459 }
695 return true; 460 return true;
696 } 461 }
697 462
698 G4bool Nucleus::decayMe() { 463 G4bool Nucleus::decayMe() {
699 // Do the phase-space decay only if Z=0 or << 464 // Do the phase-space decay only if Z=0 or Z=A
700 if(theA<=1 || (theZ!=0 && (theA+theS)!=the << 465 if(theA<=1 || (theZ!=0 && theA!=theZ))
701 return false; 466 return false;
702 467
703 ParticleList decayProducts = ClusterDecay: 468 ParticleList decayProducts = ClusterDecay::decay(this);
704 for(ParticleIter j=decayProducts.begin(), << 469 for(ParticleIter j = decayProducts.begin(); j!=decayProducts.end(); ++j)
705 (*j)->setBiasCollisionVector(this->getBi <<
706 theStore->addToOutgoing(*j); 470 theStore->addToOutgoing(*j);
707 } << 471
708 <<
709 return true; 472 return true;
710 } 473 }
711 474
712 void Nucleus::emitInsidePions() { 475 void Nucleus::emitInsidePions() {
713 /* Forcing emissions of all pions in the n 476 /* Forcing emissions of all pions in the nucleus. This probably violates
714 * energy conservation (although the compu 477 * energy conservation (although the computation of the recoil kinematics
715 * might sweep this under the carpet). 478 * might sweep this under the carpet).
716 */ 479 */
717 INCL_WARN("Forcing emissions of all pions << 480 WARN("Forcing emissions of all pions in the nucleus." << std::endl);
718 481
719 // Emit the pions with this kinetic energy 482 // Emit the pions with this kinetic energy
720 const G4double tinyPionEnergy = 0.1; // Me 483 const G4double tinyPionEnergy = 0.1; // MeV
721 484
722 // Push out the emitted pions 485 // Push out the emitted pions
723 ParticleList const &inside = theStore->get << 486 ParticleList inside = theStore->getParticles();
724 ParticleList toEject; << 487 for(ParticleIter i = inside.begin(); i != inside.end(); ++i) {
725 for(ParticleIter i=inside.begin(), e=insid <<
726 if((*i)->isPion()) { 488 if((*i)->isPion()) {
727 Particle * const thePion = *i; << 489 (*i)->setEmissionTime(theStore->getBook()->getCurrentTime());
728 INCL_DEBUG("Forcing emission of the fo <<
729 << thePion->print() << '\n' <<
730 thePion->setEmissionTime(theStore->get <<
731 // Correction for real masses 490 // Correction for real masses
732 const G4double theQValueCorrection = t << 491 const G4double theQValueCorrection = (*i)->getEmissionQValueCorrection(theA,theZ);
733 const G4double kineticEnergyOutside = << 492 const G4double kineticEnergyOutside = (*i)->getKineticEnergy() - (*i)->getPotentialEnergy() + theQValueCorrection;
734 thePion->setTableMass(); << 493 (*i)->setTableMass();
735 if(kineticEnergyOutside > 0.0) 494 if(kineticEnergyOutside > 0.0)
736 thePion->setEnergy(thePion->getMass( << 495 (*i)->setEnergy((*i)->getMass()+kineticEnergyOutside);
737 else 496 else
738 thePion->setEnergy(thePion->getMass( << 497 (*i)->setEnergy((*i)->getMass()+tinyPionEnergy);
739 thePion->adjustMomentumFromEnergy(); << 498 (*i)->adjustMomentumFromEnergy();
740 thePion->setPotentialEnergy(0.); << 499 (*i)->setPotentialEnergy(0.);
741 theZ -= thePion->getZ(); << 500 theZ -= (*i)->getZ();
742 toEject.push_back(thePion); << 501 theStore->particleHasBeenEjected((*i)->getID());
>> 502 theStore->addToOutgoing(*i);
743 } 503 }
744 } 504 }
745 for(ParticleIter i=toEject.begin(), e=toEj <<
746 theStore->particleHasBeenEjected(*i); <<
747 theStore->addToOutgoing(*i); <<
748 (*i)->setParticleBias(Particle::getTotal <<
749 } <<
750 } 505 }
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 506
760 // Emit the strange particles with this ki << 507 G4bool Nucleus::isEventTransparent() const {
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 508
795 G4int Nucleus::emitInsideLambda() { << 509 // Forced transparent
796 /* Forcing emissions of all Lambda in the << 510 if(forceTransparent)
797 * This probably violates energy conservat << 511 return true;
798 * (although the computation of the recoil <<
799 * might sweep this under the carpet). <<
800 */ <<
801 INCL_DEBUG("Forcing emissions of all Lambd <<
802 512
803 // Emit the Lambda with this kinetic energ << 513 ParticleList const &pL = theStore->getOutgoingParticles();
804 const G4double tinyEnergy = 0.1; // MeV << 514 G4int outZ = 0, outA = 0;
805 515
806 // Push out the emitted Lambda << 516 // If any of the particles has undergone a collision, the event is not a
807 ParticleList const &inside = theStore->get << 517 // transparent.
808 ParticleList toEject; << 518 for(ParticleIter p = pL.begin(); p != pL.end(); ++p ) {
809 for(ParticleIter i=inside.begin(), e=insid << 519 if( (*p)->getNumberOfCollisions() != 0 ) return false;
810 if((*i)->isLambda()) { << 520 if( (*p)->getNumberOfDecays() != 0 ) return false;
811 Particle * const theLambda = *i; << 521 outZ += (*p)->getZ();
812 INCL_DEBUG("Forcing emission of the fo << 522 outA += (*p)->getA();
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 } 523 }
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 524
849 // Push out the emitted kaon << 525 // Add the geometrical spectators to the Z and A count
850 ParticleList const &inside = theStore->get << 526 if(theProjectileRemnant) {
851 ParticleList toEject; << 527 outZ += theProjectileRemnant->getZ();
852 for(ParticleIter i=inside.begin(), e=insid << 528 outA += theProjectileRemnant->getA();
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 } 529 }
878 theNKaon -= 1; <<
879 return toEject.size() != 0; <<
880 } <<
881 <<
882 G4bool Nucleus::isEventTransparent() const { <<
883 530
884 Book const &theBook = theStore->getBook(); << 531 if(outZ!=projectileZ || outA!=projectileA) return false;
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 532
891 return false; << 533 return true;
892 534
893 } 535 }
894 536
895 void Nucleus::computeOneNucleonRecoilKinemat 537 void Nucleus::computeOneNucleonRecoilKinematics() {
896 // We should be here only if the nucleus c 538 // We should be here only if the nucleus contains only one nucleon
897 // assert(theStore->getParticles().size()==1); 539 // assert(theStore->getParticles().size()==1);
898 540
>> 541 ERROR("Computing one-nucleon recoil kinematics. We should never be here nowadays, cascade should stop earlier than this." << std::endl);
>> 542
899 // No excitation energy! 543 // No excitation energy!
900 theExcitationEnergy = 0.0; 544 theExcitationEnergy = 0.0;
901 545
902 // Move the nucleon to the outgoing list 546 // Move the nucleon to the outgoing list
903 Particle *remN = theStore->getParticles(). 547 Particle *remN = theStore->getParticles().front();
904 theA -= remN->getA(); 548 theA -= remN->getA();
905 theZ -= remN->getZ(); 549 theZ -= remN->getZ();
906 theS -= remN->getS(); << 550 theStore->particleHasBeenEjected(remN->getID());
907 theStore->particleHasBeenEjected(remN); <<
908 theStore->addToOutgoing(remN); 551 theStore->addToOutgoing(remN);
909 remN->setEmissionTime(theStore->getBook(). << 552 remN->setEmissionTime(theStore->getBook()->getCurrentTime());
910 553
911 // Treat the special case of a remaining d 554 // Treat the special case of a remaining delta
912 if(remN->isDelta()) { 555 if(remN->isDelta()) {
913 IAvatar *decay = new DecayAvatar(remN, 0 556 IAvatar *decay = new DecayAvatar(remN, 0.0, NULL);
914 FinalState *fs = decay->getFinalState(); 557 FinalState *fs = decay->getFinalState();
915 // Eject the pion 558 // Eject the pion
916 ParticleList const &created = fs->getCre << 559 ParticleList created = fs->getCreatedParticles();
917 for(ParticleIter j=created.begin(), e=cr << 560 for(ParticleIter j = created.begin(); j != created.end(); ++j)
918 theStore->addToOutgoing(*j); 561 theStore->addToOutgoing(*j);
919 delete fs; 562 delete fs;
920 delete decay; 563 delete decay;
921 } 564 }
922 565
923 // Do different things depending on how ma 566 // Do different things depending on how many outgoing particles we have
924 ParticleList const &outgoing = theStore->g << 567 ParticleList outgoing = theStore->getOutgoingParticles();
925 if(outgoing.size() == 2) { 568 if(outgoing.size() == 2) {
926 569
927 INCL_DEBUG("Two particles in the outgoin << 570 DEBUG("Two particles in the outgoing channel, applying exact two-body kinematics" << std::endl);
928 571
929 // Can apply exact 2-body kinematics her 572 // Can apply exact 2-body kinematics here. Keep the CM emission angle of
930 // the first particle. 573 // the first particle.
931 Particle *p1 = outgoing.front(), *p2 = o 574 Particle *p1 = outgoing.front(), *p2 = outgoing.back();
932 const ThreeVector aBoostVector = incomin 575 const ThreeVector aBoostVector = incomingMomentum / initialEnergy;
933 // Boost to the initial CM 576 // Boost to the initial CM
934 p1->boost(aBoostVector); 577 p1->boost(aBoostVector);
935 const G4double sqrts = std::sqrt(initial 578 const G4double sqrts = std::sqrt(initialEnergy*initialEnergy - incomingMomentum.mag2());
936 const G4double pcm = KinematicsUtils::mo 579 const G4double pcm = KinematicsUtils::momentumInCM(sqrts, p1->getMass(), p2->getMass());
937 const G4double scale = pcm/(p1->getMomen 580 const G4double scale = pcm/(p1->getMomentum().mag());
938 // Reset the momenta 581 // Reset the momenta
939 p1->setMomentum(p1->getMomentum()*scale) 582 p1->setMomentum(p1->getMomentum()*scale);
940 p2->setMomentum(-p1->getMomentum()); 583 p2->setMomentum(-p1->getMomentum());
941 p1->adjustEnergyFromMomentum(); 584 p1->adjustEnergyFromMomentum();
942 p2->adjustEnergyFromMomentum(); 585 p2->adjustEnergyFromMomentum();
943 // Unboost 586 // Unboost
944 p1->boost(-aBoostVector); 587 p1->boost(-aBoostVector);
945 p2->boost(-aBoostVector); 588 p2->boost(-aBoostVector);
946 589
947 } else { 590 } else {
948 591
949 INCL_DEBUG("Trying to adjust final-state << 592 DEBUG("Trying to adjust final-state momenta to achieve energy and momentum conservation" << std::endl);
950 593
951 const G4int maxIterations=8; 594 const G4int maxIterations=8;
952 G4double totalEnergy, energyScale; 595 G4double totalEnergy, energyScale;
953 G4double val=1.E+100, oldVal=1.E+100, ol 596 G4double val=1.E+100, oldVal=1.E+100, oldOldVal=1.E+100, oldOldOldVal;
954 ThreeVector totalMomentum, deltaP; 597 ThreeVector totalMomentum, deltaP;
955 std::vector<ThreeVector> minMomenta; // 598 std::vector<ThreeVector> minMomenta; // use it to store the particle momenta that minimize the merit function
956 599
957 // Reserve the vector size 600 // Reserve the vector size
958 minMomenta.reserve(outgoing.size()); 601 minMomenta.reserve(outgoing.size());
959 602
960 // Compute the initial total momentum 603 // Compute the initial total momentum
961 totalMomentum.setX(0.0); 604 totalMomentum.setX(0.0);
962 totalMomentum.setY(0.0); 605 totalMomentum.setY(0.0);
963 totalMomentum.setZ(0.0); 606 totalMomentum.setZ(0.0);
964 for(ParticleIter i=outgoing.begin(), e=o << 607 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
965 totalMomentum += (*i)->getMomentum(); 608 totalMomentum += (*i)->getMomentum();
966 609
967 // Compute the initial total energy 610 // Compute the initial total energy
968 totalEnergy = 0.0; 611 totalEnergy = 0.0;
969 for(ParticleIter i=outgoing.begin(), e=o << 612 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
970 totalEnergy += (*i)->getEnergy(); 613 totalEnergy += (*i)->getEnergy();
971 614
972 // Iterative algorithm starts here: 615 // Iterative algorithm starts here:
973 for(G4int iterations=0; iterations < max 616 for(G4int iterations=0; iterations < maxIterations; ++iterations) {
974 617
975 // Save the old merit-function values 618 // Save the old merit-function values
976 oldOldOldVal = oldOldVal; 619 oldOldOldVal = oldOldVal;
977 oldOldVal = oldVal; 620 oldOldVal = oldVal;
978 oldVal = val; 621 oldVal = val;
979 622
980 if(iterations%2 == 0) { 623 if(iterations%2 == 0) {
981 INCL_DEBUG("Momentum step" << '\n'); << 624 DEBUG("Momentum step" << std::endl);
982 // Momentum step: modify all the par 625 // Momentum step: modify all the particle momenta
983 deltaP = incomingMomentum - totalMom 626 deltaP = incomingMomentum - totalMomentum;
984 G4double pOldTot = 0.0; 627 G4double pOldTot = 0.0;
985 for(ParticleIter i=outgoing.begin(), << 628 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
986 pOldTot += (*i)->getMomentum().mag 629 pOldTot += (*i)->getMomentum().mag();
987 for(ParticleIter i=outgoing.begin(), << 630 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
988 const ThreeVector mom = (*i)->getM 631 const ThreeVector mom = (*i)->getMomentum();
989 (*i)->setMomentum(mom + deltaP*mom 632 (*i)->setMomentum(mom + deltaP*mom.mag()/pOldTot);
990 (*i)->adjustEnergyFromMomentum(); 633 (*i)->adjustEnergyFromMomentum();
991 } 634 }
992 } else { 635 } else {
993 INCL_DEBUG("Energy step" << '\n'); << 636 DEBUG("Energy step" << std::endl);
994 // Energy step: modify all the parti 637 // Energy step: modify all the particle momenta
995 energyScale = initialEnergy/totalEne 638 energyScale = initialEnergy/totalEnergy;
996 for(ParticleIter i=outgoing.begin(), << 639 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
997 const ThreeVector mom = (*i)->getM 640 const ThreeVector mom = (*i)->getMomentum();
998 G4double pScale = ((*i)->getEnergy 641 G4double pScale = ((*i)->getEnergy()*energyScale - std::pow((*i)->getMass(),2))/mom.mag2();
999 if(pScale>0) { 642 if(pScale>0) {
1000 (*i)->setEnergy((*i)->getEnergy 643 (*i)->setEnergy((*i)->getEnergy()*energyScale);
1001 (*i)->adjustMomentumFromEnergy( 644 (*i)->adjustMomentumFromEnergy();
1002 } 645 }
1003 } 646 }
1004 } 647 }
1005 648
1006 // Compute the current total momentum 649 // Compute the current total momentum and energy
1007 totalMomentum.setX(0.0); 650 totalMomentum.setX(0.0);
1008 totalMomentum.setY(0.0); 651 totalMomentum.setY(0.0);
1009 totalMomentum.setZ(0.0); 652 totalMomentum.setZ(0.0);
1010 totalEnergy = 0.0; 653 totalEnergy = 0.0;
1011 for(ParticleIter i=outgoing.begin(), << 654 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i) {
1012 totalMomentum += (*i)->getMomentum( 655 totalMomentum += (*i)->getMomentum();
1013 totalEnergy += (*i)->getEnergy(); 656 totalEnergy += (*i)->getEnergy();
1014 } 657 }
1015 658
1016 // Merit factor 659 // Merit factor
1017 val = std::pow(totalEnergy - initialE 660 val = std::pow(totalEnergy - initialEnergy,2) +
1018 0.25*(totalMomentum - incomingMomen 661 0.25*(totalMomentum - incomingMomentum).mag2();
1019 INCL_DEBUG("Merit function: val=" << << 662 DEBUG("Merit function: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << std::endl);
1020 663
1021 // Store the minimum 664 // Store the minimum
1022 if(val < oldVal) { 665 if(val < oldVal) {
1023 INCL_DEBUG("New minimum found, stor << 666 DEBUG("New minimum found, storing the particle momenta" << std::endl);
1024 minMomenta.clear(); 667 minMomenta.clear();
1025 for(ParticleIter i=outgoing.begin() << 668 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i)
1026 minMomenta.push_back((*i)->getMom 669 minMomenta.push_back((*i)->getMomentum());
1027 } 670 }
1028 671
1029 // Stop the algorithm if the search d 672 // Stop the algorithm if the search diverges
1030 if(val > oldOldVal && oldVal > oldOld 673 if(val > oldOldVal && oldVal > oldOldOldVal) {
1031 INCL_DEBUG("Search is diverging, br << 674 DEBUG("Search is diverging, breaking out of the iteration loop: val=" << val << ", oldVal=" << oldVal << ", oldOldVal=" << oldOldVal << ", oldOldOldVal=" << oldOldOldVal << std::endl);
1032 break; 675 break;
1033 } 676 }
1034 } 677 }
1035 678
1036 // We should have made at least one suc 679 // We should have made at least one successful iteration here
1037 // assert(minMomenta.size()==outgoing.size()) 680 // assert(minMomenta.size()==outgoing.size());
1038 681
1039 // Apply the optimal momenta 682 // Apply the optimal momenta
1040 INCL_DEBUG("Applying the solution" << ' << 683 DEBUG("Applying the solution" << std::endl);
1041 std::vector<ThreeVector>::const_iterato 684 std::vector<ThreeVector>::const_iterator v = minMomenta.begin();
1042 for(ParticleIter i=outgoing.begin(), e= << 685 for(ParticleIter i=outgoing.begin(); i!=outgoing.end(); ++i, ++v) {
1043 (*i)->setMomentum(*v); 686 (*i)->setMomentum(*v);
1044 (*i)->adjustEnergyFromMomentum(); 687 (*i)->adjustEnergyFromMomentum();
1045 INCL_DATABLOCK((*i)->print()); << 688 DATABLOCK((*i)->print());
1046 } 689 }
1047 690
1048 } 691 }
1049 692
1050 } 693 }
1051 694
1052 void Nucleus::fillEventInfo(EventInfo *even 695 void Nucleus::fillEventInfo(EventInfo *eventInfo) {
1053 eventInfo->nParticles = 0; 696 eventInfo->nParticles = 0;
1054 G4bool isNucleonAbsorption = false; 697 G4bool isNucleonAbsorption = false;
>> 698
1055 G4bool isPionAbsorption = false; 699 G4bool isPionAbsorption = false;
1056 // It is possible to have pion absorption 700 // It is possible to have pion absorption event only if the
1057 // projectile is pion. 701 // projectile is pion.
1058 if(eventInfo->projectileType == PiPlus || 702 if(eventInfo->projectileType == PiPlus ||
1059 eventInfo->projectileType == PiMinus | 703 eventInfo->projectileType == PiMinus ||
1060 eventInfo->projectileType == PiZero) { 704 eventInfo->projectileType == PiZero) {
1061 isPionAbsorption = true; 705 isPionAbsorption = true;
1062 } 706 }
1063 707
1064 // Forced CN 708 // Forced CN
1065 eventInfo->forcedCompoundNucleus = tryCN; 709 eventInfo->forcedCompoundNucleus = tryCN;
1066 710
1067 // Outgoing particles 711 // Outgoing particles
1068 ParticleList const &outgoingParticles = g << 712 ParticleList outgoingParticles = getStore()->getOutgoingParticles();
1069 713
1070 // Check if we have a nucleon absorption 714 // Check if we have a nucleon absorption event: nucleon projectile
1071 // and no ejected particles. 715 // and no ejected particles.
1072 if(outgoingParticles.size() == 0 && 716 if(outgoingParticles.size() == 0 &&
1073 (eventInfo->projectileType == Proton | 717 (eventInfo->projectileType == Proton ||
1074 eventInfo->projectileType == Neutron)) { << 718 eventInfo->projectileType == Neutron)) {
1075 isNucleonAbsorption = true; 719 isNucleonAbsorption = true;
1076 } 720 }
1077 721
1078 // Reset the remnant counter 722 // Reset the remnant counter
1079 eventInfo->nRemnants = 0; 723 eventInfo->nRemnants = 0;
1080 eventInfo->history.clear(); 724 eventInfo->history.clear();
1081 725
1082 for(ParticleIter i=outgoingParticles.begi << 726 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i ) {
>> 727 // If the particle is a cluster and has excitation energy, treat it as a cluster
>> 728 if((*i)->isCluster()) {
>> 729 Cluster const * const c = dynamic_cast<Cluster *>(*i);
>> 730 // assert(c);
>> 731 #ifdef INCLXX_IN_GEANT4_MODE
>> 732 if(!c)
>> 733 continue;
>> 734 #endif
>> 735 const G4double eStar = c->getExcitationEnergy();
>> 736 if(std::abs(eStar)>1E-10) {
>> 737 if(eStar<0.) {
>> 738 WARN("Negative excitation energy in outgoing cluster! EStar = " << eStar << std::endl);
>> 739 }
>> 740 eventInfo->ARem[eventInfo->nRemnants] = c->getA();
>> 741 eventInfo->ZRem[eventInfo->nRemnants] = c->getZ();
>> 742 eventInfo->EStarRem[eventInfo->nRemnants] = eStar;
>> 743 ThreeVector remnantSpin = c->getSpin();
>> 744 Float_t remnantSpinMag;
>> 745 if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
>> 746 remnantSpinMag = (G4int) (remnantSpin.mag()/PhysicalConstants::hc + 0.5);
>> 747 } else { // odd-A nucleus
>> 748 remnantSpinMag = ((G4int) (remnantSpin.mag()/PhysicalConstants::hc)) + 0.5;
>> 749 }
>> 750 remnantSpin *= remnantSpinMag/remnantSpin.mag();
>> 751 eventInfo->JRem[eventInfo->nRemnants] = remnantSpinMag;
>> 752 eventInfo->jxRem[eventInfo->nRemnants] = remnantSpin.getX();
>> 753 eventInfo->jyRem[eventInfo->nRemnants] = remnantSpin.getY();
>> 754 eventInfo->jzRem[eventInfo->nRemnants] = remnantSpin.getZ();
>> 755 eventInfo->EKinRem[eventInfo->nRemnants] = c->getKineticEnergy();
>> 756 ThreeVector mom = c->getMomentum();
>> 757 eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
>> 758 eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
>> 759 eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
>> 760 eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
>> 761 eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
>> 762 eventInfo->nRemnants++;
>> 763 continue; // don't add it as a particle
>> 764 }
>> 765 }
>> 766
1083 // We have a pion absorption event only 767 // We have a pion absorption event only if the projectile is
1084 // pion and there are no ejected pions. 768 // pion and there are no ejected pions.
1085 if(isPionAbsorption) { 769 if(isPionAbsorption) {
1086 if((*i)->isPion()) { 770 if((*i)->isPion()) {
1087 isPionAbsorption = false; 771 isPionAbsorption = false;
1088 } 772 }
1089 } 773 }
1090 << 774 eventInfo->A[eventInfo->nParticles] = (*i)->getA();
1091 eventInfo->ParticleBias[eventInfo->nPar << 775 eventInfo->Z[eventInfo->nParticles] = (*i)->getZ();
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 776 eventInfo->emissionTime[eventInfo->nParticles] = (*i)->getEmissionTime();
1103 eventInfo->EKin[eventInfo->nParticles] 777 eventInfo->EKin[eventInfo->nParticles] = (*i)->getKineticEnergy();
1104 ThreeVector mom = (*i)->getMomentum(); 778 ThreeVector mom = (*i)->getMomentum();
1105 eventInfo->px[eventInfo->nParticles] = 779 eventInfo->px[eventInfo->nParticles] = mom.getX();
1106 eventInfo->py[eventInfo->nParticles] = 780 eventInfo->py[eventInfo->nParticles] = mom.getY();
1107 eventInfo->pz[eventInfo->nParticles] = 781 eventInfo->pz[eventInfo->nParticles] = mom.getZ();
1108 eventInfo->theta[eventInfo->nParticles] 782 eventInfo->theta[eventInfo->nParticles] = Math::toDegrees(mom.theta());
1109 eventInfo->phi[eventInfo->nParticles] = 783 eventInfo->phi[eventInfo->nParticles] = Math::toDegrees(mom.phi());
1110 eventInfo->origin[eventInfo->nParticles 784 eventInfo->origin[eventInfo->nParticles] = -1;
1111 #ifdef INCLXX_IN_GEANT4_MODE <<
1112 eventInfo->parentResonancePDGCode[event <<
1113 eventInfo->parentResonanceID[eventInfo- <<
1114 #endif <<
1115 eventInfo->history.push_back(""); 785 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++; 786 eventInfo->nParticles++;
1125 } 787 }
1126 eventInfo->nucleonAbsorption = isNucleonA 788 eventInfo->nucleonAbsorption = isNucleonAbsorption;
1127 eventInfo->pionAbsorption = isPionAbsorpt 789 eventInfo->pionAbsorption = isPionAbsorption;
1128 eventInfo->nCascadeParticles = eventInfo- 790 eventInfo->nCascadeParticles = eventInfo->nParticles;
1129 791
1130 // Projectile-like remnant characteristic << 792 // Remnant characteristics
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()) { 793 if(hasRemnant()) {
1169 #ifdef INCLXX_IN_GEANT4_MODE << 794 eventInfo->ARem[eventInfo->nRemnants] = getA();
1170 eventInfo->ARem[eventInfo->nRemnants] = << 795 eventInfo->ZRem[eventInfo->nRemnants] = getZ();
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 796 eventInfo->EStarRem[eventInfo->nRemnants] = getExcitationEnergy();
1179 if(eventInfo->EStarRem[eventInfo->nRemn 797 if(eventInfo->EStarRem[eventInfo->nRemnants]<0.) {
1180 INCL_WARN("Negative excitation energy in << 798 WARN("Negative excitation energy! EStarRem = " << eventInfo->EStarRem[eventInfo->nRemnants] << std::endl);
1181 } 799 }
1182 const ThreeVector &spin = getSpin(); <<
1183 if(eventInfo->ARem[eventInfo->nRemnants 800 if(eventInfo->ARem[eventInfo->nRemnants]%2==0) { // even-A nucleus
1184 eventInfo->JRem[eventInfo->nRemnants] = ( << 801 eventInfo->JRem[eventInfo->nRemnants] = (G4int) (getSpin().mag()/PhysicalConstants::hc + 0.5);
1185 } else { // odd-A nucleus 802 } else { // odd-A nucleus
1186 eventInfo->JRem[eventInfo->nRemnants] = ( << 803 eventInfo->JRem[eventInfo->nRemnants] = ((G4int) (getSpin().mag()/PhysicalConstants::hc)) + 0.5;
1187 } 804 }
1188 eventInfo->EKinRem[eventInfo->nRemnants 805 eventInfo->EKinRem[eventInfo->nRemnants] = getKineticEnergy();
1189 const ThreeVector &mom = getMomentum(); << 806 ThreeVector mom = getMomentum();
1190 eventInfo->pxRem[eventInfo->nRemnants] 807 eventInfo->pxRem[eventInfo->nRemnants] = mom.getX();
1191 eventInfo->pyRem[eventInfo->nRemnants] 808 eventInfo->pyRem[eventInfo->nRemnants] = mom.getY();
1192 eventInfo->pzRem[eventInfo->nRemnants] 809 eventInfo->pzRem[eventInfo->nRemnants] = mom.getZ();
1193 eventInfo->jxRem[eventInfo->nRemnants] <<
1194 eventInfo->jyRem[eventInfo->nRemnants] <<
1195 eventInfo->jzRem[eventInfo->nRemnants] <<
1196 eventInfo->thetaRem[eventInfo->nRemnant 810 eventInfo->thetaRem[eventInfo->nRemnants] = Math::toDegrees(mom.theta());
1197 eventInfo->phiRem[eventInfo->nRemnants] 811 eventInfo->phiRem[eventInfo->nRemnants] = Math::toDegrees(mom.phi());
1198 eventInfo->nRemnants++; 812 eventInfo->nRemnants++;
1199 } 813 }
1200 814
1201 // Global counters, flags, etc. 815 // Global counters, flags, etc.
1202 Book const &theBook = theStore->getBook() << 816 eventInfo->nCollisions = getStore()->getBook()->getAcceptedCollisions();
1203 eventInfo->nCollisions = theBook.getAccep << 817 eventInfo->nBlockedCollisions = getStore()->getBook()->getBlockedCollisions();
1204 eventInfo->nBlockedCollisions = theBook.g << 818 eventInfo->nDecays = getStore()->getBook()->getAcceptedDecays();
1205 eventInfo->nDecays = theBook.getAcceptedD << 819 eventInfo->nBlockedDecays = getStore()->getBook()->getBlockedDecays();
1206 eventInfo->nBlockedDecays = theBook.getBl << 820 eventInfo->firstCollisionTime = getStore()->getBook()->getFirstCollisionTime();
1207 eventInfo->firstCollisionTime = theBook.g << 821 eventInfo->firstCollisionXSec = getStore()->getBook()->getFirstCollisionXSec();
1208 eventInfo->firstCollisionXSec = theBook.g << 822 eventInfo->nReflectionAvatars = getStore()->getBook()->getAvatars(SurfaceAvatarType);
1209 eventInfo->firstCollisionSpectatorPositio << 823 eventInfo->nCollisionAvatars = getStore()->getBook()->getAvatars(CollisionAvatarType);
1210 eventInfo->firstCollisionSpectatorMomentu << 824 eventInfo->nDecayAvatars = getStore()->getBook()->getAvatars(DecayAvatarType);
1211 eventInfo->firstCollisionIsElastic = theB <<
1212 eventInfo->nReflectionAvatars = theBook.g <<
1213 eventInfo->nCollisionAvatars = theBook.ge <<
1214 eventInfo->nDecayAvatars = theBook.getAva <<
1215 eventInfo->nEnergyViolationInteraction = <<
1216 } 825 }
1217 826
1218 <<
1219 <<
1220 Nucleus::ConservationBalance Nucleus::getCo 827 Nucleus::ConservationBalance Nucleus::getConservationBalance(const EventInfo &theEventInfo, const G4bool afterRecoil) const {
1221 ConservationBalance theBalance; 828 ConservationBalance theBalance;
1222 // Initialise balance variables with the 829 // Initialise balance variables with the incoming values
1223 INCL_DEBUG("theEventInfo " << theEventInf <<
1224 theBalance.Z = theEventInfo.Zp + theEvent 830 theBalance.Z = theEventInfo.Zp + theEventInfo.Zt;
1225 theBalance.A = theEventInfo.Ap + theEvent 831 theBalance.A = theEventInfo.Ap + theEventInfo.At;
1226 theBalance.S = theEventInfo.Sp + theEvent << 832
1227 INCL_DEBUG("theBalance Z and A " << theBa <<
1228 theBalance.energy = getInitialEnergy(); 833 theBalance.energy = getInitialEnergy();
1229 theBalance.momentum = getIncomingMomentum 834 theBalance.momentum = getIncomingMomentum();
1230 835
1231 // Process outgoing particles 836 // Process outgoing particles
1232 ParticleList const &outgoingParticles = t << 837 ParticleList outgoingParticles = theStore->getOutgoingParticles();
1233 for(ParticleIter i=outgoingParticles.begi << 838 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i ) {
1234 theBalance.Z -= (*i)->getZ(); 839 theBalance.Z -= (*i)->getZ();
1235 theBalance.A -= (*i)->getA(); 840 theBalance.A -= (*i)->getA();
1236 theBalance.S -= (*i)->getS(); <<
1237 // For outgoing clusters, the total ene 841 // For outgoing clusters, the total energy automatically includes the
1238 // excitation energy 842 // excitation energy
1239 theBalance.energy -= (*i)->getEnergy(); 843 theBalance.energy -= (*i)->getEnergy(); // Note that outgoing particles should have the real mass
1240 theBalance.momentum -= (*i)->getMomentu 844 theBalance.momentum -= (*i)->getMomentum();
1241 } 845 }
1242 846
1243 // Projectile-like remnant contribution, << 847 // Remnant contribution, if present
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()) { 848 if(hasRemnant()) {
1256 theBalance.Z -= getZ(); 849 theBalance.Z -= getZ();
1257 theBalance.A -= getA(); 850 theBalance.A -= getA();
1258 theBalance.S -= getS(); << 851 theBalance.energy -= ParticleTable::getTableMass(getA(),getZ()) +
1259 theBalance.energy -= ParticleTable::get <<
1260 getExcitationEnergy(); 852 getExcitationEnergy();
1261 if(afterRecoil) 853 if(afterRecoil)
1262 theBalance.energy -= getKineticEnergy 854 theBalance.energy -= getKineticEnergy();
1263 theBalance.momentum -= getMomentum(); 855 theBalance.momentum -= getMomentum();
1264 } 856 }
1265 857
1266 return theBalance; 858 return theBalance;
1267 } 859 }
1268 860
1269 void Nucleus::useFusionKinematics() { 861 void Nucleus::useFusionKinematics() {
1270 setEnergy(initialEnergy); 862 setEnergy(initialEnergy);
1271 setMomentum(incomingMomentum); 863 setMomentum(incomingMomentum);
1272 setSpin(incomingAngularMomentum); 864 setSpin(incomingAngularMomentum);
1273 theExcitationEnergy = std::sqrt(theEnergy 865 theExcitationEnergy = std::sqrt(theEnergy*theEnergy-theMomentum.mag2()) - getTableMass();
1274 setMass(getTableMass() + theExcitationEne 866 setMass(getTableMass() + theExcitationEnergy);
1275 } 867 }
1276 868
1277 void Nucleus::finalizeProjectileRemnant(con 869 void Nucleus::finalizeProjectileRemnant(const G4double anEmissionTime) {
1278 // Deal with the projectile remnant 870 // Deal with the projectile remnant
1279 const G4int prA = theProjectileRemnant->g << 871 if(theProjectileRemnant->getA()>1) {
1280 if(prA>=1) { <<
1281 // Set the mass 872 // Set the mass
1282 const G4double aMass = theProjectileRem 873 const G4double aMass = theProjectileRemnant->getInvariantMass();
1283 theProjectileRemnant->setMass(aMass); 874 theProjectileRemnant->setMass(aMass);
1284 875
1285 // Compute the excitation energy from t 876 // Compute the excitation energy from the invariant mass
1286 const G4double anExcitationEnergy = aMa 877 const G4double anExcitationEnergy = aMass
1287 - ParticleTable::getTableMass(prA, th << 878 - ParticleTable::getTableMass(theProjectileRemnant->getA(), theProjectileRemnant->getZ());
1288 879
1289 // Set the excitation energy 880 // Set the excitation energy
1290 theProjectileRemnant->setExcitationEner 881 theProjectileRemnant->setExcitationEnergy(anExcitationEnergy);
1291 882
1292 // No spin! << 883 // Set the spin
1293 theProjectileRemnant->setSpin(ThreeVect << 884 theProjectileRemnant->setSpin(DeJongSpin::shoot(theProjectileRemnant->getNumberStoredComponents(), theProjectileRemnant->getA()));
1294 885
1295 // Set the emission time 886 // Set the emission time
1296 theProjectileRemnant->setEmissionTime(a 887 theProjectileRemnant->setEmissionTime(anEmissionTime);
1297 } << 888
>> 889 // Put it in the outgoing list
>> 890 theStore->addToOutgoing(theProjectileRemnant);
>> 891
>> 892 // NULL theProjectileRemnant
>> 893 theProjectileRemnant = NULL;
>> 894 } else if(theProjectileRemnant->getA()==1) {
>> 895 // Put the nucleon in the outgoing list
>> 896 Particle *theNucleon = theProjectileRemnant->getParticles().front();
>> 897 theStore->addToOutgoing(theNucleon);
>> 898 // Delete the remnant
>> 899 deleteProjectileRemnant();
>> 900 } else
>> 901 deleteProjectileRemnant();
1298 } 902 }
1299 903
1300 } 904 }
1301 905
1302 #endif 906 #endif
1303 907