Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // INCL++ intra-nuclear cascade model 27 // Alain Boudard, CEA-Saclay, France 28 // Joseph Cugnon, University of Liege, Belgium 29 // Jean-Christophe David, CEA-Saclay, France 30 // Pekka Kaitaniemi, CEA-Saclay, France, and H 31 // Sylvie Leray, CEA-Saclay, France 32 // Davide Mancusi, CEA-Saclay, France 33 // 34 #define INCLXX_IN_GEANT4_MODE 1 35 36 #include "globals.hh" 37 38 /** \file G4INCLProjectileRemnant.cc 39 * \brief Class for constructing a projectile- 40 * 41 * \date 20 March 2012 42 * \author Davide Mancusi 43 */ 44 45 #include "G4INCLProjectileRemnant.hh" 46 #include <algorithm> 47 #include <numeric> 48 49 namespace G4INCL { 50 51 void ProjectileRemnant::reset() { 52 deleteParticles(); 53 thePosition = ThreeVector(); 54 theMomentum = ThreeVector(); 55 theEnergy = 0.0; 56 thePotentialEnergy = 0.0; 57 theA = 0; 58 theZ = 0; 59 nCollisions = 0; 60 61 for(std::map<long, Particle*>::const_itera 62 Particle *p = new Particle(*(i->second)) 63 EnergyLevelMap::iterator energyIter = th 64 // assert(energyIter!=theInitialEnergyLevels.e 65 const G4double energyLevel = energyIter- 66 theInitialEnergyLevels.erase(energyIter) 67 theInitialEnergyLevels[p->getID()] = ene 68 addParticle(p); 69 } 70 if(theA>0) 71 thePosition /= theA; 72 setTableMass(); 73 INCL_DEBUG("ProjectileRemnant object was r 74 } 75 76 void ProjectileRemnant::removeParticle(Parti 77 // assert(p->isNucleon() || p->isLambda()); 78 79 INCL_DEBUG("The following Particle is abou 80 << '\n' << p->print() 81 << "theProjectileCorrection=" << thePr 82 // Update A, Z, S, momentum, and energy of 83 theA -= p->getA(); 84 theZ -= p->getZ(); 85 theS -= p->getS(); 86 87 ThreeVector const &oldMomentum = p->getMom 88 const G4double oldEnergy = p->getEnergy(); 89 Cluster::removeParticle(p); 90 91 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA 92 ThreeVector theTotalMomentum; 93 G4double theTotalEnergy = 0.; 94 const G4double theThreshold = 0.1; 95 #endif 96 97 if(getA()>0) { // if there are any particl 98 // assert((unsigned int)getA()==particles.size 99 100 const G4double theProjectileCorrectionPe 101 102 // Update the kinematics of the componen 103 for(ParticleIter i=particles.begin(), e= 104 (*i)->setEnergy((*i)->getEnergy() + th 105 (*i)->setMass((*i)->getInvariantMass() 106 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEA 107 theTotalMomentum += (*i)->getMomentum( 108 theTotalEnergy += (*i)->getEnergy(); 109 #endif 110 } 111 } 112 113 theMomentum -= oldMomentum; 114 theEnergy -= oldEnergy - theProjectileCorr 115 116 // assert(std::abs((theTotalMomentum-theMoment 117 // assert(std::abs(theTotalEnergy-theEnergy)<t 118 INCL_DEBUG("After Particle removal, the Pr 119 << '\n' << print()); 120 } 121 122 ParticleList ProjectileRemnant::addDynamical 123 // Try as hard as possible to add back all 124 // Don't add spectators that lead to negat 125 // iterate over the spectators as many tim 126 // absolutely sure that all of them were r 127 unsigned int accepted; 128 unsigned long loopCounter = 0; 129 const unsigned long maxLoopCounter = 10000 130 do { 131 accepted = 0; 132 ParticleList toBeAdded = pL; 133 for(ParticleIter p=toBeAdded.begin(), e= 134 G4bool isAccepted = addDynamicalSpecta 135 if(isAccepted) { 136 pL.remove(*p); 137 accepted++; 138 } 139 } 140 ++loopCounter; 141 } while(loopCounter<maxLoopCounter && acce 142 return pL; 143 } 144 145 ParticleList ProjectileRemnant::addAllDynami 146 // Put all the spectators in the projectil 147 ThreeVector theNewMomentum = theMomentum; 148 G4double theNewEnergy = theEnergy; 149 G4int theNewA = theA; 150 G4int theNewZ = theZ; 151 G4int theNewS = theS; 152 for(ParticleIter p=pL.begin(), e=pL.end(); 153 // assert((*p)->isNucleonorLambda()); 154 // Add the initial (off-shell) momentum 155 theNewMomentum += getStoredMomentum(*p); 156 theNewEnergy += (*p)->getEnergy(); 157 theNewA += (*p)->getA(); 158 theNewZ += (*p)->getZ(); 159 theNewS += (*p)->getS(); 160 } 161 162 // Check that the excitation energy of the 163 const G4double theNewMass = ParticleTable: 164 const G4double theNewExcitationEnergy = co 165 const G4double theNewEffectiveMass = theNe 166 167 // If this condition is satisfied, there i 168 // "most" method 169 if(theNewEnergy<theNewEffectiveMass) { 170 INCL_WARN("Could not add all the dynamic 171 << " Falling back to the \"most\" m 172 return addMostDynamicalSpectators(pL); 173 } 174 175 // Add all the participants to the project 176 for(ParticleIter p=pL.begin(), e=pL.end(); 177 particles.push_back(*p); 178 } 179 180 // Rescale the momentum of the projectile 181 // correct value 182 const G4double scalingFactorSquared = (the 183 const G4double scalingFactor = std::sqrt(s 184 INCL_DEBUG("Scaling factor for the project 185 186 theA = theNewA; 187 theZ = theNewZ; 188 theS = theNewS; 189 theMomentum = theNewMomentum * scalingFact 190 theEnergy = theNewEnergy; 191 192 return ParticleList(); 193 } 194 195 ParticleList ProjectileRemnant::addMostDynam 196 // Try as hard as possible to add back all 197 // Don't add spectators that lead to negat 198 // adding all of them, and repeatedly remo 199 // the excitation energy becomes non-negat 200 201 // Put all the spectators in the projectil 202 ThreeVector theNewMomentum = theMomentum; 203 G4double theNewEnergy = theEnergy; 204 G4int theNewA = theA; 205 G4int theNewZ = theZ; 206 G4int theNewS = theS; 207 for(ParticleIter p=pL.begin(), e=pL.end(); 208 // assert((*p)->isNucleonorLambda()); 209 // Add the initial (off-shell) momentum 210 theNewMomentum += getStoredMomentum(*p); 211 theNewEnergy += (*p)->getEnergy(); 212 theNewA += (*p)->getA(); 213 theNewZ += (*p)->getZ(); 214 theNewS += (*p)->getS(); 215 } 216 217 // Check that the excitation energy of the 218 const G4double theNewMass = ParticleTable: 219 const G4double theNewInvariantMassSquared 220 221 G4bool positiveExcitationEnergy = false; 222 if(theNewInvariantMassSquared>=0.) { 223 const G4double theNewInvariantMass = std 224 positiveExcitationEnergy = (theNewInvari 225 } 226 227 // Keep removing nucleons from the project 228 // non-negative excitation energy. 229 ParticleList rejected; 230 while(!positiveExcitationEnergy && !pL.emp 231 G4double maxExcitationEnergy = -1.E30; 232 ParticleMutableIter best = pL.end(); 233 ThreeVector bestMomentum; 234 G4double bestEnergy = -1.; 235 G4int bestA = -1, bestZ = -1, bestS = 0; 236 for(ParticleList::iterator p=pL.begin(), 237 // Subtract the initial (off-shell) mo 238 // projectile remnant 239 const ThreeVector theNewerMomentum = t 240 const G4double theNewerEnergy = theNew 241 const G4int theNewerA = theNewA - (*p) 242 const G4int theNewerZ = theNewZ - (*p) 243 const G4int theNewerS = theNewS - (*p) 244 245 const G4double theNewerMass = Particle 246 const G4double theNewerInvariantMassSq 247 248 if(theNewerInvariantMassSquared>=-1.e- 249 const G4double theNewerInvariantMass 250 const G4double theNewerExcitationEne 251 // Pick the nucleon that maximises t 252 // ProjectileRemnant 253 if(theNewerExcitationEnergy>maxExcit 254 best = p; 255 maxExcitationEnergy = theNewerExci 256 bestMomentum = theNewerMomentum; 257 bestEnergy = theNewerEnergy; 258 bestA = theNewerA; 259 bestZ = theNewerZ; 260 bestS = theNewerS; 261 } 262 } 263 } 264 265 // If we couldn't even calculate the exc 266 if(best==pL.end()) 267 return pL; 268 269 rejected.push_back(*best); 270 pL.erase(best); 271 theNewMomentum = bestMomentum; 272 theNewEnergy = bestEnergy; 273 theNewA = bestA; 274 theNewZ = bestZ; 275 theNewS = bestS; 276 277 if(maxExcitationEnergy>0.) { 278 // Stop here 279 positiveExcitationEnergy = true; 280 } 281 } 282 283 // Add the accepted participants to the pr 284 for(ParticleIter p=pL.begin(), e=pL.end(); 285 particles.push_back(*p); 286 } 287 theA = theNewA; 288 theZ = theNewZ; 289 theS = theNewS; 290 theMomentum = theNewMomentum; 291 theEnergy = theNewEnergy; 292 293 return rejected; 294 } 295 296 G4bool ProjectileRemnant::addDynamicalSpecta 297 // assert(p->isNucleon()); 298 299 // Add the initial (off-shell) momentum an 300 ThreeVector const &oldMomentum = getStored 301 const ThreeVector theNewMomentum = theMome 302 const G4double oldEnergy = p->getEnergy(); 303 const G4double theNewEnergy = theEnergy + 304 305 // Check that the excitation energy of the 306 const G4double theNewMass = ParticleTable: 307 const G4double theNewInvariantMassSquared 308 309 if(theNewInvariantMassSquared<0.) 310 return false; 311 312 const G4double theNewInvariantMass = std:: 313 314 if(theNewInvariantMass-theNewMass<-1.e-5) 315 return false; // negative excitation ene 316 317 // Add the spectator to the projectile rem 318 theA += p->getA(); 319 theZ += p->getZ(); 320 theMomentum = theNewMomentum; 321 theEnergy = theNewEnergy; 322 particles.push_back(p); 323 return true; 324 } 325 326 G4double ProjectileRemnant::computeExcitatio 327 const EnergyLevels theEnergyLevels = getPr 328 return computeExcitationEnergy(theEnergyLe 329 } 330 331 G4double ProjectileRemnant::computeExcitatio 332 const EnergyLevels theEnergyLevels = getPr 333 return computeExcitationEnergy(theEnergyLe 334 } 335 336 G4double ProjectileRemnant::computeExcitatio 337 // The ground-state energy is the sum of t 338 // energies. 339 // For the last nucleon, return 0 so that 340 // shell. 341 const std::size_t theNewA = levels.size(); 342 // assert(theNewA>0); 343 if(theNewA==1) 344 return 0.; 345 346 const G4double groundState = theGroundStat 347 348 // Compute the sum of the presently occupi 349 const G4double excitedState = std::accumul 350 levels.cbegin(), 351 levels.cend(), 352 0.); 353 354 return excitedState-groundState; 355 } 356 357 ProjectileRemnant::EnergyLevels ProjectileRe 358 EnergyLevels theEnergyLevels; 359 for(ParticleIter p=particles.begin(), e=pa 360 if((*p)->getID()!=exceptID) { 361 EnergyLevelMap::const_iterator i = the 362 // assert(i!=theInitialEnergyLevels.end()); 363 theEnergyLevels.push_back(i->second); 364 } 365 } 366 // assert(theEnergyLevels.size()==particles.si 367 return theEnergyLevels; 368 } 369 370 ProjectileRemnant::EnergyLevels ProjectileRe 371 EnergyLevels theEnergyLevels; 372 for(ParticleIter p=particles.begin(), e=pa 373 EnergyLevelMap::const_iterator i = theIn 374 // assert(i!=theInitialEnergyLevels.end()); 375 theEnergyLevels.push_back(i->second); 376 } 377 for(ParticleIter p=pL.begin(), e=pL.end(); 378 EnergyLevelMap::const_iterator i = theIn 379 // assert(i!=theInitialEnergyLevels.end()); 380 theEnergyLevels.push_back(i->second); 381 } 382 383 // assert(theEnergyLevels.size()==particles.si 384 return theEnergyLevels; 385 } 386 387 } 388 389