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 G4INCLCluster_hh
00038 #define G4INCLCluster_hh 1
00039
00040 #include "G4INCLParticle.hh"
00041 #include "G4INCLNuclearDensityFactory.hh"
00042 #include "G4INCLParticleSampler.hh"
00043
00044 namespace G4INCL {
00045
00050 class Cluster : public Particle {
00051 public:
00052
00058 Cluster(const G4int Z, const G4int A, const G4bool createParticleSampler=true) :
00059 Particle(),
00060 theExcitationEnergy(0.),
00061 theSpin(0.,0.,0.),
00062 theParticleSampler(NULL)
00063 {
00064 setType(Composite);
00065 theZ = Z;
00066 theA = A;
00067 setINCLMass();
00068 if(createParticleSampler)
00069 theParticleSampler = NuclearDensityFactory::createParticleSampler(A,Z);
00070 }
00071
00075 template<class Iterator>
00076 Cluster(Iterator begin, Iterator end) :
00077 Particle(),
00078 theExcitationEnergy(0.),
00079 theSpin(0.,0.,0.),
00080 theParticleSampler(NULL)
00081 {
00082 setType(Composite);
00083 for(Iterator i = begin; i != end; ++i) {
00084 addParticle(*i);
00085 }
00086 thePosition /= theA;
00087 setINCLMass();
00088 adjustMomentumFromEnergy();
00089 }
00090
00091 virtual ~Cluster() {
00092 delete theParticleSampler;
00093 }
00094
00096 Cluster(const Cluster &rhs) :
00097 Particle(rhs),
00098 theExcitationEnergy(rhs.theExcitationEnergy)
00099 {
00100 deleteParticles();
00101 for(ParticleIter p=rhs.particles.begin(); p!=rhs.particles.end(); ++p) {
00102 particles.push_back(new Particle(**p));
00103 }
00104 }
00105
00107 Cluster &operator=(const Cluster &rhs) {
00108 Cluster temporaryCluster(rhs);
00109 Particle::operator=(temporaryCluster);
00110 swap(temporaryCluster);
00111 return *this;
00112 }
00113
00115 void swap(Cluster &rhs) {
00116 Particle::swap(rhs);
00117 std::swap(theExcitationEnergy, rhs.theExcitationEnergy);
00118 std::swap(theSpin, rhs.theSpin);
00119
00120
00121 std::swap(particles, rhs.particles);
00122 std::swap(theParticleSampler, rhs.theParticleSampler);
00123 }
00124
00125 ParticleSpecies getSpecies() const {
00126 return ParticleSpecies(theA, theZ);
00127 }
00128
00129 void deleteParticles() {
00130 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00131 delete (*p);
00132 }
00133 clearParticles();
00134 }
00135
00136 void clearParticles() { particles.clear(); }
00137
00139 void setZ(const G4int Z) { theZ = Z; }
00140
00142 void setA(const G4int A) { theA = A; }
00143
00145 G4double getExcitationEnergy() const { return theExcitationEnergy; }
00146
00148 void setExcitationEnergy(const G4double e) { theExcitationEnergy=e; }
00149
00154 inline virtual G4double getTableMass() const { return getRealMass(); }
00155
00159 ParticleList const &getParticles() const { return particles; }
00160
00162 void removeParticle(Particle * const p) { particles.remove(p); }
00163
00168 void addParticle(Particle * const p) {
00169
00170 particles.push_back(p);
00171 theEnergy += p->getEnergy();
00172 thePotentialEnergy += p->getPotentialEnergy();
00173 theMomentum += p->getMomentum();
00174 thePosition += p->getPosition();
00175 theA += p->getA();
00176 theZ += p->getZ();
00177 nCollisions += p->getNumberOfCollisions();
00178 }
00179
00181 void addParticles(ParticleList const &pL) {
00182 for(ParticleIter p=pL.begin(); p!=pL.end(); ++p)
00183 addParticle(*p);
00184 }
00185
00187 ParticleList getParticleList() const { return particles; }
00188
00189 std::string print() const {
00190 std::stringstream ss;
00191 ss << "Cluster (ID = " << ID << ") type = ";
00192 ss << ParticleTable::getName(theType);
00193 ss << std::endl
00194 << " A = " << theA << std::endl
00195 << " Z = " << theZ << std::endl
00196 << " mass = " << getMass() << std::endl
00197 << " energy = " << theEnergy << std::endl
00198 << " momentum = "
00199 << theMomentum.print()
00200 << std::endl
00201 << " position = "
00202 << thePosition.print()
00203 << std::endl
00204 << "Contains the following particles:"
00205 << std::endl;
00206 for(ParticleIter i=particles.begin(); i!=particles.end(); ++i)
00207 ss << (*i)->print();
00208 ss << std::endl;
00209 return ss.str();
00210 }
00211
00213 virtual void initializeParticles();
00214
00221 void internalBoostToCM() {
00222
00223
00224 ThreeVector theCMPosition, theTotalMomentum;
00225 G4double theTotalEnergy = 0.0;
00226 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00227 theCMPosition += (*p)->getPosition();
00228 theTotalMomentum += (*p)->getMomentum();
00229 theTotalEnergy += (*p)->getEnergy();
00230 }
00231 theCMPosition /= theA;
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 const G4double rescaling = std::sqrt(((G4double)theA)/((G4double)(theA-1)));
00242
00243
00244 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00245
00246
00247
00248
00249 (*p)->setMomentum(((*p)->getMomentum()-theTotalMomentum/theA)*rescaling);
00250
00251
00252 (*p)->setPosition(((*p)->getPosition()-theCMPosition)*rescaling);
00253 }
00254
00255
00256 thePosition.setX(0.0);
00257 thePosition.setY(0.0);
00258 thePosition.setZ(0.0);
00259 theMomentum.setX(0.0);
00260 theMomentum.setY(0.0);
00261 theMomentum.setZ(0.0);
00262 theEnergy = getMass();
00263
00264 DEBUG("Cluster boosted to internal CM:" << std::endl << print());
00265
00266 }
00267
00273 void putParticlesOffShell() {
00274
00275 const G4double theDynamicalPotential = computeDynamicalPotential();
00276 DEBUG("The dynamical potential is " << theDynamicalPotential << " MeV" << std::endl);
00277
00278 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00279 const G4double energy = (*p)->getEnergy() - theDynamicalPotential;
00280 const ThreeVector &momentum = (*p)->getMomentum();
00281
00282
00283 (*p)->setEnergy(energy);
00284 (*p)->setMass(std::sqrt(energy*energy - momentum.mag2()));
00285 DEBUG("Cluster components are now off shell:" << std::endl
00286 << print());
00287 }
00288 }
00289
00295 void setPosition(const ThreeVector &position) {
00296 ThreeVector shift(position-thePosition);
00297 Particle::setPosition(position);
00298 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00299 (*p)->setPosition((*p)->getPosition()+shift);
00300 }
00301 }
00302
00311 void boost(const ThreeVector &aBoostVector) {
00312 Particle::boost(aBoostVector);
00313 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00314 (*p)->boost(aBoostVector);
00315
00316 (*p)->lorentzContract(aBoostVector,thePosition);
00317 }
00318
00319 DEBUG("Cluster was boosted with (bx,by,bz)=("
00320 << aBoostVector.getX() << ", " << aBoostVector.getY() << ", " << aBoostVector.getZ() << "):"
00321 << std::endl << print());
00322
00323 }
00324
00333 void freezeInternalMotion() {
00334 const ThreeVector &normMomentum = theMomentum / getMass();
00335 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00336 const G4double pMass = (*p)->getMass();
00337 const ThreeVector frozenMomentum = normMomentum * pMass;
00338 const G4double frozenEnergy = std::sqrt(frozenMomentum.mag2()+pMass*pMass);
00339 (*p)->setFrozenMomentum(frozenMomentum);
00340 (*p)->setFrozenEnergy(frozenEnergy);
00341 (*p)->freezePropagation();
00342 }
00343 }
00344
00352 virtual void rotate(const G4double angle, const ThreeVector &axis) {
00353 Particle::rotate(angle, axis);
00354 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00355 (*p)->rotate(angle, axis);
00356 }
00357 }
00358
00360 virtual void makeProjectileSpectator() {
00361 Particle::makeProjectileSpectator();
00362 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00363 (*p)->makeProjectileSpectator();
00364 }
00365 }
00366
00368 virtual void makeTargetSpectator() {
00369 Particle::makeTargetSpectator();
00370 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00371 (*p)->makeTargetSpectator();
00372 }
00373 }
00374
00376 virtual void makeParticipant() {
00377 Particle::makeParticipant();
00378 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00379 (*p)->makeParticipant();
00380 }
00381 }
00382
00384 ThreeVector const &getSpin() const { return theSpin; }
00385
00387 void setSpin(const ThreeVector &j) { theSpin = j; }
00388
00390 G4INCL::ThreeVector getAngularMomentum() const {
00391 return Particle::getAngularMomentum() + getSpin();
00392 }
00393
00394 private:
00401 G4double computeDynamicalPotential() {
00402 G4double theDynamicalPotential = 0.0;
00403 for(ParticleIter p=particles.begin(); p!=particles.end(); ++p) {
00404 theDynamicalPotential += (*p)->getEnergy();
00405 }
00406 theDynamicalPotential -= getTableMass();
00407 theDynamicalPotential /= theA;
00408
00409 return theDynamicalPotential;
00410 }
00411
00412 protected:
00413 ParticleList particles;
00414 G4double theExcitationEnergy;
00415 ThreeVector theSpin;
00416 ParticleSampler *theParticleSampler;
00417
00418 };
00419
00420 }
00421
00422 #endif