Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/inclxx/incl_physics/src/G4INCLStandardPropagationModel.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/G4INCLStandardPropagationModel.cc (Version 11.3.0) and /processes/hadronic/models/inclxx/incl_physics/src/G4INCLStandardPropagationModel.cc (Version 11.2.2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // INCL++ intra-nuclear cascade model              26 // INCL++ intra-nuclear cascade model
 27 // Alain Boudard, CEA-Saclay, France               27 // Alain Boudard, CEA-Saclay, France
 28 // Joseph Cugnon, University of Liege, Belgium     28 // Joseph Cugnon, University of Liege, Belgium
 29 // Jean-Christophe David, CEA-Saclay, France       29 // Jean-Christophe David, CEA-Saclay, France
 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H     30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
 31 // Sylvie Leray, CEA-Saclay, France                31 // Sylvie Leray, CEA-Saclay, France
 32 // Davide Mancusi, CEA-Saclay, France              32 // Davide Mancusi, CEA-Saclay, France
 33 //                                                 33 //
 34 #define INCLXX_IN_GEANT4_MODE 1                    34 #define INCLXX_IN_GEANT4_MODE 1
 35                                                    35 
 36 #include "globals.hh"                              36 #include "globals.hh"
 37                                                    37 
 38 /*                                                 38 /*
 39  * StandardPropagationModel.cpp                    39  * StandardPropagationModel.cpp
 40  *                                                 40  *
 41  *  \date 4 juin 2009                              41  *  \date 4 juin 2009
 42  * \author Pekka Kaitaniemi                        42  * \author Pekka Kaitaniemi
 43  */                                                43  */
 44                                                    44 
 45 #include "G4INCLStandardPropagationModel.hh"       45 #include "G4INCLStandardPropagationModel.hh"
 46 #include "G4INCLPbarAtrestEntryChannel.hh"         46 #include "G4INCLPbarAtrestEntryChannel.hh"
 47 #include "G4INCLSurfaceAvatar.hh"                  47 #include "G4INCLSurfaceAvatar.hh"
 48 #include "G4INCLBinaryCollisionAvatar.hh"          48 #include "G4INCLBinaryCollisionAvatar.hh"
 49 #include "G4INCLDecayAvatar.hh"                    49 #include "G4INCLDecayAvatar.hh"
 50 #include "G4INCLCrossSections.hh"                  50 #include "G4INCLCrossSections.hh"
 51 #include "G4INCLRandom.hh"                         51 #include "G4INCLRandom.hh"
 52 #include <iostream>                                52 #include <iostream>
 53 #include <algorithm>                               53 #include <algorithm>
 54 #include "G4INCLLogger.hh"                         54 #include "G4INCLLogger.hh"
 55 #include "G4INCLGlobals.hh"                        55 #include "G4INCLGlobals.hh"
 56 #include "G4INCLKinematicsUtils.hh"                56 #include "G4INCLKinematicsUtils.hh"
 57 #include "G4INCLCoulombDistortion.hh"              57 #include "G4INCLCoulombDistortion.hh"
 58 #include "G4INCLDeltaDecayChannel.hh"              58 #include "G4INCLDeltaDecayChannel.hh"
 59 #include "G4INCLSigmaZeroDecayChannel.hh"          59 #include "G4INCLSigmaZeroDecayChannel.hh"
 60 #include "G4INCLPionResonanceDecayChannel.hh"      60 #include "G4INCLPionResonanceDecayChannel.hh"
 61 #include "G4INCLParticleEntryAvatar.hh"            61 #include "G4INCLParticleEntryAvatar.hh"
 62 #include "G4INCLIntersection.hh"                   62 #include "G4INCLIntersection.hh"
 63 #include <vector>                                  63 #include <vector>
 64                                                    64 
 65 namespace G4INCL {                                 65 namespace G4INCL {
 66                                                    66 
 67     StandardPropagationModel::StandardPropagat     67     StandardPropagationModel::StandardPropagationModel(LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType, const G4double hTime)
 68       :theNucleus(0), maximumTime(70.0), curre     68       :theNucleus(0), maximumTime(70.0), currentTime(0.0),
 69       hadronizationTime(hTime),                    69       hadronizationTime(hTime),
 70       firstAvatar(true),                           70       firstAvatar(true),
 71       theLocalEnergyType(localEnergyType),         71       theLocalEnergyType(localEnergyType),
 72       theLocalEnergyDeltaType(localEnergyDelta     72       theLocalEnergyDeltaType(localEnergyDeltaType)
 73     {                                              73     {
 74     }                                              74     }
 75                                                    75 
 76     StandardPropagationModel::~StandardPropaga     76     StandardPropagationModel::~StandardPropagationModel()
 77     {                                              77     {
 78       delete theNucleus;                           78       delete theNucleus;
 79     }                                              79     }
 80                                                    80 
 81     G4INCL::Nucleus* StandardPropagationModel:     81     G4INCL::Nucleus* StandardPropagationModel::getNucleus()
 82     {                                              82     {
 83       return theNucleus;                           83       return theNucleus;
 84     }                                              84     }
 85                                                    85 
 86 //D                                                86 //D
 87                                                    87 
 88     G4double StandardPropagationModel::shoot(P     88     G4double StandardPropagationModel::shoot(ParticleSpecies const &projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
 89       if(projectileSpecies.theType==Composite)     89       if(projectileSpecies.theType==Composite){
 90         return shootComposite(projectileSpecie     90         return shootComposite(projectileSpecies, kineticEnergy, impactParameter, phi);
 91       }                                            91       }
 92       else if(projectileSpecies.theType==antiP     92       else if(projectileSpecies.theType==antiProton && theNucleus->getAnnihilationType()!=Def){
 93         return shootAtrest(projectileSpecies.t     93         return shootAtrest(projectileSpecies.theType, kineticEnergy);
 94       }                                            94       }
 95       else{                                        95       else{
 96         return shootParticle(projectileSpecies     96         return shootParticle(projectileSpecies.theType, kineticEnergy, impactParameter, phi);
 97       }                                            97       }
 98     }                                              98     }
 99                                                    99 
100 //D                                               100 //D
101                                                   101 
102     G4double StandardPropagationModel::shootAt    102     G4double StandardPropagationModel::shootAtrest(ParticleType const t, const G4double kineticEnergy) {
103       theNucleus->setParticleNucleusCollision(    103       theNucleus->setParticleNucleusCollision(); 
104       currentTime = 0.0;                          104       currentTime = 0.0;
105                                                   105 
106       // Create final state particles             106       // Create final state particles
107       const G4double projectileMass = Particle    107       const G4double projectileMass = ParticleTable::getTableParticleMass(t);
108       G4double energy = kineticEnergy + projec    108       G4double energy = kineticEnergy + projectileMass;
109       G4double momentumZ = std::sqrt(energy*en    109       G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
110       ThreeVector momentum(0.0, 0.0, momentumZ    110       ThreeVector momentum(0.0, 0.0, momentumZ);
111       Particle *pb = new G4INCL::Particle(t, e    111       Particle *pb = new G4INCL::Particle(t, energy, momentum, ThreeVector());
112       PbarAtrestEntryChannel *obj = new PbarAt    112       PbarAtrestEntryChannel *obj = new PbarAtrestEntryChannel(theNucleus, pb); 
113       ParticleList fslist = obj->makeMesonStar    113       ParticleList fslist = obj->makeMesonStar();
114       const G4bool isProton = obj->ProtonIsThe    114       const G4bool isProton = obj->ProtonIsTheVictim();
115       delete pb;                                  115       delete pb;
116                                                   116 
117       //set Stopping time according to highest    117       //set Stopping time according to highest meson energy of the star
118       G4double temfin;                            118       G4double temfin;
119       G4double TLab;                              119       G4double TLab;
120       std::vector<G4double> energies;             120       std::vector<G4double> energies;
121       std::vector<G4double> projections;          121       std::vector<G4double> projections;
122       ThreeVector ab, cd;                         122       ThreeVector ab, cd;
123                                                   123 
124       for(ParticleIter pit = fslist.begin(), e    124       for(ParticleIter pit = fslist.begin(), e = fslist.end(); pit!=e; ++pit){
125         energies.push_back((*pit)->getKineticE    125         energies.push_back((*pit)->getKineticEnergy());
126         ab = (*pit)->boostVector();               126         ab = (*pit)->boostVector();
127         cd = (*pit)->getPosition();               127         cd = (*pit)->getPosition();
128         projections.push_back(ab.dot(cd)); //p    128         projections.push_back(ab.dot(cd)); //projection length
129       }// make vector of energies                 129       }// make vector of energies
130                                                   130       
131       temfin = 30.18 * std::pow(theNucleus->ge    131       temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
132       TLab = *max_element(energies.begin(), en    132       TLab = *max_element(energies.begin(), energies.end()); //choose max energy
133                                                   133 
134       // energy-dependent stopping time above     134       // energy-dependent stopping time above 2 AGeV
135       if(TLab>2000.)                              135       if(TLab>2000.)
136         temfin *= (5.8E4-TLab)/5.6E4;             136         temfin *= (5.8E4-TLab)/5.6E4;
137                                                   137 
138       maximumTime = temfin;                       138       maximumTime = temfin;  
139                                                   139 
140       // If the incoming particle is slow, use    140       // If the incoming particle is slow, use a larger stopping time
141       const G4double rMax = theNucleus->getUni    141       const G4double rMax = theNucleus->getUniverseRadius();
142       const G4double distance = 2.*rMax;          142       const G4double distance = 2.*rMax;
143       const G4double maxMesonVelocityProjectio    143       const G4double maxMesonVelocityProjection = *max_element(energies.begin(), energies.end());
144       const G4double traversalTime = distance     144       const G4double traversalTime = distance / maxMesonVelocityProjection;
145       if(maximumTime < traversalTime)             145       if(maximumTime < traversalTime)
146         maximumTime = traversalTime;              146         maximumTime = traversalTime;
147       INCL_DEBUG("Cascade stopping time is " <    147       INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
148                                                   148 
149                                                   149 
150       // Fill in the relevant kinematic variab    150       // Fill in the relevant kinematic variables
151       theNucleus->setIncomingAngularMomentum(G    151       theNucleus->setIncomingAngularMomentum(G4INCL::ThreeVector(0., 0., 0.));
152       theNucleus->setIncomingMomentum(G4INCL::    152       theNucleus->setIncomingMomentum(G4INCL::ThreeVector(0., 0., 0.));
153       if(isProton){                               153       if(isProton){
154         theNucleus->setInitialEnergy(pb->getMa    154         theNucleus->setInitialEnergy(pb->getMass()
155           + ParticleTable::getTableMass(theNuc    155           + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ() + 1,theNucleus->getS()));
156       }                                           156       }
157       else{                                       157       else{
158         theNucleus->setInitialEnergy(pb->getMa    158         theNucleus->setInitialEnergy(pb->getMass()
159           + ParticleTable::getTableMass(theNuc    159           + ParticleTable::getTableMass(theNucleus->getA() + 1,theNucleus->getZ(),theNucleus->getS()));
160       }                                           160       }
161       //kinetic energy excluded from the balan    161       //kinetic energy excluded from the balance
162                                                   162 
163       for(ParticleIter p = fslist.begin(), e =    163       for(ParticleIter p = fslist.begin(), e = fslist.end(); p!=e; ++p){
164         (*p)->makeProjectileSpectator();          164         (*p)->makeProjectileSpectator();                                        
165       }                                           165       }
166                                                   166       
167       generateAllAvatars();                       167       generateAllAvatars();
168       firstAvatar = false;                        168       firstAvatar = false;
169                                                   169 
170       // Get the entry avatars for mesons         170       // Get the entry avatars for mesons
171       IAvatarList theAvatarList = obj->bringMe    171       IAvatarList theAvatarList = obj->bringMesonStar(fslist, theNucleus);
172       delete obj;                                 172       delete obj;
173       theNucleus->getStore()->addParticleEntry    173       theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
174       INCL_DEBUG("Avatars added" << '\n');        174       INCL_DEBUG("Avatars added" << '\n');
175                                                   175       
176       return 99.;                                 176       return 99.;
177     }                                             177     } 
178                                                   178        
179 //D                                               179 //D
180                                                   180       
181     G4double StandardPropagationModel::shootPa    181     G4double StandardPropagationModel::shootParticle(ParticleType const type, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
182       theNucleus->setParticleNucleusCollision(    182       theNucleus->setParticleNucleusCollision();
183       currentTime = 0.0;                          183       currentTime = 0.0;
184                                                   184 
185       // Create the projectile particle           185       // Create the projectile particle
186       const G4double projectileMass = Particle    186       const G4double projectileMass = ParticleTable::getTableParticleMass(type);
187       G4double energy = kineticEnergy + projec    187       G4double energy = kineticEnergy + projectileMass;
188       G4double momentumZ = std::sqrt(energy*en    188       G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
189       ThreeVector momentum(0.0, 0.0, momentumZ    189       ThreeVector momentum(0.0, 0.0, momentumZ);
190       Particle *p= new G4INCL::Particle(type,     190       Particle *p= new G4INCL::Particle(type, energy, momentum, ThreeVector());
191                                                   191 
192       G4double temfin;                            192       G4double temfin;
193       G4double TLab;                              193       G4double TLab;
194       if( p->isMeson()) {                         194       if( p->isMeson()) {
195         temfin = 30.18 * std::pow(theNucleus->    195         temfin = 30.18 * std::pow(theNucleus->getA(), 0.17);
196         TLab = p->getKineticEnergy();             196         TLab = p->getKineticEnergy();
197       } else {                                    197       } else {
198         temfin = 29.8 * std::pow(theNucleus->g    198         temfin = 29.8 * std::pow(theNucleus->getA(), 0.16);
199         TLab = p->getKineticEnergy()/p->getA()    199         TLab = p->getKineticEnergy()/p->getA();
200       }                                           200       }
201                                                   201 
202       // energy-dependent stopping time above     202       // energy-dependent stopping time above 2 AGeV
203       if(TLab>2000.)                              203       if(TLab>2000.)
204         temfin *= (5.8E4-TLab)/5.6E4;             204         temfin *= (5.8E4-TLab)/5.6E4;
205                                                   205 
206       maximumTime = temfin;                       206       maximumTime = temfin;
207                                                   207 
208       // If the incoming particle is slow, use    208       // If the incoming particle is slow, use a larger stopping time
209       const G4double rMax = theNucleus->getUni    209       const G4double rMax = theNucleus->getUniverseRadius();
210       const G4double distance = 2.*rMax;          210       const G4double distance = 2.*rMax;
211       const G4double projectileVelocity = p->b    211       const G4double projectileVelocity = p->boostVector().mag();
212       const G4double traversalTime = distance     212       const G4double traversalTime = distance / projectileVelocity;
213       if(maximumTime < traversalTime)             213       if(maximumTime < traversalTime)
214         maximumTime = traversalTime;              214         maximumTime = traversalTime;
215       INCL_DEBUG("Cascade stopping time is " <    215       INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
216                                                   216 
217       // If Coulomb is activated, do not proce    217       // If Coulomb is activated, do not process events with impact
218       // parameter larger than the maximum imp    218       // parameter larger than the maximum impact parameter, taking into
219       // account Coulomb distortion.              219       // account Coulomb distortion.
220       if(impactParameter>CoulombDistortion::ma    220       if(impactParameter>CoulombDistortion::maxImpactParameter(p->getSpecies(), kineticEnergy, theNucleus)) {
221         INCL_DEBUG("impactParameter>CoulombDis    221         INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << '\n');
222         delete p;                                 222         delete p;
223         return -1.;                               223         return -1.;
224       }                                           224       }
225                                                   225 
226       ThreeVector position(impactParameter * s    226       ThreeVector position(impactParameter * std::cos(phi),
227           impactParameter * std::sin(phi),        227           impactParameter * std::sin(phi),
228           0.);                                    228           0.);
229       p->setPosition(position);                   229       p->setPosition(position);
230                                                   230 
231       // Fill in the relevant kinematic variab    231       // Fill in the relevant kinematic variables
232       theNucleus->setIncomingAngularMomentum(p    232       theNucleus->setIncomingAngularMomentum(p->getAngularMomentum());
233       theNucleus->setIncomingMomentum(p->getMo    233       theNucleus->setIncomingMomentum(p->getMomentum());
234       theNucleus->setInitialEnergy(p->getEnerg    234       theNucleus->setInitialEnergy(p->getEnergy()
235           + ParticleTable::getTableMass(theNuc    235           + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ(),theNucleus->getS()));
236                                                   236 
237       // Reset the particle kinematics to the     237       // Reset the particle kinematics to the INCL values
238       p->setINCLMass();                           238       p->setINCLMass();
239       p->setEnergy(p->getMass() + kineticEnerg    239       p->setEnergy(p->getMass() + kineticEnergy);
240       p->adjustMomentumFromEnergy();              240       p->adjustMomentumFromEnergy();
241                                                   241 
242       p->makeProjectileSpectator();               242       p->makeProjectileSpectator();
243       generateAllAvatars();                       243       generateAllAvatars();
244       firstAvatar = false;                        244       firstAvatar = false;
245                                                   245 
246       // Get the entry avatars from Coulomb an    246       // Get the entry avatars from Coulomb and put them in the Store
247       ParticleEntryAvatar *theEntryAvatar = Co    247       ParticleEntryAvatar *theEntryAvatar = CoulombDistortion::bringToSurface(p, theNucleus);
248       if(theEntryAvatar) {                        248       if(theEntryAvatar) {
249         theNucleus->getStore()->addParticleEnt    249         theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
250                                                   250 
251         return p->getTransversePosition().mag(    251         return p->getTransversePosition().mag();
252       } else {                                    252       } else {
253         delete p;                                 253         delete p;
254         return -1.;                               254         return -1.;
255       }                                           255       }
256     }                                             256     }
257                                                   257 
258     G4double StandardPropagationModel::shootCo    258     G4double StandardPropagationModel::shootComposite(ParticleSpecies const &species, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
259       theNucleus->setNucleusNucleusCollision()    259       theNucleus->setNucleusNucleusCollision();
260       currentTime = 0.0;                          260       currentTime = 0.0;
261                                                   261 
262       // Create the ProjectileRemnant object      262       // Create the ProjectileRemnant object
263       ProjectileRemnant *pr = new ProjectileRe    263       ProjectileRemnant *pr = new ProjectileRemnant(species, kineticEnergy);
264                                                   264 
265       // Same stopping time as for nucleon-nuc    265       // Same stopping time as for nucleon-nucleus
266       maximumTime = 29.8 * std::pow(theNucleus    266       maximumTime = 29.8 * std::pow(theNucleus->getA(), 0.16);
267                                                   267 
268       // If the incoming cluster is slow, use     268       // If the incoming cluster is slow, use a larger stopping time
269       const G4double rms = ParticleTable::getL    269       const G4double rms = ParticleTable::getLargestNuclearRadius(pr->getA(), pr->getZ());
270       const G4double rMax = theNucleus->getUni    270       const G4double rMax = theNucleus->getUniverseRadius();
271       const G4double distance = 2.*rMax + 2.72    271       const G4double distance = 2.*rMax + 2.725*rms;
272       const G4double projectileVelocity = pr->    272       const G4double projectileVelocity = pr->boostVector().mag();
273       const G4double traversalTime = distance     273       const G4double traversalTime = distance / projectileVelocity;
274       if(maximumTime < traversalTime)             274       if(maximumTime < traversalTime)
275         maximumTime = traversalTime;              275         maximumTime = traversalTime;
276       INCL_DEBUG("Cascade stopping time is " <    276       INCL_DEBUG("Cascade stopping time is " << maximumTime << '\n');
277                                                   277 
278       // If Coulomb is activated, do not proce    278       // If Coulomb is activated, do not process events with impact
279       // parameter larger than the maximum imp    279       // parameter larger than the maximum impact parameter, taking into
280       // account Coulomb distortion.              280       // account Coulomb distortion.
281       if(impactParameter>CoulombDistortion::ma    281       if(impactParameter>CoulombDistortion::maxImpactParameter(pr,theNucleus)) {
282         INCL_DEBUG("impactParameter>CoulombDis    282         INCL_DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << '\n');
283         delete pr;                                283         delete pr;
284         return -1.;                               284         return -1.;
285       }                                           285       }
286                                                   286 
287       // Position the cluster at the right imp    287       // Position the cluster at the right impact parameter
288       ThreeVector position(impactParameter * s    288       ThreeVector position(impactParameter * std::cos(phi),
289           impactParameter * std::sin(phi),        289           impactParameter * std::sin(phi),
290           0.);                                    290           0.);
291       pr->setPosition(position);                  291       pr->setPosition(position);
292                                                   292 
293       // Fill in the relevant kinematic variab    293       // Fill in the relevant kinematic variables
294       theNucleus->setIncomingAngularMomentum(p    294       theNucleus->setIncomingAngularMomentum(pr->getAngularMomentum());
295       theNucleus->setIncomingMomentum(pr->getM    295       theNucleus->setIncomingMomentum(pr->getMomentum());
296       theNucleus->setInitialEnergy(pr->getEner    296       theNucleus->setInitialEnergy(pr->getEnergy()
297           + ParticleTable::getTableMass(theNuc    297           + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ(),theNucleus->getS()));
298                                                   298 
299       generateAllAvatars();                       299       generateAllAvatars();
300       firstAvatar = false;                        300       firstAvatar = false;
301                                                   301 
302       // Get the entry avatars from Coulomb       302       // Get the entry avatars from Coulomb
303       IAvatarList theAvatarList                   303       IAvatarList theAvatarList
304         = CoulombDistortion::bringToSurface(pr    304         = CoulombDistortion::bringToSurface(pr, theNucleus);
305                                                   305 
306       if(theAvatarList.empty()) {                 306       if(theAvatarList.empty()) {
307         INCL_DEBUG("No ParticleEntryAvatar fou    307         INCL_DEBUG("No ParticleEntryAvatar found, transparent event" << '\n');
308         delete pr;                                308         delete pr;
309         return -1.;                               309         return -1.;
310       }                                           310       }
311                                                   311 
312       /* Store the internal kinematics of the     312       /* Store the internal kinematics of the projectile remnant.
313        *                                          313        *
314        * Note that this is at variance with th    314        * Note that this is at variance with the Fortran version, which stores
315        * the initial kinematics of the particl    315        * the initial kinematics of the particles *after* putting the spectators
316        * on mass shell, but *before* removing     316        * on mass shell, but *before* removing the same energy from the
317        * participants. Due to the different co    317        * participants. Due to the different code flow, doing so in the C++
318        * version leads to wrong excitation ene    318        * version leads to wrong excitation energies for the forced compound
319        * nucleus.                                 319        * nucleus.
320        */                                         320        */
321       pr->storeComponents();                      321       pr->storeComponents();
322                                                   322 
323       // Tell the Nucleus about the Projectile    323       // Tell the Nucleus about the ProjectileRemnant
324       theNucleus->setProjectileRemnant(pr);       324       theNucleus->setProjectileRemnant(pr);
325                                                   325 
326       // Register the ParticleEntryAvatars        326       // Register the ParticleEntryAvatars
327       theNucleus->getStore()->addParticleEntry    327       theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
328                                                   328 
329       return pr->getTransversePosition().mag()    329       return pr->getTransversePosition().mag();
330     }                                             330     }
331                                                   331 
332     G4double StandardPropagationModel::getStop    332     G4double StandardPropagationModel::getStoppingTime() {
333       return maximumTime;                         333       return maximumTime;
334     }                                             334     }
335                                                   335 
336     void StandardPropagationModel::setStopping    336     void StandardPropagationModel::setStoppingTime(G4double time) {
337 // assert(time>0.0);                              337 // assert(time>0.0);
338       maximumTime = time;                         338       maximumTime = time;
339     }                                             339     }
340                                                   340 
341     G4double StandardPropagationModel::getCurr    341     G4double StandardPropagationModel::getCurrentTime() {
342       return currentTime;                         342       return currentTime;
343     }                                             343     }
344                                                   344 
345     void StandardPropagationModel::setNucleus(    345     void StandardPropagationModel::setNucleus(G4INCL::Nucleus *nucleus)
346     {                                             346     {
347       theNucleus = nucleus;                       347       theNucleus = nucleus;
348     }                                             348     }
349                                                   349 
350     void StandardPropagationModel::registerAva    350     void StandardPropagationModel::registerAvatar(G4INCL::IAvatar *anAvatar)
351     {                                             351     {
352       if(anAvatar) theNucleus->getStore()->add    352       if(anAvatar) theNucleus->getStore()->add(anAvatar);
353     }                                             353     }
354                                                   354 
355     IAvatar *StandardPropagationModel::generat    355     IAvatar *StandardPropagationModel::generateBinaryCollisionAvatar(Particle * const p1, Particle * const p2) {
356       // Is either particle a participant?        356       // Is either particle a participant?
357       if(!p1->isParticipant() && !p2->isPartic    357       if(!p1->isParticipant() && !p2->isParticipant() && p1->getParticipantType()==p2->getParticipantType()) return NULL;
358                                                   358 
359       // Is it a pi-resonance collision (we do    359       // Is it a pi-resonance collision (we don't treat them)?
360       if((p1->isResonance() && p2->isPion()) |    360       if((p1->isResonance() && p2->isPion()) || (p1->isPion() && p2->isResonance()))
361         return NULL;                              361         return NULL;
362                                                   362 
363       // Is it a photon collision (we don't tr    363       // Is it a photon collision (we don't treat them)?
364       if(p1->isPhoton() || p2->isPhoton())        364       if(p1->isPhoton() || p2->isPhoton())
365         return NULL;                              365         return NULL;
366                                                   366 
367       // Will the avatar take place between no    367       // Will the avatar take place between now and the end of the cascade?
368       G4double minDistOfApproachSquared = 0.0;    368       G4double minDistOfApproachSquared = 0.0;
369       G4double t = getTime(p1, p2, &minDistOfA    369       G4double t = getTime(p1, p2, &minDistOfApproachSquared);
370       if(t>maximumTime || t<currentTime+hadron    370       if(t>maximumTime || t<currentTime+hadronizationTime) return NULL;
371                                                   371 
372       // Local energy. Jump through some hoops    372       // Local energy. Jump through some hoops to calculate the cross section
373       // at the collision point, and clean up     373       // at the collision point, and clean up after yourself afterwards.
374       G4bool hasLocalEnergy;                      374       G4bool hasLocalEnergy;
375       if(p1->isPion() || p2->isPion())            375       if(p1->isPion() || p2->isPion())
376         hasLocalEnergy = ((theLocalEnergyDelta    376         hasLocalEnergy = ((theLocalEnergyDeltaType == FirstCollisionLocalEnergy &&
377               theNucleus->getStore()->getBook(    377               theNucleus->getStore()->getBook().getAcceptedCollisions()==0) ||
378             theLocalEnergyDeltaType == AlwaysL    378             theLocalEnergyDeltaType == AlwaysLocalEnergy);
379       else                                        379       else
380         hasLocalEnergy = ((theLocalEnergyType     380         hasLocalEnergy = ((theLocalEnergyType == FirstCollisionLocalEnergy &&
381               theNucleus->getStore()->getBook(    381               theNucleus->getStore()->getBook().getAcceptedCollisions()==0) ||
382             theLocalEnergyType == AlwaysLocalE    382             theLocalEnergyType == AlwaysLocalEnergy);
383       const G4bool p1HasLocalEnergy = (hasLoca    383       const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->isMeson() && !p1->isAntiNucleon());
384       const G4bool p2HasLocalEnergy = (hasLoca    384       const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->isMeson() && !p2->isAntiNucleon());
385                                                   385 
386       if(p1HasLocalEnergy) {                      386       if(p1HasLocalEnergy) {
387         backupParticle1 = *p1;                    387         backupParticle1 = *p1;
388         p1->propagate(t - currentTime);           388         p1->propagate(t - currentTime);
389         if(p1->getPosition().mag() > theNucleu    389         if(p1->getPosition().mag() > theNucleus->getSurfaceRadius(p1)) {
390           *p1 = backupParticle1;                  390           *p1 = backupParticle1;
391           return NULL;                            391           return NULL;
392         }                                         392         }
393         KinematicsUtils::transformToLocalEnerg    393         KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p1);
394       }                                           394       } 
395       if(p2HasLocalEnergy) {                      395       if(p2HasLocalEnergy) {
396         backupParticle2 = *p2;                    396         backupParticle2 = *p2;
397         p2->propagate(t - currentTime);           397         p2->propagate(t - currentTime);
398         if(p2->getPosition().mag() > theNucleu    398         if(p2->getPosition().mag() > theNucleus->getSurfaceRadius(p2)) {
399           *p2 = backupParticle2;                  399           *p2 = backupParticle2;
400           if(p1HasLocalEnergy) {                  400           if(p1HasLocalEnergy) {
401             *p1 = backupParticle1;                401             *p1 = backupParticle1;
402           }                                       402           }
403           return NULL;                            403           return NULL;
404         }                                         404         }
405         KinematicsUtils::transformToLocalEnerg    405         KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p2);
406       }                                           406       }
407                                                   407 
408       // Compute the total cross section          408       // Compute the total cross section
409       const G4double totalCrossSection = Cross    409       const G4double totalCrossSection = CrossSections::total(p1, p2);
410       const G4double squareTotalEnergyInCM = K    410       const G4double squareTotalEnergyInCM = KinematicsUtils::squareTotalEnergyInCM(p1,p2);
411                                                   411 
412       // Restore particles to their state befo    412       // Restore particles to their state before the local-energy tweak
413       if(p1HasLocalEnergy) {                      413       if(p1HasLocalEnergy) {
414         *p1 = backupParticle1;                    414         *p1 = backupParticle1;
415       }                                           415       }
416       if(p2HasLocalEnergy) {                      416       if(p2HasLocalEnergy) {
417         *p2 = backupParticle2;                    417         *p2 = backupParticle2;
418       }                                           418       }
419                                                   419 
420       // Is the CM energy > cutNN? (no cutNN o    420       // Is the CM energy > cutNN? (no cutNN on the first collision)
421       if(theNucleus->getStore()->getBook().get    421       if(theNucleus->getStore()->getBook().getAcceptedCollisions()>0
422           && p1->isNucleon() && p2->isNucleon(    422           && p1->isNucleon() && p2->isNucleon()
423           && squareTotalEnergyInCM < BinaryCol    423           && squareTotalEnergyInCM < BinaryCollisionAvatar::getCutNNSquared()) return NULL;
424                                                   424 
425       // Do the particles come close enough to    425       // Do the particles come close enough to each other?
426       if(Math::tenPi*minDistOfApproachSquared     426       if(Math::tenPi*minDistOfApproachSquared > totalCrossSection) return NULL;
427                                                   427 
428       // Bomb out if the two collision partner    428       // Bomb out if the two collision partners are the same particle
429 // assert(p1->getID() != p2->getID());            429 // assert(p1->getID() != p2->getID());
430                                                   430 
431       // Return a new avatar, then!               431       // Return a new avatar, then!
432       return new G4INCL::BinaryCollisionAvatar    432       return new G4INCL::BinaryCollisionAvatar(t, totalCrossSection, theNucleus, p1, p2);
433     }                                             433     }
434                                                   434 
435     G4double StandardPropagationModel::getRefl    435     G4double StandardPropagationModel::getReflectionTime(G4INCL::Particle const * const aParticle) {
436       Intersection theIntersection(               436       Intersection theIntersection(
437           IntersectionFactory::getLaterTraject    437           IntersectionFactory::getLaterTrajectoryIntersection(
438             aParticle->getPosition(),             438             aParticle->getPosition(),
439             aParticle->getPropagationVelocity(    439             aParticle->getPropagationVelocity(),
440             theNucleus->getSurfaceRadius(aPart    440             theNucleus->getSurfaceRadius(aParticle)));
441       G4double time;                              441       G4double time;
442       if(theIntersection.exists) {                442       if(theIntersection.exists) {
443         time = currentTime + theIntersection.t    443         time = currentTime + theIntersection.time;
444       } else {                                    444       } else {
445         INCL_ERROR("Imaginary reflection time     445         INCL_ERROR("Imaginary reflection time for particle: " << '\n'
446           << aParticle->print());                 446           << aParticle->print());
447         time = 10000.0;                           447         time = 10000.0;
448       }                                           448       }
449       return time;                                449       return time;
450     }                                             450     }
451                                                   451 
452     G4double StandardPropagationModel::getTime    452     G4double StandardPropagationModel::getTime(G4INCL::Particle const * const particleA,
453                G4INCL::Particle const * const     453                G4INCL::Particle const * const particleB, G4double *minDistOfApproach) const
454     {                                             454     {
455       G4double time;                              455       G4double time;
456       G4INCL::ThreeVector t13 = particleA->get    456       G4INCL::ThreeVector t13 = particleA->getPropagationVelocity();
457       t13 -= particleB->getPropagationVelocity    457       t13 -= particleB->getPropagationVelocity();
458       G4INCL::ThreeVector distance = particleA    458       G4INCL::ThreeVector distance = particleA->getPosition();
459       distance -= particleB->getPosition();       459       distance -= particleB->getPosition();
460       const G4double t7 = t13.dot(distance);      460       const G4double t7 = t13.dot(distance);
461       const G4double dt = t13.mag2();             461       const G4double dt = t13.mag2();
462       if(dt <= 1.0e-10) {                         462       if(dt <= 1.0e-10) {
463         (*minDistOfApproach) = 100000.0;          463         (*minDistOfApproach) = 100000.0;
464         return currentTime + 100000.0;            464         return currentTime + 100000.0;
465       } else {                                    465       } else {
466         time = -t7/dt;                            466         time = -t7/dt;
467       }                                           467       }
468       (*minDistOfApproach) = distance.mag2() +    468       (*minDistOfApproach) = distance.mag2() + time * t7;
469       return currentTime + time;                  469       return currentTime + time;
470     }                                             470     }
471                                                   471 
472     void StandardPropagationModel::generateUpd    472     void StandardPropagationModel::generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles) {
473                                                   473 
474       // Loop over all the updated particles      474       // Loop over all the updated particles
475       for(ParticleIter updated=updatedParticle    475       for(ParticleIter updated=updatedParticles.begin(), e=updatedParticles.end(); updated!=e; ++updated)
476       {                                           476       {
477         // Loop over all the particles            477         // Loop over all the particles
478         for(ParticleIter particle=particles.be    478         for(ParticleIter particle=particles.begin(), end=particles.end(); particle!=end; ++particle)
479         {                                         479         {
480           /* Consider the generation of a coll    480           /* Consider the generation of a collision avatar only if (*particle)
481            * is not one of the updated particl    481            * is not one of the updated particles.
482            * The criterion makes sure that you    482            * The criterion makes sure that you don't generate avatars between
483            * updated particles. */                483            * updated particles. */
484           if(updatedParticles.contains(*partic    484           if(updatedParticles.contains(*particle)) continue;
485                                                   485 
486           registerAvatar(generateBinaryCollisi    486           registerAvatar(generateBinaryCollisionAvatar(*particle,*updated));
487         }                                         487         }
488       }                                           488       }
489     }                                             489     }
490                                                   490 
491     void StandardPropagationModel::generateCol    491     void StandardPropagationModel::generateCollisions(const ParticleList &particles) {
492       // Loop over all the particles              492       // Loop over all the particles
493       for(ParticleIter p1=particles.begin(), e    493       for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1) {
494         // Loop over the rest of the particles    494         // Loop over the rest of the particles
495         for(ParticleIter p2 = p1 + 1; p2 != pa    495         for(ParticleIter p2 = p1 + 1; p2 != particles.end(); ++p2) {
496           registerAvatar(generateBinaryCollisi    496           registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
497         }                                         497         }
498       }                                           498       }
499     }                                             499     }
500                                                   500 
501     void StandardPropagationModel::generateCol    501     void StandardPropagationModel::generateCollisions(const ParticleList &particles, const ParticleList &except) {
502                                                   502 
503       const G4bool haveExcept = !except.empty(    503       const G4bool haveExcept = !except.empty();
504                                                   504 
505       // Loop over all the particles              505       // Loop over all the particles
506       for(ParticleIter p1=particles.begin(), e    506       for(ParticleIter p1=particles.begin(), e=particles.end(); p1!=e; ++p1)
507       {                                           507       {
508         // Loop over the rest of the particles    508         // Loop over the rest of the particles
509         ParticleIter p2 = p1;                     509         ParticleIter p2 = p1;
510         for(++p2; p2 != particles.end(); ++p2)    510         for(++p2; p2 != particles.end(); ++p2)
511         {                                         511         {
512           // Skip the collision if both partic    512           // Skip the collision if both particles must be excluded
513           if(haveExcept && except.contains(*p1    513           if(haveExcept && except.contains(*p1) && except.contains(*p2)) continue;
514                                                   514 
515           registerAvatar(generateBinaryCollisi    515           registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
516         }                                         516         }
517       }                                           517       }
518                                                   518 
519     }                                             519     }
520                                                   520 
521     void StandardPropagationModel::updateAvata    521     void StandardPropagationModel::updateAvatars(const ParticleList &particles) {
522                                                   522 
523       for(ParticleIter iter=particles.begin(),    523       for(ParticleIter iter=particles.begin(), e=particles.end(); iter!=e; ++iter) {
524         G4double time = this->getReflectionTim    524         G4double time = this->getReflectionTime(*iter);
525         if(time <= maximumTime) registerAvatar    525         if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*iter, time, theNucleus));
526       }                                           526       }
527       ParticleList const &p = theNucleus->getS    527       ParticleList const &p = theNucleus->getStore()->getParticles();
528       generateUpdatedCollisions(particles, p);    528       generateUpdatedCollisions(particles, p);             // Predict collisions with spectators and participants
529     }                                             529     }
530                                                   530 
531     void StandardPropagationModel::generateAll    531     void StandardPropagationModel::generateAllAvatars() {
532       ParticleList const &particles = theNucle    532       ParticleList const &particles = theNucleus->getStore()->getParticles();
533 // assert(!particles.empty());                    533 // assert(!particles.empty());
534       G4double time;                              534       G4double time;
535       for(ParticleIter i=particles.begin(), e=    535       for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
536         time = this->getReflectionTime(*i);       536         time = this->getReflectionTime(*i);
537         if(time <= maximumTime) registerAvatar    537         if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
538       }                                           538       }
539       generateCollisions(particles);              539       generateCollisions(particles);
540       generateDecays(particles);                  540       generateDecays(particles);
541     }                                             541     }
542                                                   542 
543 #ifdef INCL_REGENERATE_AVATARS                    543 #ifdef INCL_REGENERATE_AVATARS
544     void StandardPropagationModel::generateAll    544     void StandardPropagationModel::generateAllAvatarsExceptUpdated(FinalState const * const fs) {
545       ParticleList const &particles = theNucle    545       ParticleList const &particles = theNucleus->getStore()->getParticles();
546 // assert(!particles.empty());                    546 // assert(!particles.empty());
547       for(ParticleIter i=particles.begin(), e=    547       for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
548         G4double time = this->getReflectionTim    548         G4double time = this->getReflectionTime(*i);
549         if(time <= maximumTime) registerAvatar    549         if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
550       }                                           550       }
551       ParticleList except = fs->getModifiedPar    551       ParticleList except = fs->getModifiedParticles();
552       ParticleList const &entering = fs->getEn    552       ParticleList const &entering = fs->getEnteringParticles();
553       except.insert(except.end(), entering.beg    553       except.insert(except.end(), entering.begin(), entering.end());
554       generateCollisions(particles,except);       554       generateCollisions(particles,except);
555       generateDecays(particles);                  555       generateDecays(particles);
556     }                                             556     }
557 #endif                                            557 #endif
558                                                   558 
559     void StandardPropagationModel::generateDec    559     void StandardPropagationModel::generateDecays(const ParticleList &particles) {
560       for(ParticleIter i=particles.begin(), e=    560       for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i) {
561          if((*i)->isDelta()) {                    561          if((*i)->isDelta()) {
562              G4double decayTime = DeltaDecayCh    562              G4double decayTime = DeltaDecayChannel::computeDecayTime((*i)); // time in fm/c
563            G4double time = currentTime + decay    563            G4double time = currentTime + decayTime;
564            if(time <= maximumTime) {              564            if(time <= maximumTime) {
565              registerAvatar(new DecayAvatar((*    565              registerAvatar(new DecayAvatar((*i), time, theNucleus));
566            }                                      566            }
567          }                                        567          }
568          else if((*i)->getType() == SigmaZero)    568          else if((*i)->getType() == SigmaZero) {
569            G4double decayTime = SigmaZeroDecay    569            G4double decayTime = SigmaZeroDecayChannel::computeDecayTime((*i)); // time in fm/c
570            G4double time = currentTime + decay    570            G4double time = currentTime + decayTime;
571            if(time <= maximumTime) {              571            if(time <= maximumTime) {
572              registerAvatar(new DecayAvatar((*    572              registerAvatar(new DecayAvatar((*i), time, theNucleus));
573            }                                      573            }
574          }                                        574          }
575         if((*i)->isOmega()) {                     575         if((*i)->isOmega()) {
576           G4double decayTimeOmega = PionResona    576           G4double decayTimeOmega = PionResonanceDecayChannel::computeDecayTime((*i));
577           G4double timeOmega = currentTime + d    577           G4double timeOmega = currentTime + decayTimeOmega;
578           if(timeOmega <= maximumTime) {          578           if(timeOmega <= maximumTime) {
579             registerAvatar(new DecayAvatar((*i    579             registerAvatar(new DecayAvatar((*i), timeOmega, theNucleus));
580           }                                       580           }
581         }                                         581         }
582       }                                           582       }
583     }                                             583     }
584                                                   584 
585     G4INCL::IAvatar* StandardPropagationModel:    585     G4INCL::IAvatar* StandardPropagationModel::propagate(FinalState const * const fs)
586     {                                             586     {
587       if(fs) {                                    587       if(fs) {
588         // We update only the information rela    588         // We update only the information related to particles that were updated
589         // by the previous avatar.                589         // by the previous avatar.
590 #ifdef INCL_REGENERATE_AVATARS                    590 #ifdef INCL_REGENERATE_AVATARS
591 #warning "The INCL_REGENERATE_AVATARS code has    591 #warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
592         if(!fs->getModifiedParticles().empty()    592         if(!fs->getModifiedParticles().empty() || !fs->getEnteringParticles().empty() || !fs->getCreatedParticles().empty()) {
593           // Regenerates the entire avatar lis    593           // Regenerates the entire avatar list, skipping collisions between
594           // updated particles                    594           // updated particles
595           theNucleus->getStore()->clearAvatars    595           theNucleus->getStore()->clearAvatars();
596           theNucleus->getStore()->initialisePa    596           theNucleus->getStore()->initialiseParticleAvatarConnections();
597           generateAllAvatarsExceptUpdated(fs);    597           generateAllAvatarsExceptUpdated(fs);
598         }                                         598         }
599 #else                                             599 #else
600         ParticleList const &updatedParticles =    600         ParticleList const &updatedParticles = fs->getModifiedParticles();
601         if(fs->getValidity()==PauliBlockedFS)     601         if(fs->getValidity()==PauliBlockedFS) {
602           // This final state might represents    602           // This final state might represents the outcome of a Pauli-blocked delta
603           // decay                                603           // decay
604 // assert(updatedParticles.empty() || (updated    604 // assert(updatedParticles.empty() || (updatedParticles.size()==1 && updatedParticles.front()->isResonance()));
605 // assert(fs->getEnteringParticles().empty() &    605 // assert(fs->getEnteringParticles().empty() && fs->getCreatedParticles().empty() && fs->getOutgoingParticles().empty() && fs->getDestroyedParticles().empty());
606           generateDecays(updatedParticles);       606           generateDecays(updatedParticles);
607         } else {                                  607         } else {
608           ParticleList const &entering = fs->g    608           ParticleList const &entering = fs->getEnteringParticles();
609           generateDecays(updatedParticles);       609           generateDecays(updatedParticles);
610           generateDecays(entering);               610           generateDecays(entering);
611                                                   611 
612           ParticleList const &created = fs->ge    612           ParticleList const &created = fs->getCreatedParticles();
613           if(created.empty() && entering.empty    613           if(created.empty() && entering.empty())
614             updateAvatars(updatedParticles);      614             updateAvatars(updatedParticles);
615           else {                                  615           else {
616             ParticleList updatedParticlesCopy     616             ParticleList updatedParticlesCopy = updatedParticles;
617             updatedParticlesCopy.insert(update    617             updatedParticlesCopy.insert(updatedParticlesCopy.end(), entering.begin(), entering.end());
618             updatedParticlesCopy.insert(update    618             updatedParticlesCopy.insert(updatedParticlesCopy.end(), created.begin(), created.end());
619             updateAvatars(updatedParticlesCopy    619             updateAvatars(updatedParticlesCopy);
620           }                                       620           }
621         }                                         621         }
622 #endif                                            622 #endif
623       }                                           623       }
624                                                   624 
625       G4INCL::IAvatar *theAvatar = theNucleus-    625       G4INCL::IAvatar *theAvatar = theNucleus->getStore()->findSmallestTime();
626       if(theAvatar == 0) return 0; // Avatar l    626       if(theAvatar == 0) return 0; // Avatar list is empty
627       //      theAvatar->dispose();               627       //      theAvatar->dispose();
628                                                   628 
629       if(theAvatar->getTime() < currentTime) {    629       if(theAvatar->getTime() < currentTime) {
630         INCL_ERROR("Avatar time = " << theAvat    630         INCL_ERROR("Avatar time = " << theAvatar->getTime() << ", currentTime = " << currentTime << '\n');
631         return 0;                                 631         return 0;
632       } else if(theAvatar->getTime() > current    632       } else if(theAvatar->getTime() > currentTime) {
633         theNucleus->getStore()->timeStep(theAv    633         theNucleus->getStore()->timeStep(theAvatar->getTime() - currentTime);
634                                                   634 
635         currentTime = theAvatar->getTime();       635         currentTime = theAvatar->getTime();
636         theNucleus->getStore()->getBook().setC    636         theNucleus->getStore()->getBook().setCurrentTime(currentTime);
637       }                                           637       }
638                                                   638 
639       return theAvatar;                           639       return theAvatar;
640     }                                             640     }
641 }                                                 641 }
642                                                   642