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 #include "G4INCLDecayAvatar.hh" 38 #include "G4INCLDecayAvatar.hh" 39 39 40 #include "G4INCLDeltaDecayChannel.hh" 40 #include "G4INCLDeltaDecayChannel.hh" 41 #include "G4INCLPionResonanceDecayChannel.hh" 41 #include "G4INCLPionResonanceDecayChannel.hh" 42 #include "G4INCLSigmaZeroDecayChannel.hh" 42 #include "G4INCLSigmaZeroDecayChannel.hh" 43 #include "G4INCLNeutralKaonDecayChannel.hh" 43 #include "G4INCLNeutralKaonDecayChannel.hh" 44 #include "G4INCLStrangeAbsorbtionChannel.hh" 44 #include "G4INCLStrangeAbsorbtionChannel.hh" 45 #include "G4INCLPauliBlocking.hh" 45 #include "G4INCLPauliBlocking.hh" 46 #include <sstream> 46 #include <sstream> 47 #include <string> 47 #include <string> 48 // #include <cassert> 48 // #include <cassert> 49 49 50 namespace G4INCL { 50 namespace G4INCL { 51 51 52 DecayAvatar::DecayAvatar(G4INCL::Particle *a 52 DecayAvatar::DecayAvatar(G4INCL::Particle *aParticle, G4double time, G4INCL::Nucleus *n, G4bool force) 53 : InteractionAvatar(time, n, aParticle), f 53 : InteractionAvatar(time, n, aParticle), forced(force), 54 incidentDirection(aParticle->getMomentum 54 incidentDirection(aParticle->getMomentum()) 55 { 55 { 56 setType(DecayAvatarType); 56 setType(DecayAvatarType); 57 } 57 } 58 58 59 DecayAvatar::DecayAvatar(G4INCL::Particle *a 59 DecayAvatar::DecayAvatar(G4INCL::Particle *aParticle, G4INCL::Particle *bParticle, G4double time, G4INCL::Nucleus *n, G4bool force) 60 : InteractionAvatar(time, n, aParticle, bP 60 : InteractionAvatar(time, n, aParticle, bParticle), forced(force), 61 incidentDirection(aParticle->getMomentum 61 incidentDirection(aParticle->getMomentum()) 62 { 62 { 63 setType(DecayAvatarType); 63 setType(DecayAvatarType); 64 } 64 } 65 65 66 DecayAvatar::~DecayAvatar() { 66 DecayAvatar::~DecayAvatar() { 67 67 68 } 68 } 69 69 70 G4INCL::IChannel* DecayAvatar::getChannel() 70 G4INCL::IChannel* DecayAvatar::getChannel() { 71 if(!particle2){ 71 if(!particle2){ 72 if(particle1->isDelta()) { 72 if(particle1->isDelta()) { 73 INCL_DEBUG("DeltaDecayChannel chosen. 73 INCL_DEBUG("DeltaDecayChannel chosen." << '\n'); 74 return new DeltaDecayChannel(particle 74 return new DeltaDecayChannel(particle1, incidentDirection); 75 } 75 } 76 else if(particle1->isEta() || particle1 76 else if(particle1->isEta() || particle1->isOmega()) { 77 INCL_DEBUG("PionResonanceDecayChannel 77 INCL_DEBUG("PionResonanceDecayChannel chosen." << '\n'); 78 return new PionResonanceDecayChannel( 78 return new PionResonanceDecayChannel(particle1, incidentDirection); 79 } 79 } 80 else if(particle1->getType() == SigmaZe 80 else if(particle1->getType() == SigmaZero) { 81 INCL_DEBUG("SigmaZeroDecayChannel cho 81 INCL_DEBUG("SigmaZeroDecayChannel chosen." << '\n'); 82 return new SigmaZeroDecayChannel(part 82 return new SigmaZeroDecayChannel(particle1, incidentDirection); 83 } 83 } 84 else if(particle1->getType() == KZero | 84 else if(particle1->getType() == KZero || particle1->getType() == KZeroBar) { 85 INCL_DEBUG("NeutralKaonDecayChannel c 85 INCL_DEBUG("NeutralKaonDecayChannel chosen." << '\n'); 86 return new NeutralKaonDecayChannel(pa 86 return new NeutralKaonDecayChannel(particle1); 87 } 87 } 88 } 88 } 89 else if(((particle1->isAntiKaon() || parti 89 else if(((particle1->isAntiKaon() || particle1->isSigma()) && particle2->isNucleon()) || ((particle2->isAntiKaon() || particle2->isSigma()) && particle1->isNucleon())){ 90 INCL_DEBUG("StrangeAbsorbtion." << '\n') 90 INCL_DEBUG("StrangeAbsorbtion." << '\n'); 91 return new StrangeAbsorbtionChannel(part 91 return new StrangeAbsorbtionChannel(particle1, particle2); 92 } 92 } 93 return NULL; 93 return NULL; 94 } 94 } 95 95 96 void DecayAvatar::preInteraction() { 96 void DecayAvatar::preInteraction() { 97 InteractionAvatar::preInteraction(); 97 InteractionAvatar::preInteraction(); 98 } 98 } 99 99 100 void DecayAvatar::postInteraction(FinalState 100 void DecayAvatar::postInteraction(FinalState *fs) { 101 // Make sure we have at least two particle 101 // Make sure we have at least two particles in the final state 102 // Removed because of neutral kaon decay 102 // Removed because of neutral kaon decay 103 103 104 // assert((fs->getModifiedParticles().size() + 104 // assert((fs->getModifiedParticles().size() + fs->getCreatedParticles().size() - fs->getDestroyedParticles().size() >= 2) || ((*fs->getModifiedParticles().begin())->getType() == KShort || (*fs->getModifiedParticles().begin())->getType() == KLong )); 105 //assert((fs->getModifiedParticles().size( 105 //assert((fs->getModifiedParticles().size() + fs->getCreatedParticles().size() - fs->getDestroyedParticles().size() >= 1)); 106 if(!forced) { // Normal decay 106 if(!forced) { // Normal decay 107 // Call the postInteraction method of th 107 // Call the postInteraction method of the parent class 108 // (provides Pauli blocking and enforces 108 // (provides Pauli blocking and enforces energy conservation) 109 InteractionAvatar::postInteraction(fs); 109 InteractionAvatar::postInteraction(fs); 110 110 111 if(fs->getValidity() == PauliBlockedFS) 111 if(fs->getValidity() == PauliBlockedFS) 112 /* If the decay was Pauli-blocked, mak 112 /* If the decay was Pauli-blocked, make sure the propagation model 113 * generates a new decay avatar on the 113 * generates a new decay avatar on the next call to propagate(). 114 * 114 * 115 * \bug{Note that we don't generate ne 115 * \bug{Note that we don't generate new decay avatars for deltas that 116 * could not satisfy energy conservati 116 * could not satisfy energy conservation. This is in keeping with 117 * INCL4.6, but doesn't seem to make m 117 * INCL4.6, but doesn't seem to make much sense to me (DM), as energy 118 * conservation can be impossible to s 118 * conservation can be impossible to satisfy due to weird local-energy 119 * conditions, for example, that evolv 119 * conditions, for example, that evolve with time.} 120 */ 120 */ 121 fs->addModifiedParticle(particle1); 121 fs->addModifiedParticle(particle1); 122 } else { // Forced decay 122 } else { // Forced decay 123 modified = fs->getModifiedParticles(); 123 modified = fs->getModifiedParticles(); 124 created = fs->getCreatedParticles(); 124 created = fs->getCreatedParticles(); 125 Destroyed = fs->getDestroyedParticles(); 125 Destroyed = fs->getDestroyedParticles(); 126 modifiedAndCreated = modified; 126 modifiedAndCreated = modified; 127 modifiedAndCreated.insert(modifiedAndCre 127 modifiedAndCreated.insert(modifiedAndCreated.end(), created.begin(), created.end()); 128 ModifiedAndDestroyed = modified; 128 ModifiedAndDestroyed = modified; 129 ModifiedAndDestroyed.insert(ModifiedAndD 129 ModifiedAndDestroyed.insert(ModifiedAndDestroyed.end(), Destroyed.begin(), Destroyed.end()); 130 130 131 std::vector<G4int> newBiasCollisionVecto 131 std::vector<G4int> newBiasCollisionVector; 132 newBiasCollisionVector = ModifiedAndDest 132 newBiasCollisionVector = ModifiedAndDestroyed.getParticleListBiasVector(); 133 for(ParticleIter i=modifiedAndCreated.be 133 for(ParticleIter i=modifiedAndCreated.begin(), e=modifiedAndCreated.end(); i!=e; ++i ) { 134 (*i)->setBiasCollisionVector(newBiasColl 134 (*i)->setBiasCollisionVector(newBiasCollisionVector); 135 } 135 } 136 // Try to enforce energy conservation 136 // Try to enforce energy conservation 137 fs->setTotalEnergyBeforeInteraction(oldT 137 fs->setTotalEnergyBeforeInteraction(oldTotalEnergy); 138 const G4bool success = enforceEnergyCons 138 const G4bool success = enforceEnergyConservation(fs); 139 if(!success) { 139 if(!success) { 140 INCL_DEBUG("Enforcing energy conservat 140 INCL_DEBUG("Enforcing energy conservation: failed!" << '\n'); 141 141 142 if(theNucleus) { 142 if(theNucleus) { 143 // Restore the state of the initial 143 // Restore the state of the initial particles 144 restoreParticles(); 144 restoreParticles(); 145 145 146 // Delete newly created particles 146 // Delete newly created particles 147 for(ParticleIter i=created.begin(), 147 for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i ) 148 delete *i; 148 delete *i; 149 149 150 fs->reset(); 150 fs->reset(); 151 fs->makeNoEnergyConservation(); 151 fs->makeNoEnergyConservation(); 152 fs->setTotalEnergyBeforeInteraction( 152 fs->setTotalEnergyBeforeInteraction(0.0); 153 153 154 return; // Interaction is blocked. R 154 return; // Interaction is blocked. Return an empty final state. 155 } else { 155 } else { 156 // If there is no nucleus we have to 156 // If there is no nucleus we have to continue anyway, even if energy 157 // conservation failed. We cannot af 157 // conservation failed. We cannot afford producing unphysical 158 // remnants. 158 // remnants. 159 INCL_DEBUG("No nucleus, continuing a 159 INCL_DEBUG("No nucleus, continuing anyway." << '\n'); 160 } 160 } 161 } else { 161 } else { 162 INCL_DEBUG("Enforcing energy conservat 162 INCL_DEBUG("Enforcing energy conservation: success!" << '\n'); 163 } 163 } 164 164 165 if(theNucleus) { 165 if(theNucleus) { 166 // Test CDPP blocking 166 // Test CDPP blocking 167 G4bool isCDPPBlocked = Pauli::isCDPPBl 167 G4bool isCDPPBlocked = Pauli::isCDPPBlocked(created, theNucleus); 168 168 169 if(isCDPPBlocked) { 169 if(isCDPPBlocked) { 170 INCL_DEBUG("CDPP: Blocked!" << '\n') 170 INCL_DEBUG("CDPP: Blocked!" << '\n'); 171 171 172 // Restore the state of both particl 172 // Restore the state of both particles 173 restoreParticles(); 173 restoreParticles(); 174 174 175 // Delete newly created particles 175 // Delete newly created particles 176 for(ParticleIter i=created.begin(), 176 for(ParticleIter i=created.begin(), e=created.end(); i!=e; ++i ) 177 delete *i; 177 delete *i; 178 178 179 fs->reset(); 179 fs->reset(); 180 fs->makePauliBlocked(); 180 fs->makePauliBlocked(); 181 fs->setTotalEnergyBeforeInteraction( 181 fs->setTotalEnergyBeforeInteraction(0.0); 182 182 183 return; // Interaction is blocked. R 183 return; // Interaction is blocked. Return an empty final state. 184 } 184 } 185 INCL_DEBUG("CDPP: Allowed!" << '\n'); 185 INCL_DEBUG("CDPP: Allowed!" << '\n'); 186 186 187 } 187 } 188 } 188 } 189 // If there is a nucleus, increment the co 189 // If there is a nucleus, increment the counters 190 if(theNucleus) { 190 if(theNucleus) { 191 switch(fs->getValidity()) { 191 switch(fs->getValidity()) { 192 case PauliBlockedFS: 192 case PauliBlockedFS: 193 theNucleus->getStore()->getBook().in 193 theNucleus->getStore()->getBook().incrementBlockedDecays(); 194 break; 194 break; 195 case NoEnergyConservationFS: 195 case NoEnergyConservationFS: 196 case ParticleBelowFermiFS: 196 case ParticleBelowFermiFS: 197 case ParticleBelowZeroFS: 197 case ParticleBelowZeroFS: 198 break; 198 break; 199 case ValidFS: 199 case ValidFS: 200 theNucleus->getStore()->getBook().in 200 theNucleus->getStore()->getBook().incrementAcceptedDecays(); 201 } 201 } 202 } 202 } 203 203 204 return; 204 return; 205 } 205 } 206 206 207 std::string DecayAvatar::dump() const { 207 std::string DecayAvatar::dump() const { 208 std::stringstream ss; 208 std::stringstream ss; 209 ss << "(avatar " << theTime << " 'decay" < 209 ss << "(avatar " << theTime << " 'decay" << '\n' 210 << "(list " << '\n' 210 << "(list " << '\n' 211 << particle1->dump() 211 << particle1->dump() 212 << "))" << '\n'; 212 << "))" << '\n'; 213 return ss.str(); 213 return ss.str(); 214 } 214 } 215 } 215 } 216 216