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 #include "G4INCLParticleEntryChannel.hh"
00038 #include "G4INCLRootFinder.hh"
00039 #include "G4INCLIntersection.hh"
00040 #include <algorithm>
00041
00042 namespace G4INCL {
00043
00044 ParticleEntryChannel::ParticleEntryChannel(Nucleus *n, Particle *p)
00045 :theNucleus(n), theParticle(p)
00046 {}
00047
00048 ParticleEntryChannel::~ParticleEntryChannel()
00049 {}
00050
00051 FinalState* ParticleEntryChannel::getFinalState() {
00052
00053 G4bool isNN = theNucleus->isNucleusNucleusCollision();
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079 G4double theCorrection;
00080 if(isNN) {
00081
00082 ProjectileRemnant * const projectileRemnant = theNucleus->getProjectileRemnant();
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
00115 const G4double theProjectileExcitationEnergy =
00116 (projectileRemnant->getA()-theParticle->getA()>1) ?
00117 (projectileRemnant->computeExcitationEnergy(theParticle->getID())) :
00118 0.;
00119 const G4double theProjectileEffectiveMass =
00120 ParticleTable::getTableMass(projectileRemnant->getA() - theParticle->getA(), projectileRemnant->getZ() - theParticle->getZ())
00121 + theProjectileExcitationEnergy;
00122 const ThreeVector &theProjectileMomentum = projectileRemnant->getMomentum() - theParticle->getMomentum();
00123 const G4double theProjectileEnergy = std::sqrt(theProjectileMomentum.mag2() + theProjectileEffectiveMass*theProjectileEffectiveMass);
00124 const G4double theProjectileCorrection = theProjectileEnergy - (projectileRemnant->getEnergy() - theParticle->getEnergy());
00125 theCorrection = theParticle->getEmissionQValueCorrection(
00126 theNucleus->getA() + theParticle->getA(),
00127 theNucleus->getZ() + theParticle->getZ())
00128 + theParticle->getTableMass() - theParticle->getINCLMass()
00129 + theProjectileCorrection;
00130
00131 projectileRemnant->removeParticle(theParticle, theProjectileCorrection);
00132 } else {
00133 const G4int ACN = theNucleus->getA() + theParticle->getA();
00134 const G4int ZCN = theNucleus->getZ() + theParticle->getZ();
00135
00136 theCorrection = theParticle->getEmissionQValueCorrection(ACN,ZCN);
00137 DEBUG("The following Particle enters with correction " << theCorrection
00138 << theParticle->print() << std::endl);
00139 }
00140
00141 const G4double energyBefore = theParticle->getEnergy() - theCorrection;
00142 G4bool success = particleEnters(theCorrection);
00143 FinalState *fs = new FinalState();
00144 fs->addEnteringParticle(theParticle);
00145
00146 if(!success) {
00147 fs->makeParticleBelowZero();
00148 } else if(theParticle->isNucleon() &&
00149 theParticle->getKineticEnergy()<theNucleus->getPotential()->getFermiEnergy(theParticle)) {
00150
00151
00152 fs->makeParticleBelowFermi();
00153 }
00154
00155 fs->setTotalEnergyBeforeInteraction(energyBefore);
00156 return fs;
00157 }
00158
00159 G4bool ParticleEntryChannel::particleEnters(const G4double theQValueCorrection) {
00160
00161
00162
00163 theParticle->setINCLMass();
00164
00165
00166
00167
00168 class IncomingEFunctor : public RootFunctor {
00169 public:
00170 IncomingEFunctor(Particle * const p, Nucleus const * const n, const G4double correction) :
00171 RootFunctor(0., 1E6),
00172 theParticle(p),
00173 thePotential(n->getPotential()),
00174 theEnergy(theParticle->getEnergy()),
00175 theMass(theParticle->getMass()),
00176 theQValueCorrection(correction),
00177 theMomentum(theParticle->getMomentum())
00178 {}
00179 ~IncomingEFunctor() {}
00180 G4double operator()(const G4double v) const {
00181 G4double energyInside = std::max(theMass, theEnergy + v - theQValueCorrection);
00182 theParticle->setEnergy(energyInside);
00183 theParticle->setPotentialEnergy(v);
00184
00185 theParticle->setMomentum(theMomentum);
00186 theParticle->adjustMomentumFromEnergy();
00187 return v - thePotential->computePotentialEnergy(theParticle);
00188 }
00189 void cleanUp(const G4bool ) const {}
00190 private:
00191 Particle *theParticle;
00192 NuclearPotential::INuclearPotential const *thePotential;
00193 G4double theEnergy;
00194 G4double theMass;
00195 G4double theQValueCorrection;
00196 ThreeVector theMomentum;
00197 } theIncomingEFunctor(theParticle,theNucleus,theQValueCorrection);
00198
00199 G4double v = theNucleus->getPotential()->computePotentialEnergy(theParticle);
00200 if(theParticle->getKineticEnergy()+v-theQValueCorrection<0.) {
00201 DEBUG("Particle " << theParticle->getID() << " is trying to enter below 0" << std::endl);
00202 return false;
00203 }
00204 G4bool success = RootFinder::solve(&theIncomingEFunctor, v);
00205 if(success) {
00206 std::pair<G4double,G4double> theSolution = RootFinder::getSolution();
00207 theIncomingEFunctor(theSolution.first);
00208 } else {
00209 WARN("Couldn't compute the potential for incoming particle, root-finding algorithm failed." << std::endl);
00210 }
00211 return success;
00212 }
00213
00214 }
00215