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