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