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
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
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
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
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) {
00099
00100
00101 const G4double theProjectileCorrectionPerNucleon = theProjectileCorrection / particles.size();
00102
00103
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
00118
00119 DEBUG("After Particle removal, the ProjectileRemnant looks like this:"
00120 << std::endl << print());
00121 }
00122
00123 ParticleList ProjectileRemnant::addDynamicalSpectators(ParticleList pL) {
00124
00125
00126
00127
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
00145
00146
00147
00148
00149
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
00156
00157 theNewMomentum += getStoredMomentum(*p);
00158 theNewEnergy += (*p)->getEnergy();
00159 theNewA += (*p)->getA();
00160 theNewZ += (*p)->getZ();
00161 }
00162
00163
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
00174
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
00184
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
00197
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
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
00222 positiveExcitationEnergy = true;
00223 }
00224 }
00225
00226
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
00240
00241
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
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;
00258
00259
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
00270
00271
00272
00273 if(theA==1)
00274 return 0.;
00275
00276 const G4double groundState = theGroundStateEnergies.at(theA-2);
00277
00278
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