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