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