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