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