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