34 #define INCLXX_IN_GEANT4_MODE 1
67 particle1(p1), particle2(NULL),
70 violationEFunctor(NULL)
77 particle1(p1), particle2(p2),
78 isPiN((p1->isPion() && p2->isNucleon()) || (p2->isPion() && p1->isNucleon())),
80 violationEFunctor(NULL)
97 (*backupParticle1) = (*particle1);
103 (*backupParticle2) = (*particle2);
146 const short maxIterations=50;
148 if(pos2 < r*r)
return true;
150 while( pos2 >= r*r && iterations<maxIterations )
152 pos *= std::sqrt(r*r*0.9801/pos2);
156 if( iterations < maxIterations)
167 INCL_DEBUG(
"postInteraction: final state: " <<
'\n' << fs->
print() <<
'\n');
185 if(((*i)->isPion() || (*i)->isKaon() || (*i)->isAntiKaon()) && (*i)->getPosition().mag() >
theNucleus->
getSurfaceRadius(*i)) {
186 (*i)->makeParticipant();
187 (*i)->setOutOfWell();
189 INCL_DEBUG(
"Pion was created outside its potential well." <<
'\n'
197 INCL_DEBUG(
"Enforcing energy conservation: failed!" <<
'\n');
212 INCL_DEBUG(
"Enforcing energy conservation: success!" <<
'\n');
214 INCL_DEBUG(
"postInteraction after energy conservation: final state: " <<
'\n' << fs->
print() <<
'\n');
217 for(
ParticleIter i=modified.begin(),
e=modified.end(); i!=
e; ++i )
218 if((*i)->isDelta() &&
220 INCL_DEBUG(
"Mass of the produced delta below decay threshold; forbidding collision. deltaMass=" <<
221 (*i)->getMass() <<
'\n');
284 if((*i)->isOutOfWell())
continue;
287 if( !successBringParticlesInside ) {
288 INCL_ERROR(
"Failed to bring particle inside the nucleus!" <<
'\n');
294 std::vector<G4int> newBiasCollisionVector;
296 if(std::fabs(
weight-1.) > 1
E-6){
302 (*i)->setBiasCollisionVector(newBiasCollisionVector);
303 if(!(*i)->isOutOfWell()) {
306 G4bool goesBackToSpectator =
false;
308 G4double threshold = (*i)->getPotentialEnergy();
309 if((*i)->getType()==
Proton)
311 if((*i)->getKineticEnergy() < threshold)
312 goesBackToSpectator =
true;
315 (*i)->thawPropagation();
318 if(goesBackToSpectator) {
319 INCL_DEBUG(
"The following particle goes back to spectator:" <<
'\n'
320 << (*i)->print() <<
'\n');
321 if(!(*i)->isTargetSpectator()) {
324 (*i)->makeTargetSpectator();
326 if((*i)->isTargetSpectator()) {
329 (*i)->makeParticipant();
334 for(
ParticleIter i=destroyed.begin(),
e=destroyed.end(); i!=
e; ++i )
335 if(!(*i)->isTargetSpectator())
341 (*particle1) = (*backupParticle1);
343 (*particle2) = (*backupParticle2);
363 if(manyBodyFinalState)
377 (*violationEFunctor)(theSolution.
x);
379 INCL_DEBUG(
"Couldn't enforce energy conservation after an interaction, root-finding algorithm failed." <<
'\n');
393 finalParticles(modAndCre),
394 initialEnergy(totalEnergyBeforeInteraction),
397 shouldUseLocalEnergy(localE)
408 particleMomenta.clear();
412 scaleParticleMomenta(alpha);
415 for(
ParticleIter i=finalParticles.begin(),
e=finalParticles.end(); i!=
e; ++i)
416 deltaE += (*i)->getEnergy() - (*i)->getPotentialEnergy();
417 deltaE -= initialEnergy;
423 std::vector<ThreeVector>::const_iterator iP = particleMomenta.begin();
424 for(
ParticleIter i=finalParticles.begin(),
e=finalParticles.end(); i!=
e; ++i, ++iP) {
425 (*i)->setMomentum((*iP)*alpha);
426 (*i)->adjustEnergyFromMomentum();
432 (*i)->setPotentialEnergy(0.);
436 !(*i)->isKaon() && !(*i)->isAntiKaon() && !(*i)->isSigma() && !(*i)->isLambda()) {
442 for(
G4int iterLocE=0;
446 (*i)->setEnergy(energy + locE);
447 (*i)->adjustMomentumFromEnergy();
461 for(
G4int iterLocE=0;
465 (*i)->setEnergy(energy + locE);
466 (*i)->adjustMomentumFromEnergy();
477 scaleParticleMomenta(1.);
486 initialEnergy(totalEnergyBeforeInteraction),
488 theParticle(aParticle),
489 theEnergy(theParticle->getEnergy()),
490 theMomentum(theParticle->getMomentum()),
498 setParticleEnergy(alpha);
499 return theParticle->getEnergy() - theParticle->getPotentialEnergy() - initialEnergy;
512 for(
G4int iterLocE=0;
516 G4double particleEnergy = energyThreshold + locE + alpha*(theEnergy-energyThreshold);
517 const G4double theMass2 = std::pow(particleEnergy,2.)-theMomentum.mag2();
520 theMass = std::sqrt(theMass2);
523 particleEnergy = energyThreshold;
525 theParticle->setMass(theMass);
526 theParticle->setEnergy(particleEnergy);
542 setParticleEnergy(1.);