00001 // 00002 // ******************************************************************** 00003 // * License and Disclaimer * 00004 // * * 00005 // * The Geant4 software is copyright of the Copyright Holders of * 00006 // * the Geant4 Collaboration. It is provided under the terms and * 00007 // * conditions of the Geant4 Software License, included in the file * 00008 // * LICENSE and available at http://cern.ch/geant4/license . These * 00009 // * include a list of copyright holders. * 00010 // * * 00011 // * Neither the authors of this software system, nor their employing * 00012 // * institutes,nor the agencies providing financial support for this * 00013 // * work make any representation or warranty, express or implied, * 00014 // * regarding this software system or assume any liability for its * 00015 // * use. Please see the license in the file LICENSE and URL above * 00016 // * for the full disclaimer and the limitation of liability. * 00017 // * * 00018 // * This code implementation is the result of the scientific and * 00019 // * technical work of the GEANT4 collaboration. * 00020 // * By using, copying, modifying or distributing the software (or * 00021 // * any work based on the software) you agree to acknowledge its * 00022 // * use in resulting scientific publications, and indicate your * 00023 // * acceptance of all terms of the Geant4 Software license. * 00024 // ******************************************************************** 00025 // 00026 // INCL++ intra-nuclear cascade model 00027 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics 00028 // Davide Mancusi, CEA 00029 // Alain Boudard, CEA 00030 // Sylvie Leray, CEA 00031 // Joseph Cugnon, University of Liege 00032 // 00033 #define INCLXX_IN_GEANT4_MODE 1 00034 00035 #include "globals.hh" 00036 00037 /* 00038 * G4INCLNucleus.hh 00039 * 00040 * \date Jun 5, 2009 00041 * \author Pekka Kaitaniemi 00042 */ 00043 00044 #ifndef G4INCLNUCLEUS_HH_ 00045 #define G4INCLNUCLEUS_HH_ 00046 00047 #include <list> 00048 #include <string> 00049 00050 #include "G4INCLParticle.hh" 00051 #include "G4INCLEventInfo.hh" 00052 #include "G4INCLCluster.hh" 00053 #include "G4INCLFinalState.hh" 00054 #include "G4INCLStore.hh" 00055 #include "G4INCLGlobals.hh" 00056 #include "G4INCLParticleTable.hh" 00057 #include "G4INCLConfig.hh" 00058 #include "G4INCLConfigEnums.hh" 00059 #include "G4INCLCluster.hh" 00060 #include "G4INCLProjectileRemnant.hh" 00061 00062 namespace G4INCL { 00063 00064 class Nucleus : public Cluster { 00065 public: 00066 Nucleus(G4int mass, G4int charge, Config const * const conf, const G4double universeRadius=-1.); 00067 virtual ~Nucleus(); 00068 00073 void initializeParticles(); 00074 00076 void insertParticle(Particle *p) { 00077 theZ += p->getZ(); 00078 theA += p->getA(); 00079 theStore->particleHasEntered(p); 00080 if(p->isNucleon()) { 00081 theNpInitial += Math::heaviside(ParticleTable::getIsospin(p->getType())); 00082 theNnInitial += Math::heaviside(-ParticleTable::getIsospin(p->getType())); 00083 } 00084 if(!p->isTargetSpectator()) theStore->getBook()->incrementCascading(); 00085 }; 00086 00090 void applyFinalState(FinalState *); 00091 00092 G4int getInitialA() const { return theInitialA; }; 00093 G4int getInitialZ() const { return theInitialZ; }; 00094 00098 ParticleList const &getCreatedParticles() const { return justCreated; } 00099 00103 ParticleList const &getUpdatedParticles() const { return toBeUpdated; } 00104 00106 Particle *getBlockedDelta() const { return blockedDelta; } 00107 00113 void propagateParticles(G4double step); 00114 00115 G4int getNumberOfEnteringProtons() const { return theNpInitial; }; 00116 G4int getNumberOfEnteringNeutrons() const { return theNnInitial; }; 00117 00122 G4double computeSeparationEnergyBalance() const { 00123 G4double S = 0.0; 00124 ParticleList outgoing = theStore->getOutgoingParticles(); 00125 for(ParticleIter i = outgoing.begin(); i != outgoing.end(); ++i) { 00126 const ParticleType t = (*i)->getType(); 00127 switch(t) { 00128 case Proton: 00129 case Neutron: 00130 case DeltaPlusPlus: 00131 case DeltaPlus: 00132 case DeltaZero: 00133 case DeltaMinus: 00134 S += thePotential->getSeparationEnergy(*i); 00135 break; 00136 case Composite: 00137 S += (*i)->getZ() * thePotential->getSeparationEnergy(Proton) 00138 + ((*i)->getA() - (*i)->getZ()) * thePotential->getSeparationEnergy(Neutron); 00139 break; 00140 case PiPlus: 00141 S += thePotential->getSeparationEnergy(Proton) - thePotential->getSeparationEnergy(Neutron); 00142 break; 00143 case PiMinus: 00144 S += thePotential->getSeparationEnergy(Neutron) - thePotential->getSeparationEnergy(Proton); 00145 break; 00146 default: 00147 break; 00148 } 00149 } 00150 00151 S -= theNpInitial * thePotential->getSeparationEnergy(Proton); 00152 S -= theNnInitial * thePotential->getSeparationEnergy(Neutron); 00153 return S; 00154 } 00155 00160 G4bool decayOutgoingDeltas(); 00161 00166 G4bool decayInsideDeltas(); 00167 00172 G4bool decayOutgoingClusters(); 00173 00180 G4bool decayMe(); 00181 00183 void emitInsidePions(); 00184 00186 void computeRecoilKinematics(); 00187 00192 ThreeVector computeCenterOfMass() const; 00193 00198 G4double computeTotalEnergy() const; 00199 00204 G4double computeExcitationEnergy() const; 00205 00207 void setIncomingAngularMomentum(const ThreeVector &j) { 00208 incomingAngularMomentum = j; 00209 } 00210 00212 const ThreeVector &getIncomingAngularMomentum() const { return incomingAngularMomentum; } 00213 00215 void setIncomingMomentum(const ThreeVector &p) { 00216 incomingMomentum = p; 00217 } 00218 00220 const ThreeVector &getIncomingMomentum() const { 00221 return incomingMomentum; 00222 } 00223 00225 void setInitialEnergy(const G4double e) { initialEnergy = e; } 00226 00228 G4double getInitialEnergy() const { return initialEnergy; } 00229 00234 G4double getExcitationEnergy() const { return theExcitationEnergy; } 00235 00237 inline G4bool containsDeltas() { 00238 ParticleList inside = theStore->getParticles(); 00239 for(ParticleIter i=inside.begin(); i!=inside.end(); ++i) 00240 if((*i)->isDelta()) return true; 00241 return false; 00242 } 00243 00247 std::string print(); 00248 00249 std::string dump(); 00250 00251 Store* getStore() const {return theStore; }; 00252 void setStore(Store *s) { 00253 delete theStore; 00254 theStore = s; 00255 }; 00256 00257 G4double getInitialInternalEnergy() const { return initialInternalEnergy; }; 00258 00263 G4bool isEventTransparent() const; 00264 00269 G4bool hasRemnant() const { return remnant; } 00270 00274 // void fillEventInfo(Results::EventInfo *eventInfo); 00275 void fillEventInfo(EventInfo *eventInfo); 00276 00277 G4bool getTryCompoundNucleus() { return tryCN; } 00278 00280 G4int getProjectileChargeNumber() const { return projectileZ; } 00281 00283 G4int getProjectileMassNumber() const { return projectileA; } 00284 00286 void setProjectileChargeNumber(G4int n) { projectileZ=n; } 00287 00289 void setProjectileMassNumber(G4int n) { projectileA=n; } 00290 00292 G4bool isForcedTransparent() { return forceTransparent; } 00293 00295 G4double getTransmissionBarrier(Particle const * const p) { 00296 const G4double theTransmissionRadius = theDensity->getTransmissionRadius(p); 00297 const G4double theParticleZ = p->getZ(); 00298 return PhysicalConstants::eSquared*(theZ-theParticleZ)*theParticleZ/theTransmissionRadius; 00299 } 00300 00302 struct ConservationBalance { 00303 ThreeVector momentum; 00304 G4double energy; 00305 G4int Z, A; 00306 }; 00307 00309 ConservationBalance getConservationBalance(EventInfo const &theEventInfo, const G4bool afterRecoil) const; 00310 00312 void useFusionKinematics(); 00313 00322 G4double getSurfaceRadius(Particle const * const particle) const { 00323 if(particle->isPion()) 00324 // Temporarily set RPION = RMAX 00325 return getUniverseRadius(); 00326 //return 0.5*(theDensity->getTransmissionRadius(particle)+getUniverseRadius()); 00327 else { 00328 const G4double pr = particle->getMomentum().mag()/thePotential->getFermiMomentum(particle); 00329 if(pr>=1.) 00330 return getUniverseRadius(); 00331 else 00332 return theDensity->getMaxRFromP(pr); 00333 } 00334 } 00335 00350 G4double getCoulombRadius(ParticleSpecies const &p) const { 00351 if(p.theType == Composite) { 00352 const G4int zp = p.theZ; 00353 const G4int ap = p.theA; 00354 G4double barr, radius = 0.; 00355 if(zp==1 && ap==2) { // d 00356 barr = 0.2565*Math::pow23((G4double)theA)-0.78; 00357 radius = ParticleTable::eSquared*zp*theZ/barr - 2.5; 00358 } else if((zp==1 || zp==2) && ap==3) { // t, He3 00359 barr = 0.5*(0.5009*Math::pow23((G4double)theA)-1.16); 00360 radius = ParticleTable::eSquared*theZ/barr - 0.5; 00361 } else if(zp==2 && ap==4) { // alpha 00362 barr = 0.5939*Math::pow23((G4double)theA)-1.64; 00363 radius = ParticleTable::eSquared*zp*theZ/barr - 0.5; 00364 } else if(zp>2) { 00365 radius = getUniverseRadius(); 00366 } 00367 if(radius<=0.) { 00368 ERROR("Negative Coulomb radius! Using the universe radius" << std::endl); 00369 radius = getUniverseRadius(); 00370 } 00371 return radius; 00372 } else 00373 return getUniverseRadius(); 00374 } 00375 00377 G4double getUniverseRadius() const { return theUniverseRadius; } 00378 00380 void setUniverseRadius(const G4double universeRadius) { theUniverseRadius=universeRadius; } 00381 00383 G4bool isNucleusNucleusCollision() const { return isNucleusNucleus; } 00384 00386 void setNucleusNucleusCollision() { isNucleusNucleus=true; } 00387 00389 void setParticleNucleusCollision() { isNucleusNucleus=false; } 00390 00392 void setProjectileRemnant(ProjectileRemnant * const c) { 00393 delete theProjectileRemnant; 00394 theProjectileRemnant = c; 00395 } 00396 00398 ProjectileRemnant *getProjectileRemnant() const { return theProjectileRemnant; } 00399 00401 void deleteProjectileRemnant() { 00402 delete theProjectileRemnant; 00403 theProjectileRemnant = NULL; 00404 } 00405 00407 void moveProjectileRemnantComponentsToOutgoing() { 00408 theStore->addToOutgoing(theProjectileRemnant->getParticles()); 00409 theProjectileRemnant->clearParticles(); 00410 } 00411 00420 void finalizeProjectileRemnant(const G4double emissionTime); 00421 00423 inline void updatePotentialEnergy(Particle *p) { 00424 p->setPotentialEnergy(thePotential->computePotentialEnergy(p)); 00425 } 00426 00428 NuclearDensity* getDensity() const { return theDensity; }; 00429 00431 NuclearPotential::INuclearPotential* getPotential() const { return thePotential; }; 00432 00433 private: 00439 void computeOneNucleonRecoilKinematics(); 00440 00441 private: 00442 G4int theInitialZ, theInitialA; 00444 G4int theNpInitial; 00446 G4int theNnInitial; 00447 G4double initialInternalEnergy; 00448 ThreeVector incomingAngularMomentum, incomingMomentum; 00449 ThreeVector initialCenterOfMass; 00450 G4bool remnant; 00451 00452 ParticleList toBeUpdated; 00453 ParticleList justCreated; 00454 Particle *blockedDelta; 00455 G4double initialEnergy; 00456 Store *theStore; 00457 G4bool tryCN; 00458 G4bool forceTransparent; 00459 00461 G4int projectileZ; 00463 G4int projectileA; 00464 00466 G4double theUniverseRadius; 00467 00472 G4bool isNucleusNucleus; 00473 00475 ProjectileRemnant *theProjectileRemnant; 00476 00478 NuclearDensity *theDensity; 00479 00481 NuclearPotential::INuclearPotential *thePotential; 00482 00483 }; 00484 00485 } 00486 00487 #endif /* G4INCLNUCLEUS_HH_ */