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