Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLInteractionAvatar.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/inclxx/incl_physics/src/G4INCLInteractionAvatar.cc (Version 11.3.0) and /processes/hadronic/models/inclxx/incl_physics/src/G4INCLInteractionAvatar.cc (Version 9.5)


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