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