G4INCLProjectileRemnant.cc

Go to the documentation of this file.
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 
00044 #include "G4INCLProjectileRemnant.hh"
00045 #include <algorithm>
00046 #include <numeric>
00047 
00048 namespace G4INCL {
00049 
00050   G4int shuffleComponentsHelper(G4int range) {
00051     return (G4int)(Random::shoot1()*range);
00052   }
00053 
00054   void ProjectileRemnant::reset() {
00055     deleteParticles();
00056     thePosition = ThreeVector();
00057     theMomentum = ThreeVector();
00058     theEnergy = 0.0;
00059     thePotentialEnergy = 0.0;
00060     theA = 0;
00061     theZ = 0;
00062     nCollisions = 0;
00063 
00064     for(std::map<long, Particle*>::const_iterator i=storedComponents.begin(); i!=storedComponents.end(); ++i) {
00065       Particle *p = new Particle(*(i->second));
00066       EnergyLevelMap::iterator energyIter = theInitialEnergyLevels.find(i->first);
00067 // assert(energyIter!=theInitialEnergyLevels.end());
00068       const G4double energyLevel = energyIter->second;
00069       theInitialEnergyLevels.erase(energyIter);
00070       theInitialEnergyLevels[p->getID()] = energyLevel;
00071       addParticle(p);
00072     }
00073     thePosition /= theA;
00074     setTableMass();
00075     DEBUG("ProjectileRemnant object was reset:" << std::endl << print());
00076   }
00077 
00078   void ProjectileRemnant::removeParticle(Particle * const p, const G4double theProjectileCorrection) {
00079 // assert(p->isNucleon());
00080 
00081     DEBUG("The following Particle is about to be removed from the ProjectileRemnant:"
00082         << std::endl << p->print()
00083         << "theProjectileCorrection=" << theProjectileCorrection << std::endl);
00084     // Update A, Z, momentum and energy of the projectile remnant
00085     theA -= p->getA();
00086     theZ -= p->getZ();
00087 
00088     ThreeVector const &oldMomentum = p->getMomentum();
00089     const G4double oldEnergy = p->getEnergy();
00090     Cluster::removeParticle(p);
00091 
00092 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
00093     ThreeVector theTotalMomentum;
00094     G4double theTotalEnergy = 0.;
00095     const G4double theThreshold = 0.1;
00096 #endif
00097 
00098     if(getA()>0) { // if there are any particles left
00099 // assert((unsigned int)getA()==particles.size());
00100 
00101       const G4double theProjectileCorrectionPerNucleon = theProjectileCorrection / particles.size();
00102 
00103       // Update the kinematics of the components
00104       for(ParticleIter i=particles.begin(); i!=particles.end(); ++i) {
00105         (*i)->setEnergy((*i)->getEnergy() + theProjectileCorrectionPerNucleon);
00106         (*i)->setMass((*i)->getInvariantMass());
00107 #if !defined(NDEBUG) && !defined(INCLXX_IN_GEANT4_MODE)
00108         theTotalMomentum += (*i)->getMomentum();
00109         theTotalEnergy += (*i)->getEnergy();
00110 #endif
00111       }
00112     }
00113 
00114     theMomentum -= oldMomentum;
00115     theEnergy -= oldEnergy - theProjectileCorrection;
00116 
00117 // assert(std::abs((theTotalMomentum-theMomentum).mag())<theThreshold);
00118 // assert(std::abs(theTotalEnergy-theEnergy)<theThreshold);
00119     DEBUG("After Particle removal, the ProjectileRemnant looks like this:"
00120         << std::endl << print());
00121   }
00122 
00123   ParticleList ProjectileRemnant::addDynamicalSpectators(ParticleList pL) {
00124     // Try as hard as possible to add back all the dynamical spectators.
00125     // Don't add spectators that lead to negative excitation energies, but
00126     // iterate over the spectators as many times as possible, until
00127     // absolutely sure that all of them were rejected.
00128     unsigned int accepted;
00129     do {
00130       accepted = 0;
00131       ParticleList toBeAdded = pL;
00132       for(ParticleIter p=toBeAdded.begin(); p!=toBeAdded.end(); ++p) {
00133         G4bool isAccepted = addDynamicalSpectator(*p);
00134         if(isAccepted) {
00135           pL.remove(*p);
00136           accepted++;
00137         }
00138       }
00139     } while(accepted > 0);
00140     return pL;
00141   }
00142 
00143   ParticleList ProjectileRemnant::addMostDynamicalSpectators(ParticleList pL) {
00144     // Try as hard as possible to add back all the dynamical spectators.
00145     // Don't add spectators that lead to negative excitation energies. Start by
00146     // adding all of them, and repeatedly remove the most troublesome one until
00147     // the excitation energy becomes non-negative.
00148 
00149     // Put all the spectators in the projectile
00150     ThreeVector theNewMomentum = theMomentum;
00151     G4double theNewEnergy = theEnergy;
00152     G4int theNewA = theA;
00153     G4int theNewZ = theZ;
00154     for(ParticleIter p=pL.begin(); p!=pL.end(); ++p) {
00155 // assert((*p)->isNucleon());
00156       // Add the initial (off-shell) momentum and energy to the projectile remnant
00157       theNewMomentum += getStoredMomentum(*p);
00158       theNewEnergy += (*p)->getEnergy();
00159       theNewA += (*p)->getA();
00160       theNewZ += (*p)->getZ();
00161     }
00162 
00163     // Check that the excitation energy of the new projectile remnant is non-negative
00164     const G4double theNewMass = ParticleTable::getTableMass(theNewA,theNewZ);
00165     const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
00166 
00167     G4bool positiveExcitationEnergy = false;
00168     if(theNewInvariantMassSquared>=0.) {
00169       const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
00170       positiveExcitationEnergy = (theNewInvariantMass-theNewMass>-1.e-5);
00171     }
00172 
00173     // Keep removing nucleons from the projectile remnant until we achieve a
00174     // non-negative excitation energy.
00175     ParticleList rejected;
00176     while(!positiveExcitationEnergy && !pL.empty()) {
00177       G4double maxExcitationEnergy = -1.E30;
00178       ParticleList::iterator best = pL.end();
00179       ThreeVector bestMomentum;
00180       G4double bestEnergy = -1.;
00181       G4int bestA = -1, bestZ = -1;
00182       for(ParticleList::iterator p=pL.begin(); p!=pL.end(); ++p) {
00183         // Subtract the initial (off-shell) momentum and energy from the new
00184         // projectile remnant
00185         const ThreeVector theNewerMomentum = theNewMomentum - getStoredMomentum(*p);
00186         const G4double theNewerEnergy = theNewEnergy - (*p)->getEnergy();
00187         const G4int theNewerA = theNewA - (*p)->getA();
00188         const G4int theNewerZ = theNewZ - (*p)->getZ();
00189 
00190         const G4double theNewerMass = ParticleTable::getTableMass(theNewerA,theNewerZ);
00191         const G4double theNewerInvariantMassSquared = theNewerEnergy*theNewerEnergy-theNewerMomentum.mag2();
00192 
00193         if(theNewerInvariantMassSquared>=-1.e-5) {
00194           const G4double theNewerInvariantMass = std::sqrt(std::max(0.,theNewerInvariantMassSquared));
00195           const G4double theNewerExcitationEnergy = theNewerInvariantMass-theNewerMass;
00196           // Pick the nucleon that maximises the excitation energy of the
00197           // ProjectileRemnant
00198           if(theNewerExcitationEnergy>maxExcitationEnergy) {
00199             best = p;
00200             maxExcitationEnergy = theNewerExcitationEnergy;
00201             bestMomentum = theNewerMomentum;
00202             bestEnergy = theNewerEnergy;
00203             bestA = theNewerA;
00204             bestZ = theNewerZ;
00205           }
00206         }
00207       }
00208 
00209       // If we couldn't even calculate the excitation energy, fail miserably
00210       if(best==pL.end())
00211         return pL;
00212 
00213       rejected.push_back(*best);
00214       pL.erase(best);
00215       theNewMomentum = bestMomentum;
00216       theNewEnergy = bestEnergy;
00217       theNewA = bestA;
00218       theNewZ = bestZ;
00219 
00220       if(maxExcitationEnergy>0.) {
00221         // Stop here
00222         positiveExcitationEnergy = true;
00223       }
00224     }
00225 
00226     // Add the accepted participants to the projectile remnant
00227     for(ParticleIter p=pL.begin(); p!=pL.end(); ++p) {
00228       particles.push_back(*p);
00229     }
00230     theA = theNewA;
00231     theZ = theNewZ;
00232     theMomentum = theNewMomentum;
00233     theEnergy = theNewEnergy;
00234 
00235     return rejected;
00236   }
00237 
00238   G4bool ProjectileRemnant::addDynamicalSpectator(Particle * const p) {
00239 // assert(p->isNucleon());
00240 
00241     // Add the initial (off-shell) momentum and energy to the projectile remnant
00242     ThreeVector const &oldMomentum = getStoredMomentum(p);
00243     const ThreeVector theNewMomentum = theMomentum + oldMomentum;
00244     const G4double oldEnergy = p->getEnergy();
00245     const G4double theNewEnergy = theEnergy + oldEnergy;
00246 
00247     // Check that the excitation energy of the new projectile remnant is non-negative
00248     const G4double theNewMass = ParticleTable::getTableMass(theA+p->getA(),theZ+p->getZ());
00249     const G4double theNewInvariantMassSquared = theNewEnergy*theNewEnergy-theNewMomentum.mag2();
00250 
00251     if(theNewInvariantMassSquared<0.)
00252       return false;
00253 
00254     const G4double theNewInvariantMass = std::sqrt(theNewInvariantMassSquared);
00255 
00256     if(theNewInvariantMass-theNewMass<-1.e-5)
00257       return false; // negative excitation energy here
00258 
00259     // Add the spectator to the projectile remnant
00260     theA += p->getA();
00261     theZ += p->getZ();
00262     theMomentum = theNewMomentum;
00263     theEnergy = theNewEnergy;
00264     particles.push_back(p);
00265     return true;
00266   }
00267 
00268   G4double ProjectileRemnant::computeExcitationEnergy(const long exceptID) const {
00269     // The ground-state energy is the sum of the A smallest initial projectile
00270     // energies.
00271     // For the last nucleon, return 0 so that the algorithm will just put it on
00272     // shell.
00273     if(theA==1)
00274       return 0.;
00275 
00276     const G4double groundState = theGroundStateEnergies.at(theA-2);
00277 
00278     // Compute the sum of the presently occupied energy levels
00279     const EnergyLevels theEnergyLevels = getPresentEnergyLevels(exceptID);
00280     const G4double excitedState = std::accumulate(
00281         theEnergyLevels.begin(),
00282         theEnergyLevels.end(),
00283         0.);
00284 
00285     return excitedState-groundState;
00286   }
00287 
00288 }
00289 

Generated on Mon May 27 17:48:37 2013 for Geant4 by  doxygen 1.4.7