00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00037
00038
00039
00040
00041
00042
00043
00044 #include "G4INCLStandardPropagationModel.hh"
00045 #include "G4INCLSurfaceAvatar.hh"
00046 #include "G4INCLBinaryCollisionAvatar.hh"
00047 #include "G4INCLDecayAvatar.hh"
00048 #include "G4INCLCrossSections.hh"
00049 #include "G4INCLRandom.hh"
00050 #include <iostream>
00051 #include <algorithm>
00052 #include "G4INCLLogger.hh"
00053 #include "G4INCLGlobals.hh"
00054 #include "G4INCLKinematicsUtils.hh"
00055 #include "G4INCLCoulombDistortion.hh"
00056 #include "G4INCLDeltaDecayChannel.hh"
00057 #include "G4INCLParticleEntryAvatar.hh"
00058 #include "G4INCLIntersection.hh"
00059
00060 namespace G4INCL {
00061
00062 StandardPropagationModel::StandardPropagationModel(LocalEnergyType localEnergyType, LocalEnergyType localEnergyDeltaType)
00063 :theNucleus(0), maximumTime(70.0), currentTime(0.0), firstAvatar(true),
00064 theLocalEnergyType(localEnergyType),
00065 theLocalEnergyDeltaType(localEnergyDeltaType)
00066 {
00067 }
00068
00069 StandardPropagationModel::~StandardPropagationModel()
00070 {
00071 delete theNucleus;
00072 }
00073
00074 G4INCL::Nucleus* StandardPropagationModel::getNucleus()
00075 {
00076 return theNucleus;
00077 }
00078
00079 G4double StandardPropagationModel::shoot(ParticleSpecies const projectileSpecies, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
00080 if(projectileSpecies.theType==Composite)
00081 return shootComposite(projectileSpecies, kineticEnergy, impactParameter, phi);
00082 else
00083 return shootParticle(projectileSpecies.theType, kineticEnergy, impactParameter, phi);
00084 }
00085
00086 G4double StandardPropagationModel::shootParticle(ParticleType const type, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
00087 theNucleus->setParticleNucleusCollision();
00088 currentTime = 0.0;
00089
00090
00091 const G4double projectileMass = ParticleTable::getTableParticleMass(type);
00092 G4double energy = kineticEnergy + projectileMass;
00093 G4double momentumZ = std::sqrt(energy*energy - projectileMass*projectileMass);
00094 ThreeVector momentum(0.0, 0.0, momentumZ);
00095 Particle *p= new G4INCL::Particle(type, energy, momentum, ThreeVector());
00096
00097 G4double temfin = 0.0;
00098 if( p->isNucleon() )
00099 temfin = 29.8 * std::pow(theNucleus->getA(), 0.16);
00100 else {
00101 const G4double tlab = p->getEnergy() - p->getMass();
00102 temfin = 30.18 * std::pow(theNucleus->getA(), 0.17*(1.0 - 5.7E-5*tlab));
00103 }
00104
00105 maximumTime = temfin;
00106
00107
00108 const G4double rMax = theNucleus->getUniverseRadius();
00109 const G4double distance = 2.*rMax;
00110 const G4double projectileVelocity = p->boostVector().mag();
00111 const G4double traversalTime = distance / projectileVelocity;
00112 if(maximumTime < traversalTime)
00113 maximumTime = traversalTime;
00114 DEBUG("Cascade stopping time is " << maximumTime << std::endl);
00115
00116
00117
00118
00119 if(impactParameter>CoulombDistortion::maxImpactParameter(p->getSpecies(), kineticEnergy, theNucleus)) {
00120 DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << std::endl);
00121 delete p;
00122 return -1.;
00123 }
00124
00125 ThreeVector position(impactParameter * std::cos(phi),
00126 impactParameter * std::sin(phi),
00127 0.);
00128 p->setPosition(position);
00129
00130
00131 theNucleus->setIncomingAngularMomentum(p->getAngularMomentum());
00132 theNucleus->setIncomingMomentum(p->getMomentum());
00133 theNucleus->setInitialEnergy(p->getEnergy()
00134 + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ()));
00135
00136
00137 p->setINCLMass();
00138 p->setEnergy(p->getMass() + kineticEnergy);
00139 p->adjustMomentumFromEnergy();
00140
00141 p->makeProjectileSpectator();
00142 generateAllAvatars();
00143 firstAvatar = false;
00144
00145
00146 ParticleEntryAvatar *theEntryAvatar = CoulombDistortion::bringToSurface(p, theNucleus);
00147 if(theEntryAvatar) {
00148 theNucleus->getStore()->addParticleEntryAvatar(theEntryAvatar);
00149
00150 theNucleus->setProjectileChargeNumber(p->getZ());
00151 theNucleus->setProjectileMassNumber(p->getA());
00152
00153 return p->getTransversePosition().mag();
00154 } else {
00155 delete p;
00156 return -1.;
00157 }
00158 }
00159
00160 G4double StandardPropagationModel::shootComposite(ParticleSpecies const species, const G4double kineticEnergy, const G4double impactParameter, const G4double phi) {
00161 theNucleus->setNucleusNucleusCollision();
00162 currentTime = 0.0;
00163
00164
00165 ProjectileRemnant *pr = new ProjectileRemnant(species, kineticEnergy);
00166
00167
00168 maximumTime = 29.8 * std::pow(theNucleus->getA(), 0.16);
00169
00170
00171 const G4double rms = ParticleTable::getNuclearRadius(pr->getA(), pr->getZ());
00172 const G4double rMax = theNucleus->getUniverseRadius();
00173 const G4double distance = 2.*rMax + 2.725*rms;
00174 const G4double projectileVelocity = pr->boostVector().mag();
00175 const G4double traversalTime = distance / projectileVelocity;
00176 if(maximumTime < traversalTime)
00177 maximumTime = traversalTime;
00178 DEBUG("Cascade stopping time is " << maximumTime << std::endl);
00179
00180
00181
00182
00183 if(impactParameter>CoulombDistortion::maxImpactParameter(pr,theNucleus)) {
00184 pr->deleteParticles();
00185 DEBUG("impactParameter>CoulombDistortion::maxImpactParameter" << std::endl);
00186 delete pr;
00187 return -1.;
00188 }
00189
00190
00191 ThreeVector position(impactParameter * std::cos(phi),
00192 impactParameter * std::sin(phi),
00193 0.);
00194 pr->setPosition(position);
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 pr->storeComponents();
00206
00207
00208 theNucleus->setIncomingAngularMomentum(pr->getAngularMomentum());
00209 theNucleus->setIncomingMomentum(pr->getMomentum());
00210 theNucleus->setInitialEnergy(pr->getEnergy()
00211 + ParticleTable::getTableMass(theNucleus->getA(),theNucleus->getZ()));
00212
00213 generateAllAvatars();
00214 firstAvatar = false;
00215
00216
00217 IAvatarList theAvatarList
00218 = CoulombDistortion::bringToSurface(pr, theNucleus);
00219
00220 if(theAvatarList.empty()) {
00221 DEBUG("No ParticleEntryAvatar found, transparent event" << std::endl);
00222 pr->deleteParticles();
00223 delete pr;
00224 return -1.;
00225 }
00226
00227
00228 theNucleus->setProjectileRemnant(pr);
00229
00230
00231 theNucleus->setProjectileChargeNumber(pr->getZ());
00232 theNucleus->setProjectileMassNumber(pr->getA());
00233
00234
00235 theNucleus->getStore()->addParticleEntryAvatars(theAvatarList);
00236
00237 return pr->getTransversePosition().mag();
00238 }
00239
00240 G4double StandardPropagationModel::getStoppingTime() {
00241 return maximumTime;
00242 }
00243
00244 void StandardPropagationModel::setStoppingTime(G4double time) {
00245
00246 maximumTime = time;
00247 }
00248
00249 G4double StandardPropagationModel::getCurrentTime() {
00250 return currentTime;
00251 }
00252
00253 void StandardPropagationModel::setNucleus(G4INCL::Nucleus *nucleus)
00254 {
00255 theNucleus = nucleus;
00256 }
00257
00258 void StandardPropagationModel::registerAvatar(G4INCL::IAvatar *anAvatar)
00259 {
00260 if(anAvatar) theNucleus->getStore()->add(anAvatar);
00261 }
00262
00263 IAvatar *StandardPropagationModel::generateBinaryCollisionAvatar(Particle * const p1, Particle * const p2) const {
00264
00265 if(!p1->isParticipant() && !p2->isParticipant() && p1->getParticipantType()==p2->getParticipantType()) return NULL;
00266
00267
00268 if((p1->isResonance() && p2->isPion()) || (p1->isPion() && p2->isResonance()))
00269 return NULL;
00270
00271
00272 G4double minDistOfApproachSquared = 0.0;
00273 G4double t = getTime(p1, p2, &minDistOfApproachSquared);
00274 if(t>maximumTime || t<currentTime) return NULL;
00275
00276
00277
00278 G4bool hasLocalEnergy;
00279 if(p1->isPion() || p2->isPion())
00280 hasLocalEnergy = ((theLocalEnergyDeltaType == FirstCollisionLocalEnergy &&
00281 theNucleus->getStore()->getBook()->getAcceptedCollisions()==0) ||
00282 theLocalEnergyDeltaType == AlwaysLocalEnergy);
00283 else
00284 hasLocalEnergy = ((theLocalEnergyType == FirstCollisionLocalEnergy &&
00285 theNucleus->getStore()->getBook()->getAcceptedCollisions()==0) ||
00286 theLocalEnergyType == AlwaysLocalEnergy);
00287 const G4bool p1HasLocalEnergy = (hasLocalEnergy && !p1->isPion());
00288 const G4bool p2HasLocalEnergy = (hasLocalEnergy && !p2->isPion());
00289
00290 Particle backupParticle1 = *p1;
00291 if(p1HasLocalEnergy) {
00292 p1->propagate(t - currentTime);
00293 if(p1->getPosition().mag() > theNucleus->getSurfaceRadius(p1)) {
00294 *p1 = backupParticle1;
00295 return NULL;
00296 }
00297 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p1);
00298 }
00299 Particle backupParticle2 = *p2;
00300 if(p2HasLocalEnergy) {
00301 p2->propagate(t - currentTime);
00302 if(p2->getPosition().mag() > theNucleus->getSurfaceRadius(p2)) {
00303 *p2 = backupParticle2;
00304 if(p1HasLocalEnergy) {
00305 *p1 = backupParticle1;
00306 }
00307 return NULL;
00308 }
00309 KinematicsUtils::transformToLocalEnergyFrame(theNucleus, p2);
00310 }
00311
00312
00313 const G4double totalCrossSection = CrossSections::total(p1, p2);
00314 const G4double squareTotalEnergyInCM = KinematicsUtils::squareTotalEnergyInCM(p1,p2);
00315
00316
00317 if(p1HasLocalEnergy) {
00318 *p1 = backupParticle1;
00319 }
00320 if(p2HasLocalEnergy) {
00321 *p2 = backupParticle2;
00322 }
00323
00324
00325 if(theNucleus->getStore()->getBook()->getAcceptedCollisions()>0
00326 && p1->isNucleon() && p2->isNucleon()
00327 && squareTotalEnergyInCM < BinaryCollisionAvatar::cutNNSquared) return NULL;
00328
00329
00330 if(Math::tenPi*minDistOfApproachSquared > totalCrossSection) return NULL;
00331
00332
00333
00334
00335
00336 return new G4INCL::BinaryCollisionAvatar(t, totalCrossSection, theNucleus, p1, p2);
00337 }
00338
00339 G4double StandardPropagationModel::getReflectionTime(G4INCL::Particle const * const aParticle) {
00340 Intersection theIntersection(
00341 IntersectionFactory::getLaterTrajectoryIntersection(
00342 aParticle->getPosition(),
00343 aParticle->getPropagationVelocity(),
00344 theNucleus->getSurfaceRadius(aParticle)));
00345 G4double time;
00346 if(theIntersection.exists) {
00347 time = currentTime + theIntersection.time;
00348 } else {
00349 ERROR("Imaginary reflection time for particle: " << std::endl
00350 << aParticle->print());
00351 time = 10000.0;
00352 }
00353 return time;
00354 }
00355
00356 G4double StandardPropagationModel::getTime(G4INCL::Particle const * const particleA,
00357 G4INCL::Particle const * const particleB, G4double *minDistOfApproach) const
00358 {
00359 G4double time;
00360 G4INCL::ThreeVector t13 = particleA->getPropagationVelocity();
00361 t13 -= particleB->getPropagationVelocity();
00362 G4INCL::ThreeVector distance = particleA->getPosition();
00363 distance -= particleB->getPosition();
00364 const G4double t7 = t13.dot(distance);
00365 const G4double dt = t13.mag2();
00366 if(dt <= 1.0e-10) {
00367 (*minDistOfApproach) = 100000.0;
00368 return currentTime + 100000.0;
00369 } else {
00370 time = -t7/dt;
00371 }
00372 (*minDistOfApproach) = distance.mag2() + time * t7;
00373 return currentTime + time;
00374 }
00375
00376 void StandardPropagationModel::generateUpdatedCollisions(const ParticleList &updatedParticles, const ParticleList &particles) {
00377
00378
00379 for(ParticleIter updated = updatedParticles.begin(); updated != updatedParticles.end(); ++updated)
00380 {
00381
00382 for(ParticleIter particle = particles.begin(); particle != particles.end(); ++particle)
00383 {
00384
00385
00386
00387
00388 if((*particle)->isInList(updatedParticles)) continue;
00389
00390 registerAvatar(generateBinaryCollisionAvatar(*particle,*updated));
00391 }
00392 }
00393 }
00394
00395 void StandardPropagationModel::generateCollisions(const ParticleList &particles, const ParticleList &except) {
00396
00397 G4bool haveExcept;
00398 haveExcept=(except.size()!=0);
00399
00400
00401 for(ParticleIter p1 = particles.begin(); p1 != particles.end(); ++p1)
00402 {
00403
00404 ParticleIter p2 = p1;
00405 for(++p2; p2 != particles.end(); ++p2)
00406 {
00407
00408 if(haveExcept && (*p1)->isInList(except) && (*p2)->isInList(except)) continue;
00409
00410 registerAvatar(generateBinaryCollisionAvatar(*p1,*p2));
00411 }
00412 }
00413
00414 }
00415
00416 void StandardPropagationModel::updateAvatars(const ParticleList &particles) {
00417
00418 for(ParticleIter iter = particles.begin(); iter != particles.end(); ++iter) {
00419 G4double time = this->getReflectionTime(*iter);
00420 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*iter, time, theNucleus));
00421 }
00422 ParticleList const &p = theNucleus->getStore()->getParticles();
00423 generateUpdatedCollisions(particles, p);
00424 }
00425
00426 void StandardPropagationModel::generateAllAvatars(G4bool excludeUpdated) {
00427 ParticleList particles = theNucleus->getStore()->getParticles();
00428 if(particles.empty()) { ERROR("No particles inside the nucleus!" << std::endl); }
00429 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
00430 G4double time = this->getReflectionTime(*i);
00431 if(time <= maximumTime) registerAvatar(new SurfaceAvatar(*i, time, theNucleus));
00432 }
00433 ParticleList except;
00434 if(excludeUpdated)
00435 except = theNucleus->getUpdatedParticles();
00436 generateCollisions(particles,except);
00437 generateDecays(particles);
00438 }
00439
00440 void StandardPropagationModel::generateDecays(const ParticleList &particles) {
00441 for(ParticleIter i = particles.begin(); i != particles.end(); ++i) {
00442 if((*i)->isDelta()) {
00443 G4double decayTime = DeltaDecayChannel::computeDecayTime((*i));
00444 G4double time = currentTime + decayTime;
00445 if(time <= maximumTime) {
00446 registerAvatar(new DecayAvatar((*i), time, theNucleus));
00447 }
00448 }
00449 }
00450 }
00451
00452 G4INCL::IAvatar* StandardPropagationModel::propagate()
00453 {
00454
00455
00456 #ifdef INCL_REGENERATE_AVATARS
00457 #warning "The INCL_REGENERATE_AVATARS code has not been tested in a while. Use it at your peril."
00458 if(theNucleus->getUpdatedParticles().size()!=0 || theNucleus->getCreatedParticles().size()!=0) {
00459
00460
00461 theNucleus->getStore()->clearAvatars();
00462 theNucleus->getStore()->initialiseParticleAvatarConnections();
00463 generateAllAvatars(true);
00464 }
00465 #else
00466
00467
00468 Particle * const blockedDelta = theNucleus->getBlockedDelta();
00469 ParticleList updatedParticles = theNucleus->getUpdatedParticles();
00470 if(blockedDelta)
00471 updatedParticles.push_back(blockedDelta);
00472 generateDecays(updatedParticles);
00473
00474 ParticleList needNewAvatars = theNucleus->getUpdatedParticles();
00475 ParticleList created = theNucleus->getCreatedParticles();
00476 needNewAvatars.splice(needNewAvatars.end(), created);
00477 updateAvatars(needNewAvatars);
00478 #endif
00479
00480 G4INCL::IAvatar *theAvatar = theNucleus->getStore()->findSmallestTime();
00481 if(theAvatar == 0) return 0;
00482
00483
00484 if(theAvatar->getTime() < currentTime) {
00485 ERROR("Avatar time = " << theAvatar->getTime() << ", currentTime = " << currentTime << std::endl);
00486 return 0;
00487 } else if(theAvatar->getTime() > currentTime) {
00488 theNucleus->getStore()->timeStep(theAvatar->getTime() - currentTime);
00489
00490 currentTime = theAvatar->getTime();
00491 theNucleus->getStore()->getBook()->setCurrentTime(currentTime);
00492 }
00493
00494 return theAvatar;
00495 }
00496
00497 void StandardPropagationModel::putSpectatorsOnShell(IAvatarList const &entryAvatars, ParticleList const &spectators) {
00498 G4double deltaE = 0.0;
00499 for(ParticleIter p=spectators.begin(); p!=spectators.end(); ++p) {
00500
00501 const G4double oldEnergy = (*p)->getEnergy();
00502 (*p)->setTableMass();
00503 (*p)->adjustEnergyFromMomentum();
00504 deltaE += (*p)->getEnergy() - oldEnergy;
00505 }
00506
00507 deltaE /= entryAvatars.size();
00508
00509 for(IAvatarIter a=entryAvatars.begin(); a!=entryAvatars.end(); ++a) {
00510
00511 Particle *p = (*a)->getParticles().front();
00512 ParticleType const t = p->getType();
00513
00514 const G4double energy = p->getEnergy() - deltaE
00515 - ParticleTable::getTableParticleMass(t) + ParticleTable::getINCLMass(t);
00516 p->setEnergy(energy);
00517 const G4double newMass = std::sqrt(energy*energy - p->getMomentum().mag2());
00518 p->setMass(newMass);
00519 }
00520 }
00521 }