Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4INCLParticle.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 // Pekka Kaitaniemi, CEA and Helsinki Institute of Physics
28 // Davide Mancusi, CEA
29 // Alain Boudard, CEA
30 // Sylvie Leray, CEA
31 // Joseph Cugnon, University of Liege
32 //
33 #define INCLXX_IN_GEANT4_MODE 1
34 
35 #include "globals.hh"
36 
37 /*
38  * Particle.hh
39  *
40  * \date Jun 5, 2009
41  * \author Pekka Kaitaniemi
42  */
43 
44 #ifndef PARTICLE_HH_
45 #define PARTICLE_HH_
46 
47 #include "G4INCLThreeVector.hh"
48 #include "G4INCLParticleTable.hh"
49 #include "G4INCLParticleType.hh"
50 #include "G4INCLParticleSpecies.hh"
51 #include "G4INCLLogger.hh"
52 #include <list>
53 #include <sstream>
54 #include <string>
55 #include <algorithm>
56 
57 namespace G4INCL {
58 
59  class Particle;
60 
61  typedef std::list<G4INCL::Particle*> ParticleList;
62  typedef std::list<G4INCL::Particle*>::const_iterator ParticleIter;
63 
64  class Particle {
65  public:
66  Particle();
68  Particle(ParticleType t, ThreeVector momentum, ThreeVector position);
69  virtual ~Particle() {}
70 
75  Particle(const Particle &rhs) :
76  theZ(rhs.theZ),
77  theA(rhs.theA),
79  theType(rhs.theType),
80  theEnergy(rhs.theEnergy),
86  nDecays(rhs.nDecays),
88  theHelicity(rhs.theHelicity),
89  emissionTime(rhs.emissionTime),
90  outOfWell(rhs.outOfWell),
91  theMass(rhs.theMass)
92  {
93  if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
95  else
99  else
101  // ID intentionally not copied
102  ID = nextID++;
103  }
104 
105  protected:
107  void swap(Particle &rhs) {
108  std::swap(theZ, rhs.theZ);
109  std::swap(theA, rhs.theA);
111  std::swap(theType, rhs.theType);
112  if(rhs.thePropagationEnergy == &(rhs.theFrozenEnergy))
114  else
118  if(rhs.thePropagationMomentum == &(rhs.theFrozenMomentum))
120  else
126  std::swap(nDecays, rhs.nDecays);
128  // ID intentionally not swapped
129 
130  std::swap(theHelicity, rhs.theHelicity);
131  std::swap(emissionTime, rhs.emissionTime);
132  std::swap(outOfWell, rhs.outOfWell);
133 
134  std::swap(theMass, rhs.theMass);
135  }
136 
137  public:
138 
143  Particle &operator=(const Particle &rhs) {
144  Particle temporaryParticle(rhs);
145  swap(temporaryParticle);
146  return *this;
147  }
148 
154  return theType;
155  };
156 
159  return ParticleSpecies(theType);
160  };
161 
163  theType = t;
164  switch(theType)
165  {
166  case DeltaPlusPlus:
167  theA = 1;
168  theZ = 2;
169  break;
170  case Proton:
171  case DeltaPlus:
172  theA = 1;
173  theZ = 1;
174  break;
175  case Neutron:
176  case DeltaZero:
177  theA = 1;
178  theZ = 0;
179  break;
180  case DeltaMinus:
181  theA = 1;
182  theZ = -1;
183  break;
184  case PiPlus:
185  theA = 0;
186  theZ = 1;
187  break;
188  case PiZero:
189  theA = 0;
190  theZ = 0;
191  break;
192  case PiMinus:
193  theA = 0;
194  theZ = -1;
195  break;
196  case Composite:
197  // ERROR("Trying to set particle type to Composite! Construct a Cluster object instead" << std::endl);
198  theA = 0;
199  theZ = 0;
200  break;
201  case UnknownParticle:
202  theA = 0;
203  theZ = 0;
204  ERROR("Trying to set particle type to Unknown!" << std::endl);
205  break;
206  }
207 
208  if( !isResonance() && t!=Composite )
209  setINCLMass();
210  }
211 
215  G4bool isNucleon() const {
217  return true;
218  else
219  return false;
220  };
221 
223  return theParticipantType;
224  }
225 
228  }
229 
232  }
233 
236  }
237 
240  }
241 
242  virtual void makeParticipant() {
244  }
245 
246  virtual void makeTargetSpectator() {
248  }
249 
250  virtual void makeProjectileSpectator() {
252  }
253 
255  G4bool isPion() const { return (theType == PiPlus || theType == PiZero || theType == PiMinus); }
256 
258  inline G4bool isResonance() const { return isDelta(); }
259 
261  inline G4bool isDelta() const {
262  return (theType==DeltaPlusPlus || theType==DeltaPlus ||
264  }
265 
267  G4int getA() const { return theA; }
268 
270  G4int getZ() const { return theZ; }
271 
272  G4double getBeta() const {
273  const G4double P = theMomentum.mag();
274  return P/theEnergy;
275  }
276 
284  return theMomentum / theEnergy;
285  }
286 
293  void boost(const ThreeVector &aBoostVector) {
294  const G4double beta2 = aBoostVector.mag2();
295  const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
296  const G4double bp = theMomentum.dot(aBoostVector);
297  const G4double alpha = (gamma*gamma)/(1.0 + gamma);
298 
299  theMomentum = theMomentum + aBoostVector * (alpha * bp - gamma * theEnergy);
300  theEnergy = gamma * (theEnergy - bp);
301  }
302 
311  void lorentzContract(const ThreeVector &aBoostVector, const ThreeVector &refPos) {
312  const G4double beta2 = aBoostVector.mag2();
313  const G4double gamma = 1.0 / std::sqrt(1.0 - beta2);
314  const ThreeVector theRelativePosition = thePosition - refPos;
315  const ThreeVector transversePosition = theRelativePosition - aBoostVector * (theRelativePosition.dot(aBoostVector) / aBoostVector.mag2());
316  const ThreeVector longitudinalPosition = theRelativePosition - transversePosition;
317 
318  thePosition = refPos + transversePosition + longitudinalPosition / gamma;
319  }
320 
322  inline G4double getMass() const { return theMass; }
323 
325  inline G4double getINCLMass() const {
326  switch(theType) {
327  case Proton:
328  case Neutron:
329  case PiPlus:
330  case PiMinus:
331  case PiZero:
333  break;
334 
335  case DeltaPlusPlus:
336  case DeltaPlus:
337  case DeltaZero:
338  case DeltaMinus:
339  return theMass;
340  break;
341 
342  case Composite:
344  break;
345 
346  default:
347  ERROR("Particle::getINCLMass: Unknown particle type." << std::endl);
348  return 0.0;
349  break;
350  }
351  }
352 
354  inline virtual G4double getTableMass() const {
355  switch(theType) {
356  case Proton:
357  case Neutron:
358  case PiPlus:
359  case PiMinus:
360  case PiZero:
362  break;
363 
364  case DeltaPlusPlus:
365  case DeltaPlus:
366  case DeltaZero:
367  case DeltaMinus:
368  return theMass;
369  break;
370 
371  case Composite:
373  break;
374 
375  default:
376  ERROR("Particle::getTableMass: Unknown particle type." << std::endl);
377  return 0.0;
378  break;
379  }
380  }
381 
383  inline G4double getRealMass() const {
384  switch(theType) {
385  case Proton:
386  case Neutron:
387  case PiPlus:
388  case PiMinus:
389  case PiZero:
391  break;
392 
393  case DeltaPlusPlus:
394  case DeltaPlus:
395  case DeltaZero:
396  case DeltaMinus:
397  return theMass;
398  break;
399 
400  case Composite:
402  break;
403 
404  default:
405  ERROR("Particle::getRealMass: Unknown particle type." << std::endl);
406  return 0.0;
407  break;
408  }
409  }
410 
413 
416 
419 
431  G4double getEmissionQValueCorrection(const G4int AParent, const G4int ZParent) const {
432  const G4int ADaughter = AParent - theA;
433  const G4int ZDaughter = ZParent - theZ;
434 
435  // Note the minus sign here
436  G4double theQValue;
437  if(isCluster())
438  theQValue = -ParticleTable::getTableQValue(theA, theZ, ADaughter, ZDaughter);
439  else {
440  const G4double massTableParent = ParticleTable::getTableMass(AParent,ZParent);
441  const G4double massTableDaughter = ParticleTable::getTableMass(ADaughter,ZDaughter);
442  const G4double massTableParticle = getTableMass();
443  theQValue = massTableParent - massTableDaughter - massTableParticle;
444  }
445 
446  const G4double massINCLParent = ParticleTable::getINCLMass(AParent,ZParent);
447  const G4double massINCLDaughter = ParticleTable::getINCLMass(ADaughter,ZDaughter);
448  const G4double massINCLParticle = getINCLMass();
449 
450  // The rhs corresponds to the INCL Q-value
451  return theQValue - (massINCLParent-massINCLDaughter-massINCLParticle);
452  }
453 
469  G4double getTransferQValueCorrection(const G4int AFrom, const G4int ZFrom, const G4int ATo, const G4int ZTo) const {
470  const G4int AFromDaughter = AFrom - theA;
471  const G4int ZFromDaughter = ZFrom - theZ;
472  const G4int AToDaughter = ATo + theA;
473  const G4int ZToDaughter = ZTo + theZ;
474  const G4double theQValue = ParticleTable::getTableQValue(AToDaughter,ZToDaughter,AFromDaughter,ZFromDaughter,AFrom,ZFrom);
475 
476  const G4double massINCLTo = ParticleTable::getINCLMass(ATo,ZTo);
477  const G4double massINCLToDaughter = ParticleTable::getINCLMass(AToDaughter,ZToDaughter);
478  /* Note that here we have to use the table mass in the INCL Q-value. We
479  * cannot use theMass, because at this stage the particle is probably
480  * still off-shell; and we cannot use getINCLMass(), because it leads to
481  * violations of global energy conservation.
482  */
483  const G4double massINCLParticle = getTableMass();
484 
485  // The rhs corresponds to the INCL Q-value for particle absorption
486  return theQValue - (massINCLToDaughter-massINCLTo-massINCLParticle);
487  }
488 
495  const G4double mass = std::pow(theEnergy, 2) - theMomentum.dot(theMomentum);
496  if(mass < 0.0) {
497  ERROR("E*E - p*p is negative." << std::endl);
498  return 0.0;
499  } else {
500  return std::sqrt(mass);
501  }
502  };
503 
505  inline G4double getKineticEnergy() const { return theEnergy - theMass; }
506 
508  inline G4double getPotentialEnergy() const { return thePotentialEnergy; }
509 
512 
517  {
518  return theEnergy;
519  };
520 
524  void setMass(G4double mass)
525  {
526  this->theMass = mass;
527  }
528 
532  void setEnergy(G4double energy)
533  {
534  this->theEnergy = energy;
535  };
536 
541  {
542  return theMomentum;
543  };
544 
547  {
549  };
550 
554  virtual void setMomentum(const G4INCL::ThreeVector &momentum)
555  {
556  this->theMomentum = momentum;
557  };
558 
563  {
564  return thePosition;
565  };
566 
567  virtual void setPosition(const G4INCL::ThreeVector &position)
568  {
569  this->thePosition = position;
570  };
571 
572  G4double getHelicity() { return theHelicity; };
573  void setHelicity(G4double h) { theHelicity = h; };
574 
575  void propagate(G4double step) {
576  thePosition += ((*thePropagationMomentum)*(step/(*thePropagationEnergy)));
577  };
578 
581 
584 
587 
589  G4int getNumberOfDecays() const { return nDecays; }
590 
593 
596 
605  void setOutOfWell() { outOfWell = true; }
606 
608  G4bool isOutOfWell() const { return outOfWell; }
609 
610  void setEmissionTime(G4double t) { emissionTime = t; }
611  G4double getEmissionTime() { return emissionTime; };
612 
616  }
617 
621  }
622 
625 
628 
630  G4bool isInList(ParticleList const &l) const {
631  return (std::find(l.begin(), l.end(), this)!=l.end());
632  }
633 
634  G4bool isCluster() const {
635  return (theType == Composite);
636  }
637 
639  void setFrozenMomentum(const ThreeVector &momentum) { theFrozenMomentum = momentum; }
640 
642  void setFrozenEnergy(const G4double energy) { theFrozenEnergy = energy; }
643 
646 
649 
651  ThreeVector getPropagationVelocity() const { return (*thePropagationMomentum)/(*thePropagationEnergy); }
652 
662  }
663 
673  }
674 
680  virtual void rotate(const G4double angle, const ThreeVector &axis) {
681  thePosition.rotate(angle, axis);
682  theMomentum.rotate(angle, axis);
683  theFrozenMomentum.rotate(angle, axis);
684  }
685 
686  std::string print() const {
687  std::stringstream ss;
688  ss << "Particle (ID = " << ID << ") type = ";
690  ss << std::endl
691  << " energy = " << theEnergy << std::endl
692  << " momentum = "
693  << theMomentum.print()
694  << std::endl
695  << " position = "
696  << thePosition.print()
697  << std::endl;
698  return ss.str();
699  };
700 
701  std::string dump() const {
702  std::stringstream ss;
703  ss << "(particle " << ID << " ";
705  ss << std::endl
706  << thePosition.dump()
707  << std::endl
708  << theMomentum.dump()
709  << std::endl
710  << theEnergy << ")" << std::endl;
711  return ss.str();
712  };
713 
714  long getID() const { return ID; };
715 
719  ParticleList const *getParticles() const {
720  WARN("Particle::getParticles() method was called on a Particle object" << std::endl);
721  return 0;
722  }
723 
724  protected:
738  long ID;
739 
740  private:
741  G4double theHelicity;
742  G4double emissionTime;
743  G4bool outOfWell;
744 
745  G4double theMass;
746  static long nextID;
747 
748  };
749 }
750 
751 #endif /* PARTICLE_HH_ */