| 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 // G4DynamicParticle class implementation
27 //
28 // History:
29 // - 2 December 1995, G.Cosmo - first design, based on object model.
30 // - 29 January 1996, M.Asai - first implementation.
31 // - 1996 - 2007, H.Kurashige - revisions.
32 // - 15 March 2019, M.Novak - log-kinetic energy value is computed only
33 // on demand if its stored value is not up-to-date.
34 //---------------------------------------------------------------------
35
36 #include "G4DynamicParticle.hh"
37
38 #include "G4DecayProducts.hh"
39 #include "G4IonTable.hh"
40 #include "G4LorentzVector.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4PrimaryParticle.hh"
43
44 G4Allocator<G4DynamicParticle>*& pDynamicParticleAllocator()
45 {
46 G4ThreadLocalStatic G4Allocator<G4DynamicParticle>* _instance = nullptr;
47 return _instance;
48 }
49
50 static const G4double EnergyMomentumRelationAllowance = 1.0e-2 * CLHEP::keV;
51 static const G4double EnergyMRA2 =
52 EnergyMomentumRelationAllowance * EnergyMomentumRelationAllowance;
53
54 G4DynamicParticle::G4DynamicParticle()
55 : theMomentumDirection(0.0, 0.0, 1.0), thePolarization(0.0, 0.0, 0.0)
56 {}
57
58 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition* aParticleDefinition,
59 const G4ThreeVector& aMomentumDirection,
60 G4double aKineticEnergy)
61 : theMomentumDirection(aMomentumDirection),
62 thePolarization(0.0, 0.0, 0.0),
63 theParticleDefinition(aParticleDefinition),
64 theKineticEnergy(aKineticEnergy),
65 theDynamicalMass(aParticleDefinition->GetPDGMass()),
66 theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
67 theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
68 theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment())
69 {}
70
71 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition* aParticleDefinition,
72 const G4ThreeVector& aMomentumDirection,
73 G4double aKineticEnergy, const G4double dynamicalMass)
74 : theMomentumDirection(aMomentumDirection),
75 thePolarization(0.0, 0.0, 0.0),
76 theParticleDefinition(aParticleDefinition),
77 theKineticEnergy(aKineticEnergy),
78 theDynamicalMass(aParticleDefinition->GetPDGMass()),
79 theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
80 theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
81 theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment())
82 {
83 if (std::abs(theDynamicalMass - dynamicalMass) > EnergyMomentumRelationAllowance) {
84 if (dynamicalMass > EnergyMomentumRelationAllowance)
85 theDynamicalMass = dynamicalMass;
86 else
87 theDynamicalMass = 0.0;
88 }
89 }
90
91 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition* aParticleDefinition,
92 const G4ThreeVector& aParticleMomentum)
93 : thePolarization(0.0, 0.0, 0.0),
94 theParticleDefinition(aParticleDefinition),
95 theDynamicalMass(aParticleDefinition->GetPDGMass()),
96 theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
97 theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
98 theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment())
99 {
100 SetMomentum(aParticleMomentum); // 3-dim momentum is given
101 }
102
103 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition* aParticleDefinition,
104 const G4LorentzVector& aParticleMomentum)
105 : thePolarization(0.0, 0.0, 0.0),
106 theParticleDefinition(aParticleDefinition),
107 theDynamicalMass(aParticleDefinition->GetPDGMass()),
108 theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
109 theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
110 theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment())
111 {
112 Set4Momentum(aParticleMomentum); // 4-momentum vector (Lorentz vector)
113 }
114
115 G4DynamicParticle::G4DynamicParticle(const G4ParticleDefinition* aParticleDefinition,
116 G4double totalEnergy, const G4ThreeVector& aParticleMomentum)
117 : thePolarization(0.0, 0.0, 0.0),
118 theParticleDefinition(aParticleDefinition),
119 theDynamicalMass(aParticleDefinition->GetPDGMass()),
120 theDynamicalCharge(aParticleDefinition->GetPDGCharge()),
121 theDynamicalSpin(aParticleDefinition->GetPDGSpin()),
122 theDynamicalMagneticMoment(aParticleDefinition->GetPDGMagneticMoment())
123 {
124 // total energy and 3-dim momentum are given
125 G4double pModule2 = aParticleMomentum.mag2();
126 if (pModule2 > 0.0) {
127 G4double mass2 = totalEnergy * totalEnergy - pModule2;
128 G4double PDGmass2 = (aParticleDefinition->GetPDGMass()) * (aParticleDefinition->GetPDGMass());
129 SetMomentumDirection(aParticleMomentum.unit());
130 if (mass2 < EnergyMRA2) {
131 theDynamicalMass = 0.;
132 SetKineticEnergy(totalEnergy);
133 }
134 else {
135 if (std::abs(PDGmass2 - mass2) > EnergyMRA2) {
136 theDynamicalMass = std::sqrt(mass2);
137 SetKineticEnergy(totalEnergy - theDynamicalMass);
138 }
139 else {
140 SetKineticEnergy(totalEnergy - theDynamicalMass);
141 }
142 }
143 }
144 else {
145 SetMomentumDirection(1.0, 0.0, 0.0);
146 SetKineticEnergy(0.0);
147 }
148 }
149
150 G4DynamicParticle::G4DynamicParticle(const G4DynamicParticle& right)
151 : theMomentumDirection(right.theMomentumDirection),
152 thePolarization(right.thePolarization),
153 theParticleDefinition(right.theParticleDefinition),
154 // Don't copy preassignedDecayProducts
155 primaryParticle(right.primaryParticle),
156 theKineticEnergy(right.theKineticEnergy),
157 theLogKineticEnergy(right.theLogKineticEnergy),
158 theBeta(right.theBeta),
159 theProperTime(right.theProperTime),
160 theDynamicalMass(right.theDynamicalMass),
161 theDynamicalCharge(right.theDynamicalCharge),
162 theDynamicalSpin(right.theDynamicalSpin),
163 theDynamicalMagneticMoment(right.theDynamicalMagneticMoment),
164
165 verboseLevel(right.verboseLevel),
166 thePDGcode(right.thePDGcode)
167 {
168 if (right.theElectronOccupancy != nullptr) {
169 theElectronOccupancy = new G4ElectronOccupancy(*right.theElectronOccupancy);
170 }
171 }
172
173 G4DynamicParticle::G4DynamicParticle(G4DynamicParticle&& from)
174 : theMomentumDirection(from.theMomentumDirection),
175 thePolarization(from.thePolarization),
176 theParticleDefinition(from.theParticleDefinition),
177 theElectronOccupancy(from.theElectronOccupancy),
178 // Don't move preassignedDecayProducts
179 primaryParticle(from.primaryParticle),
180 theKineticEnergy(from.theKineticEnergy),
181 theLogKineticEnergy(from.theLogKineticEnergy),
182 theBeta(from.theBeta),
183 theProperTime(from.theProperTime),
184 theDynamicalMass(from.theDynamicalMass),
185 theDynamicalCharge(from.theDynamicalCharge),
186 theDynamicalSpin(from.theDynamicalSpin),
187 theDynamicalMagneticMoment(from.theDynamicalMagneticMoment),
188
189 verboseLevel(from.verboseLevel),
190 thePDGcode(from.thePDGcode)
191 {
192 // Release the data from the source object
193 from.theParticleDefinition = nullptr;
194 from.theElectronOccupancy = nullptr;
195 from.thePreAssignedDecayProducts = nullptr;
196 from.primaryParticle = nullptr;
197 }
198
199 G4DynamicParticle::~G4DynamicParticle()
200 {
201 delete thePreAssignedDecayProducts;
202 thePreAssignedDecayProducts = nullptr;
203
204 delete theElectronOccupancy;
205 theElectronOccupancy = nullptr;
206 }
207
208 G4DynamicParticle& G4DynamicParticle::operator=(const G4DynamicParticle& right)
209 {
210 if (this != &right) {
211 theMomentumDirection = right.theMomentumDirection;
212 theParticleDefinition = right.theParticleDefinition;
213 thePolarization = right.thePolarization;
214 theKineticEnergy = right.theKineticEnergy;
215 theProperTime = right.theProperTime;
216
217 theDynamicalMass = right.theDynamicalMass;
218 theDynamicalCharge = right.theDynamicalCharge;
219 theDynamicalSpin = right.theDynamicalSpin;
220 theDynamicalMagneticMoment = right.theDynamicalMagneticMoment;
221
222 delete theElectronOccupancy;
223 if (right.theElectronOccupancy != nullptr) {
224 theElectronOccupancy = new G4ElectronOccupancy(*right.theElectronOccupancy);
225 }
226 else {
227 theElectronOccupancy = nullptr;
228 }
229
230 // thePreAssignedDecayProducts must not be copied.
231 thePreAssignedDecayProducts = nullptr;
232 thePreAssignedDecayTime = -1.0;
233
234 verboseLevel = right.verboseLevel;
235
236 // Primary particle information must be preserved
237 //*** primaryParticle = right.primaryParticle;
238
239 thePDGcode = right.thePDGcode;
240 }
241 return *this;
242 }
243
244 G4DynamicParticle& G4DynamicParticle::operator=(G4DynamicParticle&& from)
245 {
246 if (this != &from) {
247 theMomentumDirection = from.theMomentumDirection;
248 thePolarization = from.thePolarization;
249 theKineticEnergy = from.theKineticEnergy;
250 theProperTime = from.theProperTime;
251
252 theDynamicalMass = from.theDynamicalMass;
253 theDynamicalCharge = from.theDynamicalCharge;
254 theDynamicalSpin = from.theDynamicalSpin;
255 theDynamicalMagneticMoment = from.theDynamicalMagneticMoment;
256
257 delete theElectronOccupancy;
258 theElectronOccupancy = from.theElectronOccupancy;
259 from.theElectronOccupancy = nullptr;
260
261 // thePreAssignedDecayProducts must not be moved.
262 thePreAssignedDecayProducts = nullptr;
263 from.thePreAssignedDecayProducts = nullptr;
264 thePreAssignedDecayTime = -1.0;
265
266 theParticleDefinition = from.theParticleDefinition;
267 from.theParticleDefinition = nullptr;
268
269 verboseLevel = from.verboseLevel;
270
271 primaryParticle = from.primaryParticle;
272 from.primaryParticle = nullptr;
273
274 thePDGcode = from.thePDGcode;
275 }
276 return *this;
277 }
278
279 void G4DynamicParticle::SetDefinition(const G4ParticleDefinition* aParticleDefinition)
280 {
281 // remove preassigned decay
282 if (thePreAssignedDecayProducts != nullptr) {
283 #ifdef G4VERBOSE
284 if (verboseLevel > 0) {
285 G4cout << " G4DynamicParticle::SetDefinition()::"
286 << "!!! Pre-assigned decay products is attached !!!! " << G4endl;
287 G4cout << "!!! New Definition is " << aParticleDefinition->GetParticleName() << " !!! "
288 << G4endl;
289 G4cout << "!!! Pre-assigned decay products will be deleted !!!! " << G4endl;
290 }
291 #endif
292 delete thePreAssignedDecayProducts;
293 }
294 thePreAssignedDecayProducts = nullptr;
295
296 theParticleDefinition = aParticleDefinition;
297
298 // set Dynamic mass/charge
299 SetMass(theParticleDefinition->GetPDGMass());
300 theDynamicalCharge = theParticleDefinition->GetPDGCharge();
301 theDynamicalSpin = theParticleDefinition->GetPDGSpin();
302 theDynamicalMagneticMoment = theParticleDefinition->GetPDGMagneticMoment();
303
304 // Set electron orbits
305 if (theElectronOccupancy != nullptr) {
306 delete theElectronOccupancy;
307 theElectronOccupancy = nullptr;
308 }
309 }
310
311 G4bool G4DynamicParticle::operator==(const G4DynamicParticle& right) const
312 {
313 return (this == (G4DynamicParticle*)&right);
314 }
315
316 G4bool G4DynamicParticle::operator!=(const G4DynamicParticle& right) const
317 {
318 return (this != (G4DynamicParticle*)&right);
319 }
320
321 void G4DynamicParticle::AllocateElectronOccupancy()
322 {
323 if (G4IonTable::IsIon(theParticleDefinition)) {
324 // Only ions can have ElectronOccupancy
325 theElectronOccupancy = new G4ElectronOccupancy();
326 }
327 else {
328 theElectronOccupancy = nullptr;
329 }
330 }
331
332 void G4DynamicParticle::SetMomentum(const G4ThreeVector& momentum)
333 {
334 G4double pModule2 = momentum.mag2();
335 if (pModule2 > 0.0) {
336 const G4double mass = theDynamicalMass;
337 SetMomentumDirection(momentum.unit());
338 SetKineticEnergy(pModule2 / (std::sqrt(pModule2 + mass * mass) + mass));
339 }
340 else {
341 SetMomentumDirection(1.0, 0.0, 0.0);
342 SetKineticEnergy(0.0);
343 }
344 }
345
346 void G4DynamicParticle::Set4Momentum(const G4LorentzVector& momentum)
347 {
348 G4double pModule2 = momentum.vect().mag2();
349 if (pModule2 > 0.0) {
350 SetMomentumDirection(momentum.vect().unit());
351 const G4double totalenergy = momentum.t();
352 const G4double mass2 = totalenergy * totalenergy - pModule2;
353 const G4double PDGmass2 =
354 (theParticleDefinition->GetPDGMass()) * (theParticleDefinition->GetPDGMass());
355 if (mass2 < EnergyMRA2) {
356 theDynamicalMass = 0.;
357 }
358 else if (std::abs(PDGmass2 - mass2) > EnergyMRA2) {
359 theDynamicalMass = std::sqrt(mass2);
360 }
361 SetKineticEnergy(totalenergy - theDynamicalMass);
362 }
363 else {
364 SetMomentumDirection(1.0, 0.0, 0.0);
365 SetKineticEnergy(0.0);
366 }
367 }
368
369 #ifdef G4VERBOSE
370 void G4DynamicParticle::DumpInfo(G4int mode) const
371 {
372 if (theParticleDefinition == nullptr) {
373 G4cout << " G4DynamicParticle::DumpInfo() - Particle type not defined !!! " << G4endl;
374 }
375 else {
376 G4cout << " Particle type - " << theParticleDefinition->GetParticleName() << G4endl
377 << " mass: " << GetMass() / CLHEP::GeV << "[GeV]" << G4endl
378 << " charge: " << GetCharge() / CLHEP::eplus << "[e]" << G4endl
379 << " Direction x: " << GetMomentumDirection().x()
380 << ", y: " << GetMomentumDirection().y() << ", z: " << GetMomentumDirection().z()
381 << G4endl << " Total Momentum = " << GetTotalMomentum() / CLHEP::GeV << "[GeV]"
382 << G4endl << " Momentum: " << GetMomentum().x() / CLHEP::GeV << "[GeV]"
383 << ", y: " << GetMomentum().y() / CLHEP::GeV << "[GeV]"
384 << ", z: " << GetMomentum().z() / CLHEP::GeV << "[GeV]" << G4endl
385 << " Total Energy = " << GetTotalEnergy() / CLHEP::GeV << "[GeV]" << G4endl
386 << " Kinetic Energy = " << GetKineticEnergy() / CLHEP::GeV << "[GeV]" << G4endl
387 << " MagneticMoment [MeV/T]: " << GetMagneticMoment() / CLHEP::MeV * CLHEP::tesla
388 << G4endl << " ProperTime = " << GetProperTime() / CLHEP::ns << "[ns]" << G4endl;
389
390 if (mode > 0) {
391 if (theElectronOccupancy != nullptr) {
392 theElectronOccupancy->DumpInfo();
393 }
394 }
395 }
396 }
397 #else
398 void G4DynamicParticle::DumpInfo(G4int) const
399 {
400 return;
401 }
402 #endif
403
404 G4double G4DynamicParticle::GetElectronMass() const
405 {
406 return CLHEP::electron_mass_c2;
407 }
408