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 #ifndef G4INCLCascade_hh
00038 #define G4INCLCascade_hh 1
00039
00040 #include "G4INCLParticle.hh"
00041 #include "G4INCLNucleus.hh"
00042 #include "G4INCLIPropagationModel.hh"
00043 #include "G4INCLEventAction.hh"
00044 #include "G4INCLPropagationAction.hh"
00045 #include "G4INCLAvatarAction.hh"
00046 #include "G4INCLEventInfo.hh"
00047 #include "G4INCLGlobalInfo.hh"
00048 #include "G4INCLLogger.hh"
00049 #include "G4INCLConfig.hh"
00050 #include "G4INCLRootFinder.hh"
00051
00052 namespace G4INCL {
00053 class INCL {
00054 public:
00055 INCL(Config const * const config);
00056
00057 ~INCL();
00058
00059 G4bool prepareReaction(const ParticleSpecies &projectileSpecies, const G4double kineticEnergy, const G4int A, const G4int Z);
00060 G4bool initializeTarget(const G4int A, const G4int Z);
00061 inline const EventInfo &processEvent() {
00062 return processEvent(
00063 theConfig->getProjectileSpecies(),
00064 theConfig->getProjectileKineticEnergy(),
00065 theConfig->getTargetA(),
00066 theConfig->getTargetZ()
00067 );
00068 }
00069 const EventInfo &processEvent(
00070 ParticleSpecies const &projectileSpecies,
00071 const G4double kineticEnergy,
00072 const G4int targetA,
00073 const G4int targetZ
00074 );
00075
00076 void finalizeGlobalInfo();
00077 const GlobalInfo &getGlobalInfo() const { return theGlobalInfo; }
00078
00079 std::string configToString() { return theConfig->echo(); }
00080
00081 private:
00082 IPropagationModel *propagationModel;
00083 G4int theA, theZ;
00084 G4bool targetInitSuccess;
00085 G4double maxImpactParameter;
00086 G4double maxUniverseRadius;
00087 G4double maxInteractionDistance;
00088 G4double fixedImpactParameter;
00089 EventAction *eventAction;
00090 PropagationAction *propagationAction;
00091 AvatarAction *avatarAction;
00092 Config const * const theConfig;
00093 Nucleus *nucleus;
00094
00095 EventInfo theEventInfo;
00096 GlobalInfo theGlobalInfo;
00097
00099 G4int minRemnantSize;
00100
00102 class RecoilFunctor : public RootFunctor {
00103 public:
00108 RecoilFunctor(Nucleus * const n, const EventInfo &ei) :
00109 RootFunctor(0., 1E6),
00110 nucleus(n),
00111 outgoingParticles(n->getStore()->getOutgoingParticles()),
00112 theEventInfo(ei) {
00113 for(ParticleIter p=outgoingParticles.begin(); p!=outgoingParticles.end(); ++p) {
00114 particleMomenta.push_back((*p)->getMomentum());
00115 particleKineticEnergies.push_back((*p)->getKineticEnergy());
00116 }
00117 }
00118 virtual ~RecoilFunctor() {}
00119
00125 G4double operator()(const G4double x) const {
00126 scaleParticleEnergies(x);
00127 return nucleus->getConservationBalance(theEventInfo,true).energy;
00128 }
00129
00131 void cleanUp(const G4bool success) const {
00132 if(!success)
00133 scaleParticleEnergies(1.);
00134 }
00135
00136 private:
00138 Nucleus *nucleus;
00140 ParticleList const &outgoingParticles;
00141
00142 EventInfo const &theEventInfo;
00144 std::list<ThreeVector> particleMomenta;
00146 std::list<G4double> particleKineticEnergies;
00147
00152 void scaleParticleEnergies(const G4double rescale) const {
00153
00154 ThreeVector pBalance = nucleus->getIncomingMomentum();
00155 std::list<ThreeVector>::const_iterator iP = particleMomenta.begin();
00156 std::list<G4double>::const_iterator iE = particleKineticEnergies.begin();
00157 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i, ++iP, ++iE)
00158 {
00159 const G4double mass = (*i)->getMass();
00160 const G4double newKineticEnergy = (*iE) * rescale;
00161
00162 (*i)->setMomentum(*iP);
00163 (*i)->setEnergy(mass + newKineticEnergy);
00164 (*i)->adjustMomentumFromEnergy();
00165
00166 pBalance -= (*i)->getMomentum();
00167 }
00168
00169 nucleus->setMomentum(pBalance);
00170 const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
00171 const G4double pRem2 = pBalance.mag2();
00172 const G4double recoilEnergy = pRem2/
00173 (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
00174 nucleus->setEnergy(remnantMass + recoilEnergy);
00175 }
00176 };
00177
00179 class RecoilCMFunctor : public RootFunctor {
00180 public:
00185 RecoilCMFunctor(Nucleus * const n, const EventInfo &ei) :
00186 RootFunctor(0., 1E6),
00187 nucleus(n),
00188 theIncomingMomentum(nucleus->getIncomingMomentum()),
00189 outgoingParticles(n->getStore()->getOutgoingParticles()),
00190 theEventInfo(ei) {
00191 thePTBoostVector = nucleus->getIncomingMomentum()/nucleus->getInitialEnergy();
00192 for(ParticleIter p=outgoingParticles.begin(); p!=outgoingParticles.end(); ++p) {
00193 (*p)->boost(thePTBoostVector);
00194 particleCMMomenta.push_back((*p)->getMomentum());
00195 }
00196 }
00197 virtual ~RecoilCMFunctor() {}
00198
00204 G4double operator()(const G4double x) const {
00205 scaleParticleCMMomenta(x);
00206 return nucleus->getConservationBalance(theEventInfo,true).energy;
00207 }
00208
00210 void cleanUp(const G4bool success) const {
00211 if(!success)
00212 scaleParticleCMMomenta(1.);
00213 }
00214
00215 private:
00217 Nucleus *nucleus;
00219 ThreeVector thePTBoostVector;
00221 ThreeVector theIncomingMomentum;
00223 ParticleList const &outgoingParticles;
00224
00225 EventInfo const &theEventInfo;
00227 std::list<ThreeVector> particleCMMomenta;
00228
00233 void scaleParticleCMMomenta(const G4double rescale) const {
00234
00235 ThreeVector remnantMomentum = theIncomingMomentum;
00236 std::list<ThreeVector>::const_iterator iP = particleCMMomenta.begin();
00237 for( ParticleIter i = outgoingParticles.begin(); i != outgoingParticles.end(); ++i, ++iP)
00238 {
00239 (*i)->setMomentum(*iP * rescale);
00240 (*i)->adjustEnergyFromMomentum();
00241 (*i)->boost(-thePTBoostVector);
00242
00243 remnantMomentum -= (*i)->getMomentum();
00244 }
00245
00246 nucleus->setMomentum(remnantMomentum);
00247 const G4double remnantMass = ParticleTable::getTableMass(nucleus->getA(),nucleus->getZ()) + nucleus->getExcitationEnergy();
00248 const G4double pRem2 = remnantMomentum.mag2();
00249 const G4double recoilEnergy = pRem2/
00250 (std::sqrt(pRem2+remnantMass*remnantMass) + remnantMass);
00251 nucleus->setEnergy(remnantMass + recoilEnergy);
00252 }
00253 };
00254
00260 void rescaleOutgoingForRecoil();
00261
00262 #ifndef INCLXX_IN_GEANT4_MODE
00263
00272 void globalConservationChecks(G4bool afterRecoil);
00273 #endif
00274
00280 G4bool continueCascade();
00281
00299 G4int makeProjectileRemnant();
00300
00308 void makeCompoundNucleus();
00309
00311 G4bool preCascade(ParticleSpecies const projectileSpecies, const G4double kineticEnergy);
00312
00314 void cascade();
00315
00317 void postCascade();
00318
00323 void initMaxInteractionDistance(ParticleSpecies const &p, const G4double kineticEnergy);
00324
00330 void initUniverseRadius(ParticleSpecies const &p, const G4double kineticEnergy, const G4int A, const G4int Z);
00331 };
00332 }
00333
00334 #endif