ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DynamicParticle.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DynamicParticle.cc
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 //
27 //
28 //
29 // --------------------------------------------------------------
30 // GEANT 4 class implementation file
31 //
32 // History: first implementation, based on object model of
33 // 2nd December 1995, G.Cosmo
34 // ---------------- G4DynamicParticle ----------------
35 // first implementation by Makoto Asai, 29 January 1996
36 // revised by G.Cosmo, 29 February 1996
37 // revised by H.Kurashige 06 May 1996
38 // revised by Hisaya Kurashige, 27 July 1996
39 // modify thePreAssignedDecayProducts
40 // add void SetMomentum(G4ThreeVector &momentum)
41 // add void Set4Momentum(G4LorentzVector &momentum)
42 // add G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
43 // G4LorentzVector &p4vector)
44 // revised by Hisaya Kurashige, 19 Oct 1996
45 // add theKillProcess
46 // add ProperTime
47 // revised by Hisaya Kurashige, 26 Mar 1997
48 // modify destructor
49 // revised by Hisaya Kurashige, 05 June 1997
50 // modify DumpInfo()
51 // revised by Hisaya Kurashige, 5 June 1998
52 // remove theKillProcess
53 // revised by Hisaya Kurashige, 5 Mar 2001
54 // fixed SetDefinition()
55 // revised by V.Ivanchenko, 12 June 2003
56 // fixed problem of massless particles
57 // revised by V.Ivanchenko, 18 June 2003
58 // take into account the case of virtual photons
59 // revised by M.Kelsey 12 May 2010
60 // ensure that all constructors initialize all data members
61 //--------------------------------------------------------------
62 
63 #include "G4DynamicParticle.hh"
64 #include "G4DecayProducts.hh"
65 #include "G4LorentzVector.hh"
66 #include "G4ParticleDefinition.hh"
67 #include "G4IonTable.hh"
68 #include "G4PrimaryParticle.hh"
69 
70 G4Allocator<G4DynamicParticle>*& pDynamicParticleAllocator()
71 {
72  G4ThreadLocalStatic G4Allocator<G4DynamicParticle>* _instance = nullptr;
73  return _instance;
74 }
75 
76 static const G4double EnergyMomentumRelationAllowance = 1.0e-2*CLHEP::keV;
77 static const G4double EnergyMRA2 =
78  EnergyMomentumRelationAllowance*EnergyMomentumRelationAllowance;
79 
81 G4DynamicParticle::G4DynamicParticle():
82  theMomentumDirection(0.0,0.0,1.0),
83  thePolarization(0.0,0.0,0.0),
84  theParticleDefinition(nullptr),
85  theElectronOccupancy(nullptr),
86  thePreAssignedDecayProducts(nullptr),
87  primaryParticle(nullptr),
88  theKineticEnergy(0.0),
89  theLogKineticEnergy(DBL_MAX),
90  theProperTime(0.0),
91  theDynamicalMass(0.0),
92  theDynamicalCharge(0.0),
93  theDynamicalSpin(0.0),
94  theDynamicalMagneticMoment(0.0),
95  thePreAssignedDecayTime(-1.0),
96  verboseLevel(1),
97  thePDGcode(0)
98 {
99 }
100 
102 // -- constructors ----
104 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition * aParticleDefinition,
105  const G4ThreeVector& aMomentumDirection,
106  G4double aKineticEnergy):
107  theMomentumDirection(aMomentumDirection),
108  thePolarization(0.0,0.0,0.0),
109  theParticleDefinition(aParticleDefinition),
110  theElectronOccupancy(nullptr),
111  thePreAssignedDecayProducts(nullptr),
112  primaryParticle(nullptr),
113  theKineticEnergy(aKineticEnergy),
114  theLogKineticEnergy(DBL_MAX),
115  theProperTime(0.0),
116  theDynamicalMass(aParticleDefinition->GetPDGMass()),
117  theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
118  theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
119  theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment()),
120  thePreAssignedDecayTime(-1.0),
121  verboseLevel(1),
122  thePDGcode(0)
123 {
124 }
125 
127 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition * aParticleDefinition,
128  const G4ThreeVector& aMomentumDirection,
129  G4double aKineticEnergy,
130  const G4double dynamicalMass):
131  theMomentumDirection(aMomentumDirection),
132  thePolarization(0.0,0.0,0.0),
133  theParticleDefinition(aParticleDefinition),
134  theElectronOccupancy(nullptr),
135  thePreAssignedDecayProducts(nullptr),
136  primaryParticle(nullptr),
137  theKineticEnergy(aKineticEnergy),
138  theLogKineticEnergy(DBL_MAX),
139  theProperTime(0.0),
140  theDynamicalMass(aParticleDefinition->GetPDGMass()),
141  theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
142  theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
143  theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment()),
144  thePreAssignedDecayTime(-1.0),
145  verboseLevel(1),
146  thePDGcode(0)
147 {
148  if (std::abs(theDynamicalMass-dynamicalMass)> EnergyMomentumRelationAllowance) {
149  if (dynamicalMass>EnergyMomentumRelationAllowance) theDynamicalMass= dynamicalMass;
150  else theDynamicalMass= 0.0;
151  }
152 }
153 
155 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition * aParticleDefinition,
156  const G4ThreeVector& aParticleMomentum):
157  thePolarization(0.0,0.0,0.0),
158  theParticleDefinition(aParticleDefinition),
159  theElectronOccupancy(nullptr),
160  thePreAssignedDecayProducts(nullptr),
161  primaryParticle(nullptr),
162  theKineticEnergy(0.0),
163  theLogKineticEnergy(DBL_MAX),
164  theProperTime(0.0),
165  theDynamicalMass(aParticleDefinition->GetPDGMass()),
166  theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
167  theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
168  theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment()),
169  thePreAssignedDecayTime(-1.0),
170  verboseLevel(1),
171  thePDGcode(0)
172 {
173  SetMomentum(aParticleMomentum); // 3-dim momentum is given
174 }
175 
177 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition * aParticleDefinition,
178  const G4LorentzVector &aParticleMomentum):
179  thePolarization(0.0,0.0,0.0),
180  theParticleDefinition(aParticleDefinition),
181  theElectronOccupancy(nullptr),
182  thePreAssignedDecayProducts(nullptr),
183  primaryParticle(nullptr),
184  theKineticEnergy(0.0),
185  theLogKineticEnergy(DBL_MAX),
186  theProperTime(0.0),
187  theDynamicalMass(aParticleDefinition->GetPDGMass()),
188  theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
189  theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
190  theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment()),
191  thePreAssignedDecayTime(-1.0),
192  verboseLevel(1),
193  thePDGcode(0)
194 {
195  Set4Momentum(aParticleMomentum); // 4-momentum vector (Lorentz vector) is given
196 }
197 
199 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition * aParticleDefinition,
200  G4double totalEnergy,
201  const G4ThreeVector &aParticleMomentum):
202  thePolarization(0.0,0.0,0.0),
203  theParticleDefinition(aParticleDefinition),
204  theElectronOccupancy(nullptr),
205  thePreAssignedDecayProducts(nullptr),
206  primaryParticle(nullptr),
207  theKineticEnergy(0.0),
208  theLogKineticEnergy(DBL_MAX),
209  theProperTime(0.0),
210  theDynamicalMass(aParticleDefinition->GetPDGMass()),
211  theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
212  theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
213  theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment()),
214  thePreAssignedDecayTime(-1.0),
215  verboseLevel(1),
216  thePDGcode(0)
217 {
218  // total energy and 3-dim momentum are given
219  G4double pModule2 = aParticleMomentum.mag2();
220  if (pModule2 > 0.0) {
221  G4double mass2 = totalEnergy*totalEnergy - pModule2;
222  G4double PDGmass2 = (aParticleDefinition->GetPDGMass())*(aParticleDefinition->GetPDGMass());
223  SetMomentumDirection(aParticleMomentum.unit());
224  if (mass2 < EnergyMRA2) {
225  theDynamicalMass = 0.;
226  SetKineticEnergy(totalEnergy);
227  } else {
228  if (std::abs(PDGmass2-mass2)>EnergyMRA2){
229  theDynamicalMass = std::sqrt(mass2);
230  SetKineticEnergy(totalEnergy-theDynamicalMass);
231  } else {
232  SetKineticEnergy(totalEnergy-theDynamicalMass);
233  }
234  }
235  } else {
236  SetMomentumDirection(1.0,0.0,0.0);
237  SetKineticEnergy(0.0);
238  }
239 }
240 
242 G4DynamicParticle::G4DynamicParticle(const G4DynamicParticle &right):
243  theMomentumDirection(right.theMomentumDirection),
244  thePolarization(right.thePolarization),
245  theParticleDefinition(right.theParticleDefinition),
246  theElectronOccupancy(nullptr),
247  thePreAssignedDecayProducts(nullptr), // Do not copy preassignedDecayProducts
248  primaryParticle(right.primaryParticle),
249  theKineticEnergy(right.theKineticEnergy),
250  theLogKineticEnergy(right.theLogKineticEnergy),
251  theProperTime(right.theProperTime),
252  theDynamicalMass(right.theDynamicalMass),
253  theDynamicalCharge(right.theDynamicalCharge),
254  theDynamicalSpin(right.theDynamicalSpin),
255  theDynamicalMagneticMoment(right.theDynamicalMagneticMoment),
256  thePreAssignedDecayTime(-1.0),
257  verboseLevel(right.verboseLevel),
258  thePDGcode(right.thePDGcode)
259 {
260  if (right.theElectronOccupancy != nullptr) {
261  theElectronOccupancy =
262  new G4ElectronOccupancy(*right.theElectronOccupancy);
263  }
264 }
265 
267 // -- destructor ----
269 G4DynamicParticle::~G4DynamicParticle()
270 {
271  // delete thePreAssignedDecayProducts
272  if (thePreAssignedDecayProducts != nullptr) delete thePreAssignedDecayProducts;
273  thePreAssignedDecayProducts = nullptr;
274 
275  if (theElectronOccupancy != nullptr) delete theElectronOccupancy;
276  theElectronOccupancy =nullptr;
277 }
278 
280 // -- operators ----
282 G4DynamicParticle & G4DynamicParticle::operator=(const G4DynamicParticle &right)
283 {
284  if (this != &right) {
285  theMomentumDirection = right.theMomentumDirection;
286  theParticleDefinition = right.theParticleDefinition;
287  thePolarization = right.thePolarization;
288  theKineticEnergy = right.theKineticEnergy;
289  theProperTime = right.theProperTime;
290 
291  theDynamicalMass = right.theDynamicalMass;
292  theDynamicalCharge = right.theDynamicalCharge;
293  theDynamicalSpin = right.theDynamicalSpin;
294  theDynamicalMagneticMoment = right.theDynamicalMagneticMoment;
295 
296  if (theElectronOccupancy != nullptr) delete theElectronOccupancy;
297  if (right.theElectronOccupancy != nullptr){
298  theElectronOccupancy =
299  new G4ElectronOccupancy(*right.theElectronOccupancy);
300  } else {
301  theElectronOccupancy = nullptr;
302  }
303 
304  // thePreAssignedDecayProducts must not be copied.
305  thePreAssignedDecayProducts = nullptr;
306  thePreAssignedDecayTime = -1.0;
307 
308  verboseLevel = right.verboseLevel;
309 
310  // Primary particle information must be preserved
311  //*** primaryParticle = right.primaryParticle;
312 
313  thePDGcode = right.thePDGcode;
314  }
315  return *this;
316 }
317 
319 void G4DynamicParticle::SetDefinition(const G4ParticleDefinition * aParticleDefinition)
320 {
321  // remove preassigned decay
322  if (thePreAssignedDecayProducts != nullptr) {
323 #ifdef G4VERBOSE
324  if (verboseLevel>0) {
325  G4cout << " G4DynamicParticle::SetDefinition()::"
326  << "!!! Pre-assigned decay products is attached !!!! " << G4endl;
327  G4cout << "!!! New Definition is " << aParticleDefinition->GetParticleName()
328  << " !!! " << G4endl;
329  G4cout << "!!! Pre-assigned decay products will be deleted !!!! " << G4endl;
330  }
331 #endif
332  delete thePreAssignedDecayProducts;
333  }
334  thePreAssignedDecayProducts = nullptr;
335 
336  theParticleDefinition = aParticleDefinition;
337 
338  // set Dynamic mass/charge
339  theDynamicalMass = theParticleDefinition->GetPDGMass();
340  theDynamicalCharge = theParticleDefinition->GetPDGCharge();
341  theDynamicalSpin = theParticleDefinition->GetPDGSpin();
342  theDynamicalMagneticMoment = theParticleDefinition->GetPDGMagneticMoment();
343 
344  // Set electron orbits
345  if (theElectronOccupancy != nullptr) {
346  delete theElectronOccupancy;
347  theElectronOccupancy = nullptr;
348  }
349 }
350 
352 G4bool G4DynamicParticle::operator==(const G4DynamicParticle &right) const
353 {
354  return (this == (G4DynamicParticle *) &right);
355 }
356 
358 G4bool G4DynamicParticle::operator!=(const G4DynamicParticle &right) const
359 {
360  return (this != (G4DynamicParticle *) &right);
361 }
362 
364 // -- AllocateElectronOccupancy --
366 void G4DynamicParticle::AllocateElectronOccupancy()
367 {
368  if (G4IonTable::IsIon(theParticleDefinition)) {
369  // Only ions can have ElectronOccupancy
370  theElectronOccupancy = new G4ElectronOccupancy();
371 
372  } else {
373  theElectronOccupancy = nullptr;
374 
375  }
376 }
377 
379 // -- methods for setting Energy/Momentum --
381 void G4DynamicParticle::SetMomentum(const G4ThreeVector &momentum)
382 {
383  G4double pModule2 = momentum.mag2();
384  if (pModule2>0.0) {
385  G4double mass = theDynamicalMass;
386  SetMomentumDirection(momentum.unit());
387  SetKineticEnergy(pModule2/(std::sqrt(pModule2 + mass*mass)+mass));
388  } else {
389  SetMomentumDirection(1.0,0.0,0.0);
390  SetKineticEnergy(0.0);
391  }
392 }
393 
395 void G4DynamicParticle::Set4Momentum(const G4LorentzVector &momentum )
396 {
397  G4double pModule2 = momentum.vect().mag2();
398  if (pModule2>0.0) {
399  SetMomentumDirection(momentum.vect().unit());
400  G4double totalenergy = momentum.t();
401  G4double mass2 = totalenergy*totalenergy - pModule2;
402  G4double PDGmass2 = (theParticleDefinition->GetPDGMass())*(theParticleDefinition->GetPDGMass());
403  if (mass2 < EnergyMRA2) {
404  theDynamicalMass = 0.;
405  } else if (std::abs(PDGmass2-mass2)>EnergyMRA2){
406  theDynamicalMass = std::sqrt(mass2);
407  }
408  SetKineticEnergy(totalenergy-theDynamicalMass);
409  } else {
410  SetMomentumDirection(1.0,0.0,0.0);
411  SetKineticEnergy(0.0);
412  }
413 }
414 
416 // --- Dump Information --
418 #ifdef G4VERBOSE
419 void G4DynamicParticle::DumpInfo(G4int mode) const
420 {
421  if (theParticleDefinition == 0) {
422  G4cout << " G4DynamicParticle::DumpInfo():: !!!Particle type not defined !!!! " << G4endl;
423  } else {
424  G4cout << " Particle type - " << theParticleDefinition->GetParticleName() << G4endl
425  << " mass: " << GetMass()/CLHEP::GeV << "[GeV]" <<G4endl
426  << " charge: " << GetCharge()/CLHEP::eplus << "[e]" <<G4endl
427  << " Direction x: " << GetMomentumDirection().x() << ", y: "
428  << GetMomentumDirection().y() << ", z: "
429  << GetMomentumDirection().z() << G4endl
430  << " Total Momentum = " << GetTotalMomentum()/CLHEP::GeV << "[GeV]" << G4endl
431  << " Momentum: " << GetMomentum().x() /CLHEP::GeV << "[GeV]" << ", y: "
432  << GetMomentum().y() /CLHEP::GeV << "[GeV]" << ", z: "
433  << GetMomentum().z() /CLHEP::GeV << "[GeV]" << G4endl
434  << " Total Energy = " << GetTotalEnergy()/CLHEP::GeV << "[GeV]" << G4endl
435  << " Kinetic Energy = " << GetKineticEnergy() /CLHEP::GeV << "[GeV]" << G4endl
436  << " MagneticMoment [MeV/T]: " << GetMagneticMoment()/CLHEP::MeV*CLHEP::tesla << G4endl
437  << " ProperTime = " << GetProperTime() /CLHEP::ns << "[ns]" << G4endl;
438 
439  if (mode>0) {
440  if( theElectronOccupancy != nullptr) {
441  theElectronOccupancy->DumpInfo();
442  }
443  }
444  }
445 }
446 #else
447 void G4DynamicParticle::DumpInfo(G4int) const
448 {
449  return;
450 }
451 #endif
452 
454 G4double G4DynamicParticle::GetElectronMass() const
455 {
456  return CLHEP::electron_mass_c2;
457 }