Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // INCL++ intra-nuclear cascade model 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France << 27 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics 28 // Joseph Cugnon, University of Liege, Belgium << 28 // Davide Mancusi, CEA 29 // Jean-Christophe David, CEA-Saclay, France << 29 // Alain Boudard, CEA 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H << 30 // Sylvie Leray, CEA 31 // Sylvie Leray, CEA-Saclay, France << 31 // Joseph Cugnon, University of Liege 32 // Davide Mancusi, CEA-Saclay, France << 32 // >> 33 // INCL++ revision: v5.0_rc3 33 // 34 // 34 #define INCLXX_IN_GEANT4_MODE 1 35 #define INCLXX_IN_GEANT4_MODE 1 35 36 36 #include "globals.hh" 37 #include "globals.hh" 37 38 38 /* \file G4INCLInteractionAvatar.cc 39 /* \file G4INCLInteractionAvatar.cc 39 * \brief Virtual class for interaction avatar << 40 * \brief Virtual class for G4interaction avatars. 40 * 41 * 41 * This class is inherited by decay and collis 42 * This class is inherited by decay and collision avatars. The goal is to 42 * provide a uniform treatment of common physi 43 * provide a uniform treatment of common physics, such as Pauli blocking, 43 * enforcement of energy conservation, etc. 44 * enforcement of energy conservation, etc. 44 * 45 * 45 * \date Mar 1st, 2011 << 46 * Created on: Mar 1st, 2011 46 * \author Davide Mancusi << 47 * Author: Davide Mancusi 47 */ 48 */ 48 49 49 #include "G4INCLInteractionAvatar.hh" 50 #include "G4INCLInteractionAvatar.hh" 50 #include "G4INCLKinematicsUtils.hh" 51 #include "G4INCLKinematicsUtils.hh" 51 #include "G4INCLCrossSections.hh" 52 #include "G4INCLCrossSections.hh" 52 #include "G4INCLPauliBlocking.hh" 53 #include "G4INCLPauliBlocking.hh" 53 #include "G4INCLRootFinder.hh" 54 #include "G4INCLRootFinder.hh" 54 #include "G4INCLLogger.hh" 55 #include "G4INCLLogger.hh" 55 #include "G4INCLConfigEnums.hh" 56 #include "G4INCLConfigEnums.hh" 56 #include "G4INCLConfig.hh" << 57 // #include <cassert> 57 // #include <cassert> 58 58 59 namespace G4INCL { 59 namespace G4INCL { 60 60 61 const G4double InteractionAvatar::locEAccura 61 const G4double InteractionAvatar::locEAccuracy = 1.E-4; 62 const G4int InteractionAvatar::maxIterLocE = 62 const G4int InteractionAvatar::maxIterLocE = 50; 63 G4ThreadLocal Particle *InteractionAvatar::b << 64 G4ThreadLocal Particle *InteractionAvatar::b << 65 63 66 InteractionAvatar::InteractionAvatar(G4doubl 64 InteractionAvatar::InteractionAvatar(G4double time, G4INCL::Nucleus *n, G4INCL::Particle *p1) 67 : IAvatar(time), theNucleus(n), 65 : IAvatar(time), theNucleus(n), 68 particle1(p1), particle2(NULL), << 66 particle1(p1), particle2(NULL), isPiN(false) 69 isPiN(false), << 70 weight(1.), << 71 violationEFunctor(NULL) << 72 { 67 { 73 } 68 } 74 69 75 InteractionAvatar::InteractionAvatar(G4doubl 70 InteractionAvatar::InteractionAvatar(G4double time, G4INCL::Nucleus *n, G4INCL::Particle *p1, 76 G4INCL::Particle *p2) 71 G4INCL::Particle *p2) 77 : IAvatar(time), theNucleus(n), 72 : IAvatar(time), theNucleus(n), 78 particle1(p1), particle2(p2), 73 particle1(p1), particle2(p2), 79 isPiN((p1->isPion() && p2->isNucleon()) || << 74 isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())) 80 weight(1.), << 81 violationEFunctor(NULL) << 82 { 75 { 83 } 76 } 84 77 85 InteractionAvatar::~InteractionAvatar() { 78 InteractionAvatar::~InteractionAvatar() { 86 } 79 } 87 80 88 void InteractionAvatar::deleteBackupParticle << 89 delete backupParticle1; << 90 if(backupParticle2) << 91 delete backupParticle2; << 92 backupParticle1 = NULL; << 93 backupParticle2 = NULL; << 94 } << 95 << 96 void InteractionAvatar::preInteractionBlocki 81 void InteractionAvatar::preInteractionBlocking() { 97 if(backupParticle1) << 82 oldParticle1Type = particle1->getType(); 98 (*backupParticle1) = (*particle1); << 83 oldParticle1Energy = particle1->getEnergy(); 99 else << 84 oldParticle1Potential = particle1->getPotentialEnergy(); 100 backupParticle1 = new Particle(*particle << 85 oldParticle1Momentum = particle1->getMomentum(); >> 86 oldParticle1Position = particle1->getPosition(); >> 87 oldParticle1Mass = particle1->getMass(); >> 88 oldParticle1Helicity = particle1->getHelicity(); 101 89 102 if(particle2) { 90 if(particle2) { 103 if(backupParticle2) << 91 oldParticle2Type = particle2->getType(); 104 (*backupParticle2) = (*particle2); << 92 oldParticle2Energy = particle2->getEnergy(); 105 else << 93 oldParticle2Potential = particle2->getPotentialEnergy(); 106 backupParticle2 = new Particle(*partic << 94 oldParticle2Momentum = particle2->getMomentum(); 107 << 95 oldParticle2Position = particle2->getPosition(); 108 oldTotalEnergy = particle1->getEnergy() << 96 oldParticle2Mass = particle2->getMass(); >> 97 oldParticle2Helicity = particle2->getHelicity(); >> 98 oldTotalEnergy = oldParticle1Energy + oldParticle2Energy 109 - particle1->getPotentialEnergy() - pa 99 - particle1->getPotentialEnergy() - particle2->getPotentialEnergy(); 110 oldXSec = CrossSections::total(particle1 100 oldXSec = CrossSections::total(particle1, particle2); 111 } else { 101 } else { 112 oldTotalEnergy = particle1->getEnergy() << 102 oldTotalEnergy = oldParticle1Energy - particle1->getPotentialEnergy(); 113 } 103 } 114 } 104 } 115 105 116 void InteractionAvatar::preInteractionLocalE 106 void InteractionAvatar::preInteractionLocalEnergy(Particle * const p) { 117 if(!theNucleus || p->isMeson() || p->isPho << 107 if(!theNucleus || p->isPion()) return; // Local energy does not make any sense without a nucleus 118 108 119 if(shouldUseLocalEnergy()) 109 if(shouldUseLocalEnergy()) 120 KinematicsUtils::transformToLocalEnergyF 110 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p); 121 } 111 } 122 112 123 void InteractionAvatar::preInteraction() { 113 void InteractionAvatar::preInteraction() { 124 preInteractionBlocking(); 114 preInteractionBlocking(); 125 115 126 preInteractionLocalEnergy(particle1); 116 preInteractionLocalEnergy(particle1); 127 117 128 if(particle2) { 118 if(particle2) { 129 preInteractionLocalEnergy(particle2); 119 preInteractionLocalEnergy(particle2); 130 boostVector = KinematicsUtils::makeBoost << 120 if(!isPiN) { 131 particle2->boost(boostVector); << 121 boostVector = KinematicsUtils::makeBoostVector(particle1, particle2); >> 122 particle2->boost(boostVector); >> 123 } 132 } else { 124 } else { 133 boostVector = particle1->getMomentum()/p 125 boostVector = particle1->getMomentum()/particle1->getEnergy(); 134 } 126 } 135 particle1->boost(boostVector); << 127 if(!isPiN) >> 128 particle1->boost(boostVector); 136 } 129 } 137 130 138 G4bool InteractionAvatar::bringParticleInsid 131 G4bool InteractionAvatar::bringParticleInside(Particle * const p) { 139 if(!theNucleus) << 140 return false; << 141 << 142 ThreeVector pos = p->getPosition(); 132 ThreeVector pos = p->getPosition(); 143 p->rpCorrelate(); << 144 G4double pos2 = pos.mag2(); 133 G4double pos2 = pos.mag2(); 145 const G4double r = theNucleus->getSurfaceR 134 const G4double r = theNucleus->getSurfaceRadius(p); 146 short iterations=0; 135 short iterations=0; 147 const short maxIterations=50; 136 const short maxIterations=50; 148 137 149 if(pos2 < r*r) return true; 138 if(pos2 < r*r) return true; 150 139 151 while( pos2 >= r*r && iterations<maxIterat << 140 while( pos2 >= r*r && iterations<maxIterations ) 152 { 141 { 153 pos *= std::sqrt(r*r*0.9801/pos2); // 0. << 142 pos *= std::sqrt(r*r*0.99/pos2); 154 pos2 = pos.mag2(); 143 pos2 = pos.mag2(); 155 iterations++; 144 iterations++; 156 } 145 } 157 if( iterations < maxIterations) 146 if( iterations < maxIterations) 158 { 147 { 159 INCL_DEBUG("Particle position vector len << 148 DEBUG("Particle position vector length was : " << p->getPosition().mag() << ", rescaled to: " << pos.mag() << std::endl); 160 p->setPosition(pos); 149 p->setPosition(pos); 161 return true; 150 return true; 162 } 151 } 163 else 152 else 164 return false; 153 return false; 165 } 154 } 166 155 167 void InteractionAvatar::postInteraction(Fina << 156 FinalState *InteractionAvatar::postInteraction(FinalState *fs) { 168 INCL_DEBUG("postInteraction: final state: << 157 ParticleList modified = fs->getModifiedParticles(); 169 modified = fs->getModifiedParticles(); << 158 ParticleList modifiedAndCreated = modified; 170 created = fs->getCreatedParticles(); << 159 ParticleList created = fs->getCreatedParticles(); 171 Destroyed = fs->getDestroyedParticles(); << 172 modifiedAndCreated = modified; << 173 modifiedAndCreated.insert(modifiedAndCreat 160 modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end()); 174 ModifiedAndDestroyed = modified; << 175 ModifiedAndDestroyed.insert(ModifiedAndDes << 176 161 177 // Boost back to lab << 162 if(!isPiN) { 178 modifiedAndCreated.boost(-boostVector); << 163 // Boost back to lab >> 164 for( ParticleIter i = modifiedAndCreated.begin(); i != modifiedAndCreated.end(); ++i ) >> 165 (*i)->boost(-boostVector); >> 166 } 179 167 180 // If there is no Nucleus, just return 168 // If there is no Nucleus, just return 181 if(!theNucleus) return; << 169 if(!theNucleus) return fs; 182 170 183 // Mark pions and kaons that have been cre << 171 // Mark pions that have been created outside their well (we will force them 184 // to be emitted later). 172 // to be emitted later). 185 for(ParticleIter i=created.begin(), e=crea << 173 for( ParticleIter i = created.begin(); i != created.end(); ++i ) 186 if(((*i)->isPion() || (*i)->isKaon() || << 174 if((*i)->isPion() && (*i)->getPosition().mag() > theNucleus->getSurfaceRadius(*i)) { 187 (*i)->makeParticipant(); 175 (*i)->makeParticipant(); 188 (*i)->setOutOfWell(); 176 (*i)->setOutOfWell(); 189 fs->addOutgoingParticle(*i); 177 fs->addOutgoingParticle(*i); 190 INCL_DEBUG("Pion was created outside i << 178 DEBUG("Pion was created outside its potential well." << std::endl 191 << (*i)->print()); << 179 << (*i)->prG4int()); 192 } 180 } 193 181 194 // Try to enforce energy conservation 182 // Try to enforce energy conservation 195 fs->setTotalEnergyBeforeInteraction(oldTot 183 fs->setTotalEnergyBeforeInteraction(oldTotalEnergy); 196 G4bool success = enforceEnergyConservation << 184 G4bool success = true; >> 185 if(!isPiN || shouldUseLocalEnergy()) >> 186 success = enforceEnergyConservation(fs); 197 if(!success) { 187 if(!success) { 198 INCL_DEBUG("Enforcing energy conservatio << 188 DEBUG("Enforcing energy conservation: failed!" << std::endl); 199 189 200 // Restore the state of the initial part 190 // Restore the state of the initial particles 201 restoreParticles(); 191 restoreParticles(); 202 192 203 // Delete newly created particles 193 // Delete newly created particles 204 for(ParticleIter i=created.begin(), e=cr << 194 for( ParticleIter i = created.begin(); i != created.end(); ++i ) 205 delete *i; 195 delete *i; 206 196 207 fs->reset(); << 197 FinalState *fsBlocked = new FinalState; 208 fs->makeNoEnergyConservation(); << 198 delete fs; 209 fs->setTotalEnergyBeforeInteraction(0.0) << 199 fsBlocked->makeNoEnergyConservation(); >> 200 fsBlocked->setTotalEnergyBeforeInteraction(0.0); 210 201 211 return; // Interaction is blocked. Retur << 202 return fsBlocked; // Interaction is blocked. Return an empty final state. 212 } 203 } 213 INCL_DEBUG("Enforcing energy conservation: << 204 DEBUG("Enforcing energy conservation: success!" << std::endl); 214 << 215 INCL_DEBUG("postInteraction after energy c << 216 205 217 // Check that outgoing delta resonances ca 206 // Check that outgoing delta resonances can decay to pi-N 218 for(ParticleIter i=modified.begin(), e=mod << 207 for( ParticleIter i = modified.begin(); i != modified.end(); ++i ) 219 if((*i)->isDelta() && 208 if((*i)->isDelta() && 220 (*i)->getMass() < ParticleTable::min << 209 (*i)->getMass() < ParticleTable::effectiveDeltaDecayThreshold) { 221 INCL_DEBUG("Mass of the produced delta << 210 DEBUG("Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" << 222 (*i)->getMass() << '\n'); << 211 (*i)->getMass() << std::endl); 223 212 224 // Restore the state of the initial pa 213 // Restore the state of the initial particles 225 restoreParticles(); 214 restoreParticles(); 226 215 227 // Delete newly created particles 216 // Delete newly created particles 228 for(ParticleIter j=created.begin(), en << 217 for( ParticleIter i = created.begin(); i != created.end(); ++i ) 229 delete *j; << 218 delete *i; 230 219 231 fs->reset(); << 220 FinalState *fsBlocked = new FinalState; 232 fs->makeNoEnergyConservation(); << 221 delete fs; 233 fs->setTotalEnergyBeforeInteraction(0. << 222 fsBlocked->makeNoEnergyConservation(); >> 223 fsBlocked->setTotalEnergyBeforeInteraction(0.0); 234 224 235 return; // Interaction is blocked. Ret << 225 return fsBlocked; // Interaction is blocked. Return an empty final state. 236 } 226 } 237 227 238 INCL_DEBUG("Random seeds before Pauli bloc << 239 // Test Pauli blocking 228 // Test Pauli blocking 240 G4bool isBlocked = Pauli::isBlocked(modifi 229 G4bool isBlocked = Pauli::isBlocked(modifiedAndCreated, theNucleus); 241 230 242 if(isBlocked) { 231 if(isBlocked) { 243 INCL_DEBUG("Pauli: Blocked!" << '\n'); << 232 DEBUG("Pauli: Blocked!" << std::endl); 244 233 245 // Restore the state of the initial part 234 // Restore the state of the initial particles 246 restoreParticles(); 235 restoreParticles(); 247 236 248 // Delete newly created particles 237 // Delete newly created particles 249 for(ParticleIter i=created.begin(), e=cr << 238 for( ParticleIter i = created.begin(); i != created.end(); ++i ) 250 delete *i; 239 delete *i; 251 240 252 fs->reset(); << 241 FinalState *fsBlocked = new FinalState; 253 fs->makePauliBlocked(); << 242 delete fs; 254 fs->setTotalEnergyBeforeInteraction(0.0) << 243 fsBlocked->makePauliBlocked(); >> 244 fsBlocked->setTotalEnergyBeforeInteraction(0.0); 255 245 256 return; // Interaction is blocked. Retur << 246 return fsBlocked; // Interaction is blocked. Return an empty final state. 257 } 247 } 258 INCL_DEBUG("Pauli: Allowed!" << '\n'); << 248 DEBUG("Pauli: Allowed!" << std::endl); 259 249 260 // Test CDPP blocking 250 // Test CDPP blocking 261 G4bool isCDPPBlocked = Pauli::isCDPPBlocke 251 G4bool isCDPPBlocked = Pauli::isCDPPBlocked(created, theNucleus); 262 252 263 if(isCDPPBlocked) { 253 if(isCDPPBlocked) { 264 INCL_DEBUG("CDPP: Blocked!" << '\n'); << 254 DEBUG("CDPP: Blocked!" << std::endl); 265 255 266 // Restore the state of the initial part 256 // Restore the state of the initial particles 267 restoreParticles(); 257 restoreParticles(); 268 258 269 // Delete newly created particles 259 // Delete newly created particles 270 for(ParticleIter i=created.begin(), e=cr << 260 for( ParticleIter i = created.begin(); i != created.end(); ++i ) 271 delete *i; 261 delete *i; 272 262 273 fs->reset(); << 263 FinalState *fsBlocked = new FinalState; 274 fs->makePauliBlocked(); << 264 delete fs; 275 fs->setTotalEnergyBeforeInteraction(0.0) << 265 fsBlocked->makePauliBlocked(); >> 266 fsBlocked->setTotalEnergyBeforeInteraction(0.0); 276 267 277 return; // Interaction is blocked. Retur << 268 return fsBlocked; // Interaction is blocked. Return an empty final state. 278 } 269 } 279 INCL_DEBUG("CDPP: Allowed!" << '\n'); << 270 DEBUG("CDPP: Allowed!" << std::endl); 280 271 281 // If all went well, try to bring particle 272 // If all went well, try to bring particles inside the nucleus... 282 for(ParticleIter i=modifiedAndCreated.begi << 273 for( ParticleIter i = modifiedAndCreated.begin(); i != modifiedAndCreated.end(); ++i ) 283 { 274 { 284 // ...except for pions beyond their surf 275 // ...except for pions beyond their surface radius. 285 if((*i)->isOutOfWell()) continue; 276 if((*i)->isOutOfWell()) continue; 286 277 287 const G4bool successBringParticlesInside << 278 const G4bool success = bringParticleInside(*i); 288 if( !successBringParticlesInside ) { << 279 if( !success ) { 289 INCL_ERROR("Failed to bring particle i << 280 ERROR("Failed to bring particle inside the nucleus!" << std::endl); 290 } 281 } 291 } 282 } 292 283 293 // Collision accepted! 284 // Collision accepted! 294 // Biasing of the final state << 285 for( ParticleIter i = modifiedAndCreated.begin(); i != modifiedAndCreated.end(); ++i ) { 295 std::vector<G4int> newBiasCollisionVector; << 296 newBiasCollisionVector = ModifiedAndDestro << 297 if(std::fabs(weight-1.) > 1E-6){ << 298 newBiasCollisionVector.push_back(Particl << 299 Particle::FillINCLBiasVector(1./weight); << 300 weight = 1.; // useless? << 301 } << 302 for(ParticleIter i=modifiedAndCreated.begi << 303 (*i)->setBiasCollisionVector(newBiasColl << 304 if(!(*i)->isOutOfWell()) { 286 if(!(*i)->isOutOfWell()) { 305 // Decide if the particle should be ma << 287 // Decide if the particle should be made G4into a spectator 306 // (Back to spectator) << 307 G4bool goesBackToSpectator = false; 288 G4bool goesBackToSpectator = false; 308 if((*i)->isNucleon() && theNucleus->ge 289 if((*i)->isNucleon() && theNucleus->getStore()->getConfig()->getBackToSpectator()) { 309 G4double threshold = (*i)->getPotent << 290 const G4double threshold = (*i)->getPotentialEnergy() 310 if((*i)->getType()==Proton) << 291 - ParticleTable::getSeparationEnergy((*i)->getType()) 311 threshold += Math::twoThirds*theNu << 292 + theNucleus->getStore()->getConfig()->getBackToSpectatorThreshold(); 312 if((*i)->getKineticEnergy() < thresh 293 if((*i)->getKineticEnergy() < threshold) 313 goesBackToSpectator = true; 294 goesBackToSpectator = true; 314 } 295 } 315 // Thaw the particle propagation << 316 (*i)->thawPropagation(); << 317 << 318 // Increment or decrement the particip 296 // Increment or decrement the participant counters 319 if(goesBackToSpectator) { 297 if(goesBackToSpectator) { 320 INCL_DEBUG("The following particle g << 298 if((*i)->isParticipant()) { 321 << (*i)->print() << '\n'); << 299 theNucleus->getStore()->getBook()->decrementParticipants(); 322 if(!(*i)->isTargetSpectator()) { << 300 (*i)->makeSpectator(); 323 theNucleus->getStore()->getBook(). << 324 } 301 } 325 (*i)->makeTargetSpectator(); << 326 } else { 302 } else { 327 if((*i)->isTargetSpectator()) { << 303 if(!(*i)->isParticipant()) { 328 theNucleus->getStore()->getBook(). << 304 theNucleus->getStore()->getBook()->incrementParticipants(); >> 305 (*i)->makeParticipant(); 329 } 306 } 330 (*i)->makeParticipant(); << 331 } 307 } 332 } 308 } 333 } 309 } 334 ParticleList destroyed = fs->getDestroyedP 310 ParticleList destroyed = fs->getDestroyedParticles(); 335 for(ParticleIter i=destroyed.begin(), e=de << 311 for( ParticleIter i = destroyed.begin(); i != destroyed.end(); ++i ) 336 if(!(*i)->isTargetSpectator()) << 312 if((*i)->isParticipant()) 337 theNucleus->getStore()->getBook().decr << 313 theNucleus->getStore()->getBook()->decrementParticipants(); 338 return; << 314 >> 315 return fs; 339 } 316 } 340 317 341 void InteractionAvatar::restoreParticles() c 318 void InteractionAvatar::restoreParticles() const { 342 (*particle1) = (*backupParticle1); << 319 particle1->setType(oldParticle1Type); 343 if(particle2) << 320 particle1->setEnergy(oldParticle1Energy); 344 (*particle2) = (*backupParticle2); << 321 particle1->setPotentialEnergy(oldParticle1Potential); 345 } << 322 particle1->setMomentum(oldParticle1Momentum); >> 323 particle1->setPosition(oldParticle1Position); >> 324 particle1->setMass(oldParticle1Mass); >> 325 particle1->setHelicity(oldParticle1Helicity); 346 326 347 G4bool InteractionAvatar::shouldUseLocalEner << 327 if(particle2) { 348 if(!theNucleus) return false; << 328 particle2->setType(oldParticle2Type); 349 LocalEnergyType theLocalEnergyType; << 329 particle2->setEnergy(oldParticle2Energy); 350 if(theNucleus->getStore()->getConfig()->ge << 330 particle2->setPotentialEnergy(oldParticle2Potential); 351 theNucleus->getStore()->getConfig()->getP << 331 particle2->setMomentum(oldParticle2Momentum); 352 return false; << 332 particle2->setPosition(oldParticle2Position); >> 333 particle2->setMass(oldParticle2Mass); >> 334 particle2->setHelicity(oldParticle2Helicity); 353 } 335 } 354 if(getType()==DecayAvatarType || isPiN) << 355 theLocalEnergyType = theNucleus->getStor << 356 else << 357 theLocalEnergyType = theNucleus->getStor << 358 << 359 const G4bool firstAvatar = (theNucleus->ge << 360 return ((theLocalEnergyType == FirstCollis << 361 theLocalEnergyType == AlwaysLocalE << 362 } 336 } 363 337 364 G4bool InteractionAvatar::enforceEnergyConse 338 G4bool InteractionAvatar::enforceEnergyConservation(FinalState * const fs) { 365 // Set up the violationE calculation 339 // Set up the violationE calculation 366 const G4bool manyBodyFinalState = (modifie << 340 ParticleList modified = fs->getModifiedParticles(); 367 << 341 const G4bool manyBodyFinalState = (modified.size() + fs->getCreatedParticles().size() > 1); 368 if(manyBodyFinalState) 342 if(manyBodyFinalState) 369 violationEFunctor = new ViolationEMoment << 343 violationEFunctor = new ViolationEMomentumFunctor(theNucleus, fs, &boostVector, shouldUseLocalEnergy()); 370 else { 344 else { 371 if (modified.empty()) { << 345 Particle const * const p = modified.front(); 372 Particle * const p1 = created.front(); << 346 // The following condition is necessary for the functor to work 373 // The following condition is necessa << 347 // correctly. A similar condition exists in INCL4.6. 374 // correctly. A similar condition exi << 348 if(p->getMass() < ParticleTable::effectiveDeltaDecayThreshold) 375 if(p1->getMass() < ParticleTable::minD << 349 return false; 376 return false; << 350 violationEFunctor = new ViolationEEnergyFunctor(theNucleus, fs); 377 violationEFunctor = new ViolationEEner << 378 } << 379 else{ << 380 Particle * const p2 = modified.front() << 381 // The following condition is necessa << 382 // correctly. A similar condition exi << 383 if(p2->getMass() < ParticleTable::min << 384 return false; << 385 violationEFunctor = new ViolationEEne << 386 } << 387 } 351 } 388 352 389 // Apply the root-finding algorithm 353 // Apply the root-finding algorithm 390 const RootFinder::Solution theSolution = R << 354 const G4bool success = RootFinder::solve(violationEFunctor, 1.0); 391 if(theSolution.success) { // Apply the sol << 355 if(!success) { 392 (*violationEFunctor)(theSolution.x); << 356 WARN("Couldn't enforce energy conservation after an G4interaction, root-finding algorithm failed." << std::endl); 393 } else if(theNucleus){ << 394 INCL_DEBUG("Couldn't enforce energy cons << 395 theNucleus->getStore()->getBook().increm << 396 } 357 } 397 delete violationEFunctor; 358 delete violationEFunctor; 398 violationEFunctor = NULL; << 359 return success; 399 return theSolution.success; << 400 } 360 } 401 361 402 /* *** 362 /* *** *** 403 * *** InteractionAvatar::ViolationEMomentum 363 * *** InteractionAvatar::ViolationEMomentumFunctor methods *** 404 * *** 364 * *** ***/ 405 365 406 InteractionAvatar::ViolationEMomentumFunctor << 366 InteractionAvatar::ViolationEMomentumFunctor::ViolationEMomentumFunctor(Nucleus * const nucleus, FinalState const * const finalState, ThreeVector const * const boost, const G4bool localE) 407 RootFunctor(0., 1E6), << 367 : initialEnergy(finalState->getTotalEnergyBeforeInteraction()), 408 finalParticles(modAndCre), << 409 initialEnergy(totalEnergyBeforeInteraction << 410 theNucleus(nucleus), 368 theNucleus(nucleus), 411 boostVector(boost), 369 boostVector(boost), 412 shouldUseLocalEnergy(localE) 370 shouldUseLocalEnergy(localE) 413 { 371 { >> 372 // Set up the finalParticles list >> 373 finalParticles = finalState->getModifiedParticles(); >> 374 ParticleList created = finalState->getCreatedParticles(); >> 375 finalParticles.splice(finalParticles.end(), created); >> 376 414 // Store the particle momenta (necessary f 377 // Store the particle momenta (necessary for the calls to 415 // scaleParticleMomenta() to work) 378 // scaleParticleMomenta() to work) 416 for(ParticleIter i=finalParticles.begin(), << 379 particleMomenta.clear(); 417 (*i)->boost(boostVector); << 380 for(ParticleIter i=finalParticles.begin(); i!=finalParticles.end(); ++i) { >> 381 (*i)->boost(*boostVector); 418 particleMomenta.push_back((*i)->getMomen 382 particleMomenta.push_back((*i)->getMomentum()); 419 } 383 } 420 } 384 } 421 385 422 InteractionAvatar::ViolationEMomentumFunctor << 423 particleMomenta.clear(); << 424 } << 425 << 426 G4double InteractionAvatar::ViolationEMoment 386 G4double InteractionAvatar::ViolationEMomentumFunctor::operator()(const G4double alpha) const { 427 scaleParticleMomenta(alpha); 387 scaleParticleMomenta(alpha); 428 388 429 G4double deltaE = 0.0; 389 G4double deltaE = 0.0; 430 for(ParticleIter i=finalParticles.begin(), << 390 for(ParticleIter i=finalParticles.begin(); i!=finalParticles.end(); ++i) 431 deltaE += (*i)->getEnergy() - (*i)->getP 391 deltaE += (*i)->getEnergy() - (*i)->getPotentialEnergy(); 432 deltaE -= initialEnergy; 392 deltaE -= initialEnergy; 433 return deltaE; 393 return deltaE; 434 } 394 } 435 395 436 void InteractionAvatar::ViolationEMomentumFu 396 void InteractionAvatar::ViolationEMomentumFunctor::scaleParticleMomenta(const G4double alpha) const { 437 397 438 std::vector<ThreeVector>::const_iterator i << 398 std::list<ThreeVector>::const_iterator iP = particleMomenta.begin(); 439 for(ParticleIter i=finalParticles.begin(), << 399 for(ParticleIter i=finalParticles.begin(); i!=finalParticles.end(); ++i, ++iP) { 440 (*i)->setMomentum((*iP)*alpha); 400 (*i)->setMomentum((*iP)*alpha); 441 (*i)->adjustEnergyFromMomentum(); 401 (*i)->adjustEnergyFromMomentum(); 442 (*i)->rpCorrelate(); << 402 (*i)->boost(-(*boostVector)); 443 (*i)->boost(-boostVector); << 403 if(theNucleus) 444 if(theNucleus){ << 445 theNucleus->updatePotentialEnergy(*i); 404 theNucleus->updatePotentialEnergy(*i); 446 } else { << 405 else 447 (*i)->setPotentialEnergy(0.); 406 (*i)->setPotentialEnergy(0.); 448 } << 449 407 450 if(shouldUseLocalEnergy && !(*i)->isPion << 408 if(shouldUseLocalEnergy && !(*i)->isPion() && theNucleus) { // This translates AECSVT's loops 1, 3 and 4 451 !(*i)->isKaon() && !(*i)->isAntiKaon( << 409 // assert(theNucleus); // Local energy without a nucleus doesn't make sense 452 // assert(theNucleus); // Local energy without << 410 const G4double energy = (*i)->getEnergy(); // Store the energy of the particle 453 const G4double energy = (*i)->getEner << 411 G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // Initial value of local energy 454 G4double locE = KinematicsUtils::getL << 412 G4double locEOld; 455 G4double locEOld; << 413 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3; 456 G4double deltaLocE = InteractionAvata << 414 for(G4int iterLocE=0; 457 for(G4int iterLocE=0; << 458 deltaLocE>InteractionAvatar::locEA 415 deltaLocE>InteractionAvatar::locEAccuracy && iterLocE<InteractionAvatar::maxIterLocE; 459 ++iterLocE) { 416 ++iterLocE) { 460 locEOld = locE; 417 locEOld = locE; 461 (*i)->setEnergy(energy + locE); // U 418 (*i)->setEnergy(energy + locE); // Update the energy of the particle... 462 (*i)->adjustMomentumFromEnergy(); 419 (*i)->adjustMomentumFromEnergy(); 463 theNucleus->updatePotentialEnergy(*i 420 theNucleus->updatePotentialEnergy(*i); // ...update its potential energy... 464 locE = KinematicsUtils::getLocalEner 421 locE = KinematicsUtils::getLocalEnergy(theNucleus, *i); // ...and recompute locE. 465 deltaLocE = std::abs(locE-locEOld); 422 deltaLocE = std::abs(locE-locEOld); 466 } 423 } 467 } 424 } 468 << 469 //jlrs For lambdas and nuclei with masses hig << 470 if(shouldUseLocalEnergy && (*i)->isLam << 471 // assert(theNucleus); // Local energy without << 472 const G4double energy = (*i)->getEner << 473 G4double locE = KinematicsUtils::getL << 474 G4double locEOld; << 475 G4double deltaLocE = InteractionAvata << 476 for(G4int iterLocE=0; << 477 deltaLocE>InteractionAvatar::locEA << 478 ++iterLocE) { << 479 locEOld = locE; << 480 (*i)->setEnergy(energy + locE); // U << 481 (*i)->adjustMomentumFromEnergy(); << 482 theNucleus->updatePotentialEnergy(*i << 483 locE = KinematicsUtils::getLocalEner << 484 deltaLocE = std::abs(locE-locEOld); << 485 } << 486 } << 487 } 425 } 488 } 426 } 489 427 490 void InteractionAvatar::ViolationEMomentumFu 428 void InteractionAvatar::ViolationEMomentumFunctor::cleanUp(const G4bool success) const { 491 if(!success) 429 if(!success) 492 scaleParticleMomenta(1.); 430 scaleParticleMomenta(1.); 493 } 431 } 494 432 495 /* *** 433 /* *** *** 496 * *** InteractionAvatar::ViolationEEnergyFu 434 * *** InteractionAvatar::ViolationEEnergyFunctor methods *** 497 * *** 435 * *** ***/ 498 436 499 InteractionAvatar::ViolationEEnergyFunctor:: << 437 InteractionAvatar::ViolationEEnergyFunctor::ViolationEEnergyFunctor(Nucleus * const nucleus, FinalState const * const finalState) 500 RootFunctor(0., 1E6), << 438 : initialEnergy(finalState->getTotalEnergyBeforeInteraction()), 501 initialEnergy(totalEnergyBeforeInteraction << 502 theNucleus(nucleus), 439 theNucleus(nucleus), 503 theParticle(aParticle), << 440 theParticle(finalState->getModifiedParticles().front()), 504 theEnergy(theParticle->getEnergy()), 441 theEnergy(theParticle->getEnergy()), 505 theMomentum(theParticle->getMomentum()), 442 theMomentum(theParticle->getMomentum()), 506 energyThreshold(KinematicsUtils::energy(th << 443 energyThreshold(KinematicsUtils::energy(theMomentum,ParticleTable::effectiveDeltaDecayThreshold)) 507 shouldUseLocalEnergy(localE) << 508 { 444 { 509 // assert(theParticle->isDelta()); << 445 // assert(theNucleus); >> 446 // assert(finalState->getModifiedParticles().size()==1); >> 447 // assert(theParticle->isDelta()); 510 } 448 } 511 449 512 G4double InteractionAvatar::ViolationEEnergy 450 G4double InteractionAvatar::ViolationEEnergyFunctor::operator()(const G4double alpha) const { 513 setParticleEnergy(alpha); 451 setParticleEnergy(alpha); 514 return theParticle->getEnergy() - theParti 452 return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy; 515 } 453 } 516 454 517 void InteractionAvatar::ViolationEEnergyFunc 455 void InteractionAvatar::ViolationEEnergyFunctor::setParticleEnergy(const G4double alpha) const { 518 456 519 G4double locE; << 457 G4double locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle); // Initial value of local energy 520 if(shouldUseLocalEnergy) { << 521 // assert(theNucleus); // Local energy without << 522 locE = KinematicsUtils::getLocalEnergy(t << 523 } else << 524 locE = 0.; << 525 G4double locEOld; 458 G4double locEOld; 526 G4double deltaLocE = InteractionAvatar::lo 459 G4double deltaLocE = InteractionAvatar::locEAccuracy + 1E3; 527 for(G4int iterLocE=0; 460 for(G4int iterLocE=0; 528 deltaLocE>InteractionAvatar::locEAccur 461 deltaLocE>InteractionAvatar::locEAccuracy && iterLocE<InteractionAvatar::maxIterLocE; 529 ++iterLocE) { 462 ++iterLocE) { 530 locEOld = locE; 463 locEOld = locE; 531 G4double particleEnergy = energyThreshol << 464 const G4double particleEnergy = energyThreshold + alpha*(theEnergy-energyThreshold); 532 const G4double theMass2 = std::pow(parti << 465 const G4double theMass = std::sqrt(std::pow(particleEnergy,2.)-theMomentum.mag2()); 533 G4double theMass; << 534 if(theMass2>ParticleTable::minDeltaMass2 << 535 theMass = std::sqrt(theMass2); << 536 else { << 537 theMass = ParticleTable::minDeltaMass; << 538 particleEnergy = energyThreshold; << 539 } << 540 theParticle->setMass(theMass); 466 theParticle->setMass(theMass); 541 theParticle->setEnergy(particleEnergy); << 467 theParticle->setEnergy(particleEnergy + locE); // Update the energy of the particle... 542 if(theNucleus) { << 468 theParticle->adjustMomentumFromEnergy(); 543 theNucleus->updatePotentialEnergy(theP << 469 theNucleus->updatePotentialEnergy(theParticle); // ...update its potential energy... 544 if(shouldUseLocalEnergy) << 470 locE = KinematicsUtils::getLocalEnergy(theNucleus, theParticle); // ...and recompute locE. 545 locE = KinematicsUtils::getLocalEner << 546 else << 547 locE = 0.; << 548 } else << 549 locE = 0.; << 550 deltaLocE = std::abs(locE-locEOld); 471 deltaLocE = std::abs(locE-locEOld); 551 } 472 } 552 473 553 } 474 } 554 475 555 void InteractionAvatar::ViolationEEnergyFunc 476 void InteractionAvatar::ViolationEEnergyFunctor::cleanUp(const G4bool success) const { 556 if(!success) 477 if(!success) 557 setParticleEnergy(1.); 478 setParticleEnergy(1.); 558 } 479 } 559 480 560 } 481 } 561 482