Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCluster.hh
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // INCL++ intra-nuclear cascade model
27 // Alain Boudard, CEA-Saclay, France
28 // Joseph Cugnon, University of Liege, Belgium
29 // Jean-Christophe David, CEA-Saclay, France
30 // Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
31 // Sylvie Leray, CEA-Saclay, France
32 // Davide Mancusi, CEA-Saclay, France
33 //
34 #define INCLXX_IN_GEANT4_MODE 1
35 
36 #include "globals.hh"
37 
38 #ifndef G4INCLCluster_hh
39 #define G4INCLCluster_hh 1
40 
41 #include "G4INCLParticle.hh"
43 #include "G4INCLParticleSampler.hh"
44 #include "G4INCLAllocationPool.hh"
45 
46 namespace G4INCL {
47 
52  class Cluster : public Particle {
53  public:
54 
60  Cluster(const G4int Z, const G4int A, const G4bool createParticleSampler=true) :
61  Particle(),
63  theSpin(0.,0.,0.),
64  theParticleSampler(NULL)
65  {
67  theZ = Z;
68  theA = A;
69  setINCLMass();
70  if(createParticleSampler)
72  }
73 
77  template<class Iterator>
78  Cluster(Iterator begin, Iterator end) :
79  Particle(),
81  theSpin(0.,0.,0.),
82  theParticleSampler(NULL)
83  {
85  for(Iterator i = begin; i != end; ++i) {
86  addParticle(*i);
87  }
88  thePosition /= theA;
89  setINCLMass();
91  }
92 
93  virtual ~Cluster() {
94  delete theParticleSampler;
95  }
96 
98  Cluster(const Cluster &rhs) :
99  Particle(rhs),
101  theSpin(rhs.theSpin)
102  {
103  for(ParticleIter p=rhs.particles.begin(), e=rhs.particles.end(); p!=e; ++p) {
104  particles.push_back(new Particle(**p));
105  }
106  if(rhs.theParticleSampler)
108  else
109  theParticleSampler = NULL;
110  }
111 
113  Cluster &operator=(const Cluster &rhs) {
114  Cluster temporaryCluster(rhs);
115  Particle::operator=(temporaryCluster);
116  swap(temporaryCluster);
117  return *this;
118  }
119 
121  void swap(Cluster &rhs) {
122  Particle::swap(rhs);
124  std::swap(theSpin, rhs.theSpin);
125  // std::swap is overloaded by std::list and guaranteed to operate in
126  // constant time
127  std::swap(particles, rhs.particles);
128  std::swap(theParticleSampler, rhs.theParticleSampler);
129  }
130 
132  return ParticleSpecies(theA, theZ);
133  }
134 
136  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
137  delete (*p);
138  }
139  clearParticles();
140  }
141 
142  void clearParticles() { particles.clear(); }
143 
145  void setZ(const G4int Z) { theZ = Z; }
146 
148  void setA(const G4int A) { theA = A; }
149 
152 
155 
160  inline virtual G4double getTableMass() const { return getRealMass(); }
161 
165  ParticleList const &getParticles() const { return particles; }
166 
168  void removeParticle(Particle * const p) { particles.remove(p); }
169 
174  void addParticle(Particle * const p) {
175  particles.push_back(p);
176  theEnergy += p->getEnergy();
178  theMomentum += p->getMomentum();
179  thePosition += p->getPosition();
180  theA += p->getA();
181  theZ += p->getZ();
183  }
184 
187  theEnergy = 0.;
188  thePotentialEnergy = 0.;
191  theA = 0;
192  theZ = 0;
193  nCollisions = 0;
194  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
195  theEnergy += (*p)->getEnergy();
196  thePotentialEnergy += (*p)->getPotentialEnergy();
197  theMomentum += (*p)->getMomentum();
198  thePosition += (*p)->getPosition();
199  theA += (*p)->getA();
200  theZ += (*p)->getZ();
201  nCollisions += (*p)->getNumberOfCollisions();
202  }
203  }
204 
206  void addParticles(ParticleList const &pL) {
207  particles = pL;
209  }
210 
213 
214  std::string print() const {
215  std::stringstream ss;
216  ss << "Cluster (ID = " << ID << ") type = ";
218  ss << '\n'
219  << " A = " << theA << '\n'
220  << " Z = " << theZ << '\n'
221  << " mass = " << getMass() << '\n'
222  << " energy = " << theEnergy << '\n'
223  << " momentum = "
224  << theMomentum.print()
225  << '\n'
226  << " position = "
227  << thePosition.print()
228  << '\n'
229  << "Contains the following particles:"
230  << '\n';
231  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i)
232  ss << (*i)->print();
233  ss << '\n';
234  return ss.str();
235  }
236 
238  virtual void initializeParticles();
239 
247 
248  // First compute the current CM position and total momentum
249  ThreeVector theCMPosition, theTotalMomentum;
250  G4double theTotalEnergy = 0.0;
251  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
252  theCMPosition += (*p)->getPosition();
253  theTotalMomentum += (*p)->getMomentum();
254  theTotalEnergy += (*p)->getEnergy();
255  }
256  theCMPosition /= theA;
257 // assert((unsigned int)theA==particles.size());
258 
259  // Now determine the CM velocity of the particles
260  // commented out because currently unused, see below
261  // ThreeVector betaCM = theTotalMomentum / theTotalEnergy;
262 
263  // The new particle positions and momenta are scaled by a factor of
264  // \f$\sqrt{A/(A-1)}\f$, so that the resulting density distributions in
265  // the CM have the same variance as the one we started with.
266  const G4double rescaling = std::sqrt(((G4double)theA)/((G4double)(theA-1)));
267 
268  // Loop again to boost and reposition
269  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
270  // \bug{We should do the following, but the Fortran version actually
271  // does not!
272  // (*p)->boost(betaCM);
273  // Here is what the Fortran version does:}
274  (*p)->setMomentum(((*p)->getMomentum()-theTotalMomentum/theA)*rescaling);
275 
276  // Set the CM position of the particles
277  (*p)->setPosition(((*p)->getPosition()-theCMPosition)*rescaling);
278  }
279 
280  // Set the global cluster kinematic variables
281  thePosition.setX(0.0);
282  thePosition.setY(0.0);
283  thePosition.setZ(0.0);
284  theMomentum.setX(0.0);
285  theMomentum.setY(0.0);
286  theMomentum.setZ(0.0);
287  theEnergy = getMass();
288 
289  INCL_DEBUG("Cluster boosted to internal CM:" << '\n' << print());
290 
291  }
292 
299  // Compute the dynamical potential
300  const G4double theDynamicalPotential = computeDynamicalPotential();
301  INCL_DEBUG("The dynamical potential is " << theDynamicalPotential << " MeV" << '\n');
302 
303  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
304  const G4double energy = (*p)->getEnergy() - theDynamicalPotential;
305  const ThreeVector &momentum = (*p)->getMomentum();
306  // Here particles are put off-shell so that we can satisfy the energy-
307  // and momentum-conservation laws
308  (*p)->setEnergy(energy);
309  (*p)->setMass(std::sqrt(energy*energy - momentum.mag2()));
310  }
311  INCL_DEBUG("Cluster components are now off shell:" << '\n'
312  << print());
313  }
314 
321  ThreeVector shift(position-thePosition);
322  Particle::setPosition(position);
323  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
324  (*p)->setPosition((*p)->getPosition()+shift);
325  }
326  }
327 
336  void boost(const ThreeVector &aBoostVector) {
337  Particle::boost(aBoostVector);
338  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
339  (*p)->boost(aBoostVector);
340  // Apply Lorentz contraction to the particle position
341  (*p)->lorentzContract(aBoostVector,thePosition);
342  (*p)->rpCorrelate();
343  }
344 
345  INCL_DEBUG("Cluster was boosted with (bx,by,bz)=("
346  << aBoostVector.getX() << ", " << aBoostVector.getY() << ", " << aBoostVector.getZ() << "):"
347  << '\n' << print());
348 
349  }
350 
360  const ThreeVector &normMomentum = theMomentum / getMass();
361  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
362  const G4double pMass = (*p)->getMass();
363  const ThreeVector frozenMomentum = normMomentum * pMass;
364  const G4double frozenEnergy = std::sqrt(frozenMomentum.mag2()+pMass*pMass);
365  (*p)->setFrozenMomentum(frozenMomentum);
366  (*p)->setFrozenEnergy(frozenEnergy);
367  (*p)->freezePropagation();
368  }
369  }
370 
378  virtual void rotatePosition(const G4double angle, const ThreeVector &axis);
379 
387  virtual void rotateMomentum(const G4double angle, const ThreeVector &axis);
388 
390  virtual void makeProjectileSpectator() {
392  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
393  (*p)->makeProjectileSpectator();
394  }
395  }
396 
398  virtual void makeTargetSpectator() {
400  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
401  (*p)->makeTargetSpectator();
402  }
403  }
404 
406  virtual void makeParticipant() {
408  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
409  (*p)->makeParticipant();
410  }
411  }
412 
414  ThreeVector const &getSpin() const { return theSpin; }
415 
417  void setSpin(const ThreeVector &j) { theSpin = j; }
418 
422  }
423 
424  private:
431  G4double computeDynamicalPotential() {
432  G4double theDynamicalPotential = 0.0;
433  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
434  theDynamicalPotential += (*p)->getEnergy();
435  }
436  theDynamicalPotential -= getTableMass();
437  theDynamicalPotential /= theA;
438 
439  return theDynamicalPotential;
440  }
441 
442  protected:
447 
449  };
450 
451 }
452 
453 #endif
G4int getA() const
Returns the baryon number.
#define INCL_DECLARE_ALLOCATION_POOL(T)
virtual void rotateMomentum(const G4double angle, const ThreeVector &axis)
Rotate momentum of all the particles.
G4int getNumberOfCollisions() const
Return the number of collisions undergone by the particle.
virtual void makeTargetSpectator()
virtual G4double getTableMass() const
Get the real particle mass.
G4double getMass() const
Get the cached particle mass.
void putParticlesOffShell()
Put the cluster components off shell.
ParticleList const & getParticles() const
virtual void makeTargetSpectator()
Make all the components target spectators, too.
const char * p
Definition: xmltok.h:285
void setZ(const G4int Z)
Set the charge number of the cluster.
void boost(const ThreeVector &aBoostVector)
const G4INCL::ThreeVector & getMomentum() const
static G4double angle[DIM]
G4double getEnergy() const
void setINCLMass()
Set the mass of the Particle to its table mass.
virtual void makeProjectileSpectator()
Singleton for recycling allocation of instances of a given class.
void boost(const ThreeVector &aBoostVector)
Boost the cluster with the indicated velocity.
virtual ~Cluster()
void updateClusterParameters()
Set total cluster mass, energy, size, etc. from the particles.
int G4int
Definition: G4Types.hh:78
G4double mag2() const
G4double getPotentialEnergy() const
Get the particle potential energy.
void setY(G4double ay)
Set the y coordinate.
G4INCL::ThreeVector getAngularMomentum() const
Get the total angular momentum (orbital + spin)
Class for sampling particles in a nucleus.
double A(double temperature)
ThreeVector theSpin
std::string print() const
void swap(Cluster &rhs)
Helper method for the assignment operator.
bool G4bool
Definition: G4Types.hh:79
virtual void initializeParticles()
Initialise the NuclearDensity pointer and sample the particles.
virtual G4INCL::ThreeVector getAngularMomentum() const
Cluster(const Cluster &rhs)
Copy constructor.
ParticleList particles
virtual void setPosition(const G4INCL::ThreeVector &position)
G4int getZ() const
Returns the charge number.
void freezeInternalMotion()
Freeze the internal motion of the particles.
G4INCL::ThreeVector thePosition
virtual void makeParticipant()
Make all the components participants, too.
ParticleSampler * theParticleSampler
void removeParticle(Particle *const p)
Remove a particle from the cluster components.
G4double theExcitationEnergy
void setA(const G4int A)
Set the mass number of the cluster.
ParticleList getParticleList() const
Returns the list of particles that make up the cluster.
G4INCL::ParticleType theType
Cluster & operator=(const Cluster &rhs)
Assignment operator.
G4double energy(const ThreeVector &p, const G4double m)
void addParticles(ParticleList const &pL)
Add a list of particles to the cluster.
const G4INCL::ThreeVector & getPosition() const
void swap(Particle &rhs)
Helper method for the assignment operator.
virtual void rotatePosition(const G4double angle, const ThreeVector &axis)
Rotate position of all the particles.
void setType(ParticleType t)
void internalBoostToCM()
Boost to the CM of the component particles.
ParticleSpecies getSpecies() const
Get the particle species.
Cluster(Iterator begin, Iterator end)
Particle & operator=(const Particle &rhs)
Assignment operator.
void addParticle(Particle *const p)
G4INCL::ThreeVector theMomentum
ThreeVector const & getSpin() const
Get the spin of the nucleus.
G4double getX() const
std::string print() const
void setX(G4double ax)
Set the x coordinate.
G4double getExcitationEnergy() const
Get the excitation energy of the cluster.
G4double getRealMass() const
Get the real particle mass.
virtual void makeProjectileSpectator()
Make all the components projectile spectators, too.
virtual void makeParticipant()
double G4double
Definition: G4Types.hh:76
void setZ(G4double az)
Set the z coordinate.
#define INCL_DEBUG(x)
const ThreeVector & adjustMomentumFromEnergy()
Rescale the momentum to match the total energy.
void setExcitationEnergy(const G4double e)
Set the excitation energy of the cluster.
void setSpin(const ThreeVector &j)
Set the spin of the nucleus.
G4double thePotentialEnergy
void setPosition(const ThreeVector &position)
Set the position of the cluster.
ParticleList::const_iterator ParticleIter
G4double getZ() const
std::string getName(const ParticleType t)
Get the native INCL name of the particle.
G4double getY() const
Cluster(const G4int Z, const G4int A, const G4bool createParticleSampler=true)
Standard Cluster constructor.