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 /* 39 * StandardPropagationModel.cpp 40 * 41 * \date 4 juin 2009 42 * \author Pekka Kaitaniemi 43 */ 44 45 #include "G4INCLStandardPropagationModel.hh" 46 #include "G4INCLPbarAtrestEntryChannel.hh" 47 #include "G4INCLSurfaceAvatar.hh" 48 #include "G4INCLBinaryCollisionAvatar.hh" 49 #include "G4INCLDecayAvatar.hh" 50 #include "G4INCLCrossSections.hh" 51 #include "G4INCLRandom.hh" 52 #include <iostream> 53 #include <algorithm> 54 #include "G4INCLLogger.hh" 55 #include "G4INCLGlobals.hh" 56 #include "G4INCLKinematicsUtils.hh" 57 #include "G4INCLCoulombDistortion.hh" 58 #include "G4INCLDeltaDecayChannel.hh" 59 #include "G4INCLSigmaZeroDecayChannel.hh" 60 #include "G4INCLPionResonanceDecayChannel.hh" 61 #include "G4INCLParticleEntryAvatar.hh" 62 #include "G4INCLIntersection.hh" 63 #include <vector> 64 65 namespace G4INCL { 66 67 StandardPropagationModel::StandardPropagat 68 :theNucleus(0), maximumTime(70.0), curre 69 hadronizationTime(hTime), 70 firstAvatar(true), 71 theLocalEnergyType(localEnergyType), 72 theLocalEnergyDeltaType(localEnergyDelta 73 { 74 } 75 76 StandardPropagationModel::~StandardPropaga 77 { 78 delete theNucleus; 79 } 80 81 G4INCL::Nucleus* StandardPropagationModel: 82 { 83 return theNucleus; 84 } 85 86 //D 87 88 G4double StandardPropagationModel::shoot(P 89 if(projectileSpecies.theType==Composite) 90 return shootComposite(projectileSpecie 91 } 92 else if(projectileSpecies.theType==antiP 93 return shootAtrest(projectileSpecies.t 94 } 95 else{ 96 return shootParticle(projectileSpecies 97 } 98 } 99 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 182 theNucleus->setParticleNucleusCollision( 183 currentTime = 0.0; 184 185 // Create the projectile particle 186 const G4double projectileMass = Particle 187 G4double energy = kineticEnergy + projec 188 G4double momentumZ = std::sqrt(energy*en 189 ThreeVector momentum(0.0, 0.0, momentumZ 190 Particle *p= new G4INCL::Particle(type, 191 192 G4double temfin; 193 G4double TLab; 194 if( p->isMeson()) { 195 temfin = 30.18 * std::pow(theNucleus-> 196 TLab = p->getKineticEnergy(); 197 } else { 198 temfin = 29.8 * std::pow(theNucleus->g 199 TLab = p->getKineticEnergy()/p->getA() 200 } 201 202 // energy-dependent stopping time above 203 if(TLab>2000.) 204 temfin *= (5.8E4-TLab)/5.6E4; 205 206 maximumTime = temfin; 207 208 // If the incoming particle is slow, use 209 const G4double rMax = theNucleus->getUni 210 const G4double distance = 2.*rMax; 211 const G4double projectileVelocity = p->b 212 const G4double traversalTime = distance 213 if(maximumTime < traversalTime) 214 maximumTime = traversalTime; 215 INCL_DEBUG("Cascade stopping time is " < 216 217 // If Coulomb is activated, do not proce 218 // parameter larger than the maximum imp 219 // account Coulomb distortion. 220 if(impactParameter>CoulombDistortion::ma 221 INCL_DEBUG("impactParameter>CoulombDis 222 delete p; 223 return -1.; 224 } 225 226 ThreeVector position(impactParameter * s 227 impactParameter * std::sin(phi), 228 0.); 229 p->setPosition(position); 230 231 // Fill in the relevant kinematic variab 232 theNucleus->setIncomingAngularMomentum(p 233 theNucleus->setIncomingMomentum(p->getMo 234 theNucleus->setInitialEnergy(p->getEnerg 235 + ParticleTable::getTableMass(theNuc 236 237 // Reset the particle kinematics to the 238 p->setINCLMass(); 239 p->setEnergy(p->getMass() + kineticEnerg 240 p->adjustMomentumFromEnergy(); 241 242 p->makeProjectileSpectator(); 243 generateAllAvatars(); 244 firstAvatar = false; 245 246 // Get the entry avatars from Coulomb an 247 ParticleEntryAvatar *theEntryAvatar = Co 248 if(theEntryAvatar) { 249 theNucleus->getStore()->addParticleEnt 250 251 return p->getTransversePosition().mag( 252 } else { 253 delete p; 254 return -1.; 255 } 256 } 257 258 G4double StandardPropagationModel::shootCo 259 theNucleus->setNucleusNucleusCollision() 260 currentTime = 0.0; 261 262 // Create the ProjectileRemnant object 263 ProjectileRemnant *pr = new ProjectileRe 264 265 // Same stopping time as for nucleon-nuc 266 maximumTime = 29.8 * std::pow(theNucleus 267 268 // If the incoming cluster is slow, use 269 const G4double rms = ParticleTable::getL 270 const G4double rMax = theNucleus->getUni 271 const G4double distance = 2.*rMax + 2.72 272 const G4double projectileVelocity = pr-> 273 const G4double traversalTime = distance 274 if(maximumTime < traversalTime) 275 maximumTime = traversalTime; 276 INCL_DEBUG("Cascade stopping time is " < 277 278 // If Coulomb is activated, do not proce 279 // parameter larger than the maximum imp 280 // account Coulomb distortion. 281 if(impactParameter>CoulombDistortion::ma 282 INCL_DEBUG("impactParameter>CoulombDis 283 delete pr; 284 return -1.; 285 } 286 287 // Position the cluster at the right imp 288 ThreeVector position(impactParameter * s 289 impactParameter * std::sin(phi), 290 0.); 291 pr->setPosition(position); 292 293 // Fill in the relevant kinematic variab 294 theNucleus->setIncomingAngularMomentum(p 295 theNucleus->setIncomingMomentum(pr->getM 296 theNucleus->setInitialEnergy(pr->getEner 297 + ParticleTable::getTableMass(theNuc 298 299 generateAllAvatars(); 300 firstAvatar = false; 301 302 // Get the entry avatars from Coulomb 303 IAvatarList theAvatarList 304 = CoulombDistortion::bringToSurface(pr 305 306 if(theAvatarList.empty()) { 307 INCL_DEBUG("No ParticleEntryAvatar fou 308 delete pr; 309 return -1.; 310 } 311 312 /* Store the internal kinematics of the 313 * 314 * Note that this is at variance with th 315 * the initial kinematics of the particl 316 * on mass shell, but *before* removing 317 * participants. Due to the different co 318 * version leads to wrong excitation ene 319 * nucleus. 320 */ 321 pr->storeComponents(); 322 323 // Tell the Nucleus about the Projectile 324 theNucleus->setProjectileRemnant(pr); 325 326 // Register the ParticleEntryAvatars 327 theNucleus->getStore()->addParticleEntry 328 329 return pr->getTransversePosition().mag() 330 } 331 332 G4double StandardPropagationModel::getStop 333 return maximumTime; 334 } 335 336 void StandardPropagationModel::setStopping 337 // assert(time>0.0); 338 maximumTime = time; 339 } 340 341 G4double StandardPropagationModel::getCurr 342 return currentTime; 343 } 344 345 void StandardPropagationModel::setNucleus( 346 { 347 theNucleus = nucleus; 348 } 349 350 void StandardPropagationModel::registerAva 351 { 352 if(anAvatar) theNucleus->getStore()->add 353 } 354 355 IAvatar *StandardPropagationModel::generat 356 // Is either particle a participant? 357 if(!p1->isParticipant() && !p2->isPartic 358 359 // Is it a pi-resonance collision (we do 360 if((p1->isResonance() && p2->isPion()) | 361 return NULL; 362 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 368 G4double minDistOfApproachSquared = 0.0; 369 G4double t = getTime(p1, p2, &minDistOfA 370 if(t>maximumTime || t<currentTime+hadron 371 372 // Local energy. Jump through some hoops 373 // at the collision point, and clean up 374 G4bool hasLocalEnergy; 375 if(p1->isPion() || p2->isPion()) 376 hasLocalEnergy = ((theLocalEnergyDelta 377 theNucleus->getStore()->getBook( 378 theLocalEnergyDeltaType == AlwaysL 379 else 380 hasLocalEnergy = ((theLocalEnergyType 381 theNucleus->getStore()->getBook( 382 theLocalEnergyType == AlwaysLocalE 383 const G4bool p1HasLocalEnergy = (hasLoca 384 const G4bool p2HasLocalEnergy = (hasLoca 385 386 if(p1HasLocalEnergy) { 387 backupParticle1 = *p1; 388 p1->propagate(t - currentTime); 389 if(p1->getPosition().mag() > theNucleu 390 *p1 = backupParticle1; 391 return NULL; 392 } 393 KinematicsUtils::transformToLocalEnerg 394 } 395 if(p2HasLocalEnergy) { 396 backupParticle2 = *p2; 397 p2->propagate(t - currentTime); 398 if(p2->getPosition().mag() > theNucleu 399 *p2 = backupParticle2; 400 if(p1HasLocalEnergy) { 401 *p1 = backupParticle1; 402 } 403 return NULL; 404 } 405 KinematicsUtils::transformToLocalEnerg 406 } 407 408 // Compute the total cross section 409 const G4double totalCrossSection = Cross 410 const G4double squareTotalEnergyInCM = K 411 412 // Restore particles to their state befo 413 if(p1HasLocalEnergy) { 414 *p1 = backupParticle1; 415 } 416 if(p2HasLocalEnergy) { 417 *p2 = backupParticle2; 418 } 419 420 // Is the CM energy > cutNN? (no cutNN o 421 if(theNucleus->getStore()->getBook().get 422 && p1->isNucleon() && p2->isNucleon( 423 && squareTotalEnergyInCM < BinaryCol 424 425 // Do the particles come close enough to 426 if(Math::tenPi*minDistOfApproachSquared 427 428 // Bomb out if the two collision partner 429 // assert(p1->getID() != p2->getID()); 430 431 // Return a new avatar, then! 432 return new G4INCL::BinaryCollisionAvatar 433 } 434 435 G4double StandardPropagationModel::getRefl 436 Intersection theIntersection( 437 IntersectionFactory::getLaterTraject 438 aParticle->getPosition(), 439 aParticle->getPropagationVelocity( 440 theNucleus->getSurfaceRadius(aPart 441 G4double time; 442 if(theIntersection.exists) { 443 time = currentTime + theIntersection.t 444 } else { 445 INCL_ERROR("Imaginary reflection time 446 << aParticle->print()); 447 time = 10000.0; 448 } 449 return time; 450 } 451 452 G4double StandardPropagationModel::getTime 453 G4INCL::Particle const * const 454 { 455 G4double time; 456 G4INCL::ThreeVector t13 = particleA->get 457 t13 -= particleB->getPropagationVelocity 458 G4INCL::ThreeVector distance = particleA 459 distance -= particleB->getPosition(); 460 const G4double t7 = t13.dot(distance); 461 const G4double dt = t13.mag2(); 462 if(dt <= 1.0e-10) { 463 (*minDistOfApproach) = 100000.0; 464 return currentTime + 100000.0; 465 } else { 466 time = -t7/dt; 467 } 468 (*minDistOfApproach) = distance.mag2() + 469 return currentTime + time; 470 } 471 472 void StandardPropagationModel::generateUpd 473 474 // Loop over all the updated particles 475 for(ParticleIter updated=updatedParticle 476 { 477 // Loop over all the particles 478 for(ParticleIter particle=particles.be 479 { 480 /* Consider the generation of a coll 481 * is not one of the updated particl 482 * The criterion makes sure that you 483 * updated particles. */ 484 if(updatedParticles.contains(*partic 485 486 registerAvatar(generateBinaryCollisi 487 } 488 } 489 } 490 491 void StandardPropagationModel::generateCol 492 // Loop over all the particles 493 for(ParticleIter p1=particles.begin(), e 494 // Loop over the rest of the particles 495 for(ParticleIter p2 = p1 + 1; p2 != pa 496 registerAvatar(generateBinaryCollisi 497 } 498 } 499 } 500 501 void StandardPropagationModel::generateCol 502 503 const G4bool haveExcept = !except.empty( 504 505 // Loop over all the particles 506 for(ParticleIter p1=particles.begin(), e 507 { 508 // Loop over the rest of the particles 509 ParticleIter p2 = p1; 510 for(++p2; p2 != particles.end(); ++p2) 511 { 512 // Skip the collision if both partic 513 if(haveExcept && except.contains(*p1 514 515 registerAvatar(generateBinaryCollisi 516 } 517 } 518 519 } 520 521 void StandardPropagationModel::updateAvata 522 523 for(ParticleIter iter=particles.begin(), 524 G4double time = this->getReflectionTim 525 if(time <= maximumTime) registerAvatar 526 } 527 ParticleList const &p = theNucleus->getS 528 generateUpdatedCollisions(particles, p); 529 } 530 531 void StandardPropagationModel::generateAll 532 ParticleList const &particles = theNucle 533 // assert(!particles.empty()); 534 G4double time; 535 for(ParticleIter i=particles.begin(), e= 536 time = this->getReflectionTime(*i); 537 if(time <= maximumTime) registerAvatar 538 } 539 generateCollisions(particles); 540 generateDecays(particles); 541 } 542 543 #ifdef INCL_REGENERATE_AVATARS 544 void StandardPropagationModel::generateAll 545 ParticleList const &particles = theNucle 546 // assert(!particles.empty()); 547 for(ParticleIter i=particles.begin(), e= 548 G4double time = this->getReflectionTim 549 if(time <= maximumTime) registerAvatar 550 } 551 ParticleList except = fs->getModifiedPar 552 ParticleList const &entering = fs->getEn 553 except.insert(except.end(), entering.beg 554 generateCollisions(particles,except); 555 generateDecays(particles); 556 } 557 #endif 558 559 void StandardPropagationModel::generateDec 560 for(ParticleIter i=particles.begin(), e= 561 if((*i)->isDelta()) { 562 G4double decayTime = DeltaDecayCh 563 G4double time = currentTime + decay 564 if(time <= maximumTime) { 565 registerAvatar(new DecayAvatar((* 566 } 567 } 568 else if((*i)->getType() == SigmaZero) 569 G4double decayTime = SigmaZeroDecay 570 G4double time = currentTime + decay 571 if(time <= maximumTime) { 572 registerAvatar(new DecayAvatar((* 573 } 574 } 575 if((*i)->isOmega()) { 576 G4double decayTimeOmega = PionResona 577 G4double timeOmega = currentTime + d 578 if(timeOmega <= maximumTime) { 579 registerAvatar(new DecayAvatar((*i 580 } 581 } 582 } 583 } 584 585 G4INCL::IAvatar* StandardPropagationModel: 586 { 587 if(fs) { 588 // We update only the information rela 589 // by the previous avatar. 590 #ifdef INCL_REGENERATE_AVATARS 591 #warning "The INCL_REGENERATE_AVATARS code has 592 if(!fs->getModifiedParticles().empty() 593 // Regenerates the entire avatar lis 594 // updated particles 595 theNucleus->getStore()->clearAvatars 596 theNucleus->getStore()->initialisePa 597 generateAllAvatarsExceptUpdated(fs); 598 } 599 #else 600 ParticleList const &updatedParticles = 601 if(fs->getValidity()==PauliBlockedFS) 602 // This final state might represents 603 // decay 604 // assert(updatedParticles.empty() || (updated 605 // assert(fs->getEnteringParticles().empty() & 606 generateDecays(updatedParticles); 607 } else { 608 ParticleList const &entering = fs->g 609 generateDecays(updatedParticles); 610 generateDecays(entering); 611 612 ParticleList const &created = fs->ge 613 if(created.empty() && entering.empty 614 updateAvatars(updatedParticles); 615 else { 616 ParticleList updatedParticlesCopy 617 updatedParticlesCopy.insert(update 618 updatedParticlesCopy.insert(update 619 updateAvatars(updatedParticlesCopy 620 } 621 } 622 #endif 623 } 624 625 G4INCL::IAvatar *theAvatar = theNucleus- 626 if(theAvatar == 0) return 0; // Avatar l 627 // theAvatar->dispose(); 628 629 if(theAvatar->getTime() < currentTime) { 630 INCL_ERROR("Avatar time = " << theAvat 631 return 0; 632 } else if(theAvatar->getTime() > current 633 theNucleus->getStore()->timeStep(theAv 634 635 currentTime = theAvatar->getTime(); 636 theNucleus->getStore()->getBook().setC 637 } 638 639 return theAvatar; 640 } 641 } 642