| Geant4 Cross Reference |
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"
42 #include "G4INCLNuclearDensityFactory.hh"
43 #include "G4INCLParticleSampler.hh"
44 #include "G4INCLAllocationPool.hh"
45
46 namespace G4INCL {
47
48 /**
49 * Cluster is a particle (inherits from the Particle class) that is
50 * actually a collection of elementary particles.
51 */
52 class Cluster : public Particle {
53 public:
54
55 /** \brief Standard Cluster constructor
56 *
57 * This constructor should mainly be used when constructing Nucleus or
58 * when constructing Clusters to be used as composite projectiles.
59 */
60 Cluster(const G4int Z, const G4int A, const G4int S, const G4bool createParticleSampler=true) :
61 Particle(),
62 theExcitationEnergy(0.),
63 theSpin(0.,0.,0.),
64 theParticleSampler(NULL)
65 {
66 setType(Composite);
67 theZ = Z;
68 theA = A;
69 theS = S;
70 setINCLMass();
71 if(createParticleSampler)
72 theParticleSampler = new ParticleSampler(A,Z,S);
73 }
74
75 /**
76 * A cluster can be directly built from a list of particles.
77 */
78 template<class Iterator>
79 Cluster(Iterator begin, Iterator end) :
80 Particle(),
81 theExcitationEnergy(0.),
82 theSpin(0.,0.,0.),
83 theParticleSampler(NULL)
84 {
85 setType(Composite);
86 for(Iterator i = begin; i != end; ++i) {
87 addParticle(*i);
88 }
89 thePosition /= theA;
90 setINCLMass();
91 adjustMomentumFromEnergy();
92 }
93
94 virtual ~Cluster() {
95 delete theParticleSampler;
96 }
97
98 /// \brief Copy constructor
99 Cluster(const Cluster &rhs) :
100 Particle(rhs),
101 theExcitationEnergy(rhs.theExcitationEnergy),
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)
108 theParticleSampler = new ParticleSampler(rhs.theA,rhs.theZ,rhs.theS);
109 else
110 theParticleSampler = NULL;
111 }
112
113 /// \brief Assignment operator
114 Cluster &operator=(const Cluster &rhs) {
115 Cluster temporaryCluster(rhs);
116 Particle::operator=(temporaryCluster);
117 swap(temporaryCluster);
118 return *this;
119 }
120
121 /// \brief Helper method for the assignment operator
122 void swap(Cluster &rhs) {
123 Particle::swap(rhs);
124 std::swap(theExcitationEnergy, rhs.theExcitationEnergy);
125 std::swap(theSpin, rhs.theSpin);
126 // std::swap is overloaded by std::list and guaranteed to operate in
127 // constant time
128 std::swap(particles, rhs.particles);
129 std::swap(theParticleSampler, rhs.theParticleSampler);
130 }
131
132 ParticleSpecies getSpecies() const {
133 return ParticleSpecies(theA, theZ, theS);
134 }
135
136 void deleteParticles() {
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
145 /// \brief Set the charge number of the cluster
146 void setZ(const G4int Z) { theZ = Z; }
147
148 /// \brief Set the mass number of the cluster
149 void setA(const G4int A) { theA = A; }
150
151 /// \brief Set the strangess number of the cluster
152 void setS(const G4int S) { theS = S; }
153
154 /// \brief Get the excitation energy of the cluster.
155 G4double getExcitationEnergy() const { return theExcitationEnergy; }
156
157 /// \brief Set the excitation energy of the cluster.
158 void setExcitationEnergy(const G4double e) { theExcitationEnergy=e; }
159
160 /** \brief Get the real particle mass.
161 *
162 * Overloads the Particle method.
163 */
164 inline virtual G4double getTableMass() const { return getRealMass(); }
165
166 /**
167 * Get the list of particles in the cluster.
168 */
169 ParticleList const &getParticles() const { return particles; }
170
171 /// \brief Remove a particle from the cluster components.
172 void removeParticle(Particle * const p) { particles.remove(p); }
173
174 /**
175 * Add one particle to the cluster. This updates the cluster mass,
176 * energy, size, etc.
177 */
178 void addParticle(Particle * const p) {
179 particles.push_back(p);
180 theEnergy += p->getEnergy();
181 thePotentialEnergy += p->getPotentialEnergy();
182 theMomentum += p->getMomentum();
183 thePosition += p->getPosition();
184 theA += p->getA();
185 theZ += p->getZ();
186 theS += p->getS();
187 nCollisions += p->getNumberOfCollisions();
188 }
189
190 /// \brief Set total cluster mass, energy, size, etc. from the particles
191 void updateClusterParameters() {
192 theEnergy = 0.;
193 thePotentialEnergy = 0.;
194 theMomentum = ThreeVector();
195 thePosition = ThreeVector();
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
212 /// \brief Add a list of particles to the cluster
213 void addParticles(ParticleList const &pL) {
214 particles = pL;
215 updateClusterParameters();
216 }
217
218 /// \brief Returns the list of particles that make up the cluster
219 ParticleList getParticleList() const { return particles; }
220
221 std::string print() const {
222 std::stringstream ss;
223 ss << "Cluster (ID = " << ID << ") type = ";
224 ss << ParticleTable::getName(theType);
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
245 /// \brief Initialise the NuclearDensity pointer and sample the particles
246 virtual void initializeParticles();
247
248 /** \brief Boost to the CM of the component particles
249 *
250 * The position of all particles in the particles list is shifted so that
251 * their centre of mass is in the origin and their total momentum is
252 * zero.
253 */
254 void internalBoostToCM() {
255
256 // First compute the current CM position and total momentum
257 ThreeVector theCMPosition, theTotalMomentum;
258 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
259 theCMPosition += (*p)->getPosition();
260 theTotalMomentum += (*p)->getMomentum();
261 //theTotalEnergy += (*p)->getEnergy();
262 }
263 theCMPosition /= theA;
264 // assert((unsigned int)theA==particles.size());
265
266 // Now determine the CM velocity of the particles
267 // commented out because currently unused, see below
268 // ThreeVector betaCM = theTotalMomentum / theTotalEnergy;
269
270 // The new particle positions and momenta are scaled by a factor of
271 // \f$\sqrt{A/(A-1)}\f$, so that the resulting density distributions in
272 // the CM have the same variance as the one we started with.
273 const G4double rescaling = std::sqrt(((G4double)theA)/((G4double)(theA-1)));
274
275 // Loop again to boost and reposition
276 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
277 // \bug{We should do the following, but the Fortran version actually
278 // does not!
279 // (*p)->boost(betaCM);
280 // Here is what the Fortran version does:}
281 (*p)->setMomentum(((*p)->getMomentum()-theTotalMomentum/theA)*rescaling);
282
283 // Set the CM position of the particles
284 (*p)->setPosition(((*p)->getPosition()-theCMPosition)*rescaling);
285 }
286
287 // Set the global cluster kinematic variables
288 thePosition.setX(0.0);
289 thePosition.setY(0.0);
290 thePosition.setZ(0.0);
291 theMomentum.setX(0.0);
292 theMomentum.setY(0.0);
293 theMomentum.setZ(0.0);
294 theEnergy = getMass();
295
296 INCL_DEBUG("Cluster boosted to internal CM:" << '\n' << print());
297
298 }
299
300 /** \brief Put the cluster components off shell
301 *
302 * The Cluster components are put off shell in such a way that their total
303 * energy equals the cluster mass.
304 */
305 void putParticlesOffShell() {
306 // Compute the dynamical potential
307 const G4double theDynamicalPotential = computeDynamicalPotential();
308 INCL_DEBUG("The dynamical potential is " << theDynamicalPotential << " MeV" << '\n');
309
310 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
311 const G4double energy = (*p)->getEnergy() - theDynamicalPotential;
312 const ThreeVector &momentum = (*p)->getMomentum();
313 // Here particles are put off-shell so that we can satisfy the energy-
314 // and momentum-conservation laws
315 (*p)->setEnergy(energy);
316 (*p)->setMass(std::sqrt(energy*energy - momentum.mag2()));
317 }
318 INCL_DEBUG("Cluster components are now off shell:" << '\n'
319 << print());
320 }
321
322 /** \brief Set the position of the cluster
323 *
324 * This overloads the Particle method to take into account that the
325 * positions of the cluster members must be updated as well.
326 */
327 void setPosition(const ThreeVector &position) {
328 ThreeVector shift(position-thePosition);
329 Particle::setPosition(position);
330 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
331 (*p)->setPosition((*p)->getPosition()+shift);
332 }
333 }
334
335 /** \brief Boost the cluster with the indicated velocity
336 *
337 * The Cluster is boosted as a whole, just like any Particle object;
338 * moreover, the internal components (particles list) are also boosted,
339 * according to Alain Boudard's off-shell recipe.
340 *
341 * \param aBoostVector the velocity to boost to [c]
342 */
343 void boost(const ThreeVector &aBoostVector) {
344 Particle::boost(aBoostVector);
345 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
346 (*p)->boost(aBoostVector);
347 // Apply Lorentz contraction to the particle position
348 (*p)->lorentzContract(aBoostVector,thePosition);
349 (*p)->rpCorrelate();
350 }
351
352 INCL_DEBUG("Cluster was boosted with (bx,by,bz)=("
353 << aBoostVector.getX() << ", " << aBoostVector.getY() << ", " << aBoostVector.getZ() << "):"
354 << '\n' << print());
355
356 }
357
358 /** \brief Freeze the internal motion of the particles
359 *
360 * Each particle is assigned a frozen momentum four-vector determined by
361 * the collective cluster velocity. This is used for propagation, but not
362 * for dynamics. Normal propagation is restored by calling the
363 * Particle::thawPropagation() method, which should be done in
364 * InteractionAvatar::postInteraction.
365 */
366 void freezeInternalMotion() {
367 const ThreeVector &normMomentum = theMomentum / getMass();
368 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
369 const G4double pMass = (*p)->getMass();
370 const ThreeVector frozenMomentum = normMomentum * pMass;
371 const G4double frozenEnergy = std::sqrt(frozenMomentum.mag2()+pMass*pMass);
372 (*p)->setFrozenMomentum(frozenMomentum);
373 (*p)->setFrozenEnergy(frozenEnergy);
374 (*p)->freezePropagation();
375 }
376 }
377
378 /** \brief Rotate position of all the particles
379 *
380 * This includes the cluster components. Overloads Particle::rotateMomentum().
381 *
382 * \param angle the rotation angle
383 * \param axis a unit vector representing the rotation axis
384 */
385 virtual void rotatePosition(const G4double angle, const ThreeVector &axis);
386
387 /** \brief Rotate momentum of all the particles
388 *
389 * This includes the cluster components. Overloads Particle::rotateMomentum().
390 *
391 * \param angle the rotation angle
392 * \param axis a unit vector representing the rotation axis
393 */
394 virtual void rotateMomentum(const G4double angle, const ThreeVector &axis);
395
396 /// \brief Make all the components projectile spectators, too
397 virtual void makeProjectileSpectator() {
398 Particle::makeProjectileSpectator();
399 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
400 (*p)->makeProjectileSpectator();
401 }
402 }
403
404 /// \brief Make all the components target spectators, too
405 virtual void makeTargetSpectator() {
406 Particle::makeTargetSpectator();
407 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
408 (*p)->makeTargetSpectator();
409 }
410 }
411
412 /// \brief Make all the components participants, too
413 virtual void makeParticipant() {
414 Particle::makeParticipant();
415 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
416 (*p)->makeParticipant();
417 }
418 }
419
420 /// \brief Get the spin of the nucleus.
421 ThreeVector const &getSpin() const { return theSpin; }
422
423 /// \brief Set the spin of the nucleus.
424 void setSpin(const ThreeVector &j) { theSpin = j; }
425
426 /// \brief Get the total angular momentum (orbital + spin)
427 G4INCL::ThreeVector getAngularMomentum() const {
428 return Particle::getAngularMomentum() + getSpin();
429 }
430
431 private:
432 /** \brief Compute the dynamical cluster potential
433 *
434 * Alain Boudard's boost prescription for low-energy beams requires to
435 * define a "dynamical potential" that allows us to conserve momentum and
436 * energy when boosting the projectile cluster.
437 */
438 G4double computeDynamicalPotential() {
439 G4double theDynamicalPotential = 0.0;
440 for(ParticleIter p=particles.begin(), e=particles.end(); p!=e; ++p) {
441 theDynamicalPotential += (*p)->getEnergy();
442 }
443 theDynamicalPotential -= getTableMass();
444 theDynamicalPotential /= theA;
445
446 return theDynamicalPotential;
447 }
448
449 protected:
450 ParticleList particles;
451 G4double theExcitationEnergy;
452 ThreeVector theSpin;
453 ParticleSampler *theParticleSampler;
454
455 INCL_DECLARE_ALLOCATION_POOL(Cluster)
456 };
457
458 }
459
460 #endif
461