ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLCluster.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4INCLCluster.hh
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 G4int S, const G4bool createParticleSampler=true) :
61  Particle(),
63  theSpin(0.,0.,0.),
64  theParticleSampler(NULL)
65  {
67  theZ = Z;
68  theA = A;
69  theS = S;
70  setINCLMass();
71  if(createParticleSampler)
73  }
74 
78  template<class Iterator>
79  Cluster(Iterator begin, Iterator end) :
80  Particle(),
82  theSpin(0.,0.,0.),
83  theParticleSampler(NULL)
84  {
86  for(Iterator i = begin; i != end; ++i) {
87  addParticle(*i);
88  }
89  thePosition /= theA;
90  setINCLMass();
92  }
93 
94  virtual ~Cluster() {
95  delete theParticleSampler;
96  }
97 
99  Cluster(const Cluster &rhs) :
100  Particle(rhs),
102  theSpin(rhs.theSpin)
103  {
104  for(ParticleIter p=rhs.particles.begin(), e=rhs.particles.end(); p!=e; ++p) {
105  particles.push_back(new Particle(**p));
106  }
107  if(rhs.theParticleSampler)
109  else
110  theParticleSampler = NULL;
111  }
112 
114  Cluster &operator=(const Cluster &rhs) {
115  Cluster temporaryCluster(rhs);
116  Particle::operator=(temporaryCluster);
117  swap(temporaryCluster);
118  return *this;
119  }
120 
122  void swap(Cluster &rhs) {
123  Particle::swap(rhs);
125  std::swap(theSpin, rhs.theSpin);
126  // std::swap is overloaded by std::list and guaranteed to operate in
127  // constant time
130  }
131 
133  return ParticleSpecies(theA, theZ, theS);
134  }
135 
137  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
138  delete (*p);
139  }
140  clearParticles();
141  }
142 
143  void clearParticles() { particles.clear(); }
144 
146  void setZ(const G4int Z) { theZ = Z; }
147 
149  void setA(const G4int A) { theA = A; }
150 
152  void setS(const G4int S) { theS = S; }
153 
156 
159 
164  inline virtual G4double getTableMass() const { return getRealMass(); }
165 
169  ParticleList const &getParticles() const { return particles; }
170 
172  void removeParticle(Particle * const p) { particles.remove(p); }
173 
178  void addParticle(Particle * const p) {
179  particles.push_back(p);
180  theEnergy += p->getEnergy();
182  theMomentum += p->getMomentum();
183  thePosition += p->getPosition();
184  theA += p->getA();
185  theZ += p->getZ();
186  theS += p->getS();
188  }
189 
192  theEnergy = 0.;
193  thePotentialEnergy = 0.;
196  theA = 0;
197  theZ = 0;
198  theS = 0;
199  nCollisions = 0;
200  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
201  theEnergy += (*p)->getEnergy();
202  thePotentialEnergy += (*p)->getPotentialEnergy();
203  theMomentum += (*p)->getMomentum();
204  thePosition += (*p)->getPosition();
205  theA += (*p)->getA();
206  theZ += (*p)->getZ();
207  theS += (*p)->getS();
208  nCollisions += (*p)->getNumberOfCollisions();
209  }
210  }
211 
213  void addParticles(ParticleList const &pL) {
214  particles = pL;
216  }
217 
220 
221  std::string print() const {
222  std::stringstream ss;
223  ss << "Cluster (ID = " << ID << ") type = ";
225  ss << '\n'
226  << " A = " << theA << '\n'
227  << " Z = " << theZ << '\n'
228  << " S = " << theS << '\n'
229  << " mass = " << getMass() << '\n'
230  << " energy = " << theEnergy << '\n'
231  << " momentum = "
232  << theMomentum.print()
233  << '\n'
234  << " position = "
235  << thePosition.print()
236  << '\n'
237  << "Contains the following particles:"
238  << '\n';
239  for(ParticleIter i=particles.begin(), e=particles.end(); i!=e; ++i)
240  ss << (*i)->print();
241  ss << '\n';
242  return ss.str();
243  }
244 
246  virtual void initializeParticles();
247 
255 
256  // First compute the current CM position and total momentum
257  ThreeVector theCMPosition, theTotalMomentum;
258  G4double theTotalEnergy = 0.0;
259  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
260  theCMPosition += (*p)->getPosition();
261  theTotalMomentum += (*p)->getMomentum();
262  theTotalEnergy += (*p)->getEnergy();
263  }
264  theCMPosition /= theA;
265 // assert((unsigned int)theA==particles.size());
266 
267  // Now determine the CM velocity of the particles
268  // commented out because currently unused, see below
269  // ThreeVector betaCM = theTotalMomentum / theTotalEnergy;
270 
271  // The new particle positions and momenta are scaled by a factor of
272  // \f$\sqrt{A/(A-1)}\f$, so that the resulting density distributions in
273  // the CM have the same variance as the one we started with.
274  const G4double rescaling = std::sqrt(((G4double)theA)/((G4double)(theA-1)));
275 
276  // Loop again to boost and reposition
277  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
278  // \bug{We should do the following, but the Fortran version actually
279  // does not!
280  // (*p)->boost(betaCM);
281  // Here is what the Fortran version does:}
282  (*p)->setMomentum(((*p)->getMomentum()-theTotalMomentum/theA)*rescaling);
283 
284  // Set the CM position of the particles
285  (*p)->setPosition(((*p)->getPosition()-theCMPosition)*rescaling);
286  }
287 
288  // Set the global cluster kinematic variables
289  thePosition.setX(0.0);
290  thePosition.setY(0.0);
291  thePosition.setZ(0.0);
292  theMomentum.setX(0.0);
293  theMomentum.setY(0.0);
294  theMomentum.setZ(0.0);
295  theEnergy = getMass();
296 
297  INCL_DEBUG("Cluster boosted to internal CM:" << '\n' << print());
298 
299  }
300 
307  // Compute the dynamical potential
308  const G4double theDynamicalPotential = computeDynamicalPotential();
309  INCL_DEBUG("The dynamical potential is " << theDynamicalPotential << " MeV" << '\n');
310 
311  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
312  const G4double energy = (*p)->getEnergy() - theDynamicalPotential;
313  const ThreeVector &momentum = (*p)->getMomentum();
314  // Here particles are put off-shell so that we can satisfy the energy-
315  // and momentum-conservation laws
316  (*p)->setEnergy(energy);
317  (*p)->setMass(std::sqrt(energy*energy - momentum.mag2()));
318  }
319  INCL_DEBUG("Cluster components are now off shell:" << '\n'
320  << print());
321  }
322 
329  ThreeVector shift(position-thePosition);
330  Particle::setPosition(position);
331  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
332  (*p)->setPosition((*p)->getPosition()+shift);
333  }
334  }
335 
344  void boost(const ThreeVector &aBoostVector) {
345  Particle::boost(aBoostVector);
346  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
347  (*p)->boost(aBoostVector);
348  // Apply Lorentz contraction to the particle position
349  (*p)->lorentzContract(aBoostVector,thePosition);
350  (*p)->rpCorrelate();
351  }
352 
353  INCL_DEBUG("Cluster was boosted with (bx,by,bz)=("
354  << aBoostVector.getX() << ", " << aBoostVector.getY() << ", " << aBoostVector.getZ() << "):"
355  << '\n' << print());
356 
357  }
358 
368  const ThreeVector &normMomentum = theMomentum / getMass();
369  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
370  const G4double pMass = (*p)->getMass();
371  const ThreeVector frozenMomentum = normMomentum * pMass;
372  const G4double frozenEnergy = std::sqrt(frozenMomentum.mag2()+pMass*pMass);
373  (*p)->setFrozenMomentum(frozenMomentum);
374  (*p)->setFrozenEnergy(frozenEnergy);
375  (*p)->freezePropagation();
376  }
377  }
378 
386  virtual void rotatePosition(const G4double angle, const ThreeVector &axis);
387 
395  virtual void rotateMomentum(const G4double angle, const ThreeVector &axis);
396 
398  virtual void makeProjectileSpectator() {
400  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
401  (*p)->makeProjectileSpectator();
402  }
403  }
404 
406  virtual void makeTargetSpectator() {
408  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
409  (*p)->makeTargetSpectator();
410  }
411  }
412 
414  virtual void makeParticipant() {
416  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
417  (*p)->makeParticipant();
418  }
419  }
420 
422  ThreeVector const &getSpin() const { return theSpin; }
423 
425  void setSpin(const ThreeVector &j) { theSpin = j; }
426 
430  }
431 
432  private:
440  G4double theDynamicalPotential = 0.0;
441  for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
442  theDynamicalPotential += (*p)->getEnergy();
443  }
444  theDynamicalPotential -= getTableMass();
445  theDynamicalPotential /= theA;
446 
447  return theDynamicalPotential;
448  }
449 
450  protected:
455 
457  };
458 
459 }
460 
461 #endif