| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4DynamicParticle class implementation <<
27 // 23 //
28 // History: << 24 // $Id: G4DynamicParticle.cc,v 1.12 2003/03/13 09:59:39 jwellisc Exp $
29 // - 2 December 1995, G.Cosmo - first design, << 25 // GEANT4 tag $Name: geant4-05-01 $
30 // - 29 January 1996, M.Asai - first implement << 26 //
31 // - 1996 - 2007, H.Kurashige - revisions. << 27 //
32 // - 15 March 2019, M.Novak - log-kinetic en << 28 // --------------------------------------------------------------
33 // on demand if its stored << 29 // GEANT 4 class implementation file
34 //-------------------------------------------- << 30 //
35 << 31 // History: first implementation, based on object model of
>> 32 // 2nd December 1995, G.Cosmo
>> 33 // ---------------- G4DynamicParticle ----------------
>> 34 // first implementation by Makoto Asai, 29 January 1996
>> 35 // revised by G.Cosmo, 29 February 1996
>> 36 // revised by H.Kurashige 06 May 1996
>> 37 // revised by Hisaya Kurashige, 27 July 1996
>> 38 // modify thePreAssignedDecayProducts
>> 39 // add void SetMomentum(G4ThreeVector &momentum)
>> 40 // add void Set4Momentum(G4LorentzVector &momentum)
>> 41 // add G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
>> 42 // G4LorentzVector &p4vector)
>> 43 // revised by Hisaya Kurashige, 19 Oct 1996
>> 44 // add theKillProcess
>> 45 // add ProperTime
>> 46 // revised by Hisaya Kurashige, 26 Mar 1997
>> 47 // modify destructor
>> 48 // revised by Hisaya Kurashige, 05 June 1997
>> 49 // modify DumpInfo()
>> 50 // revised by Hisaya Kurashige, 5 June 1998
>> 51 // remove theKillProcess
>> 52 // revised by Hisaya Kurashige, 5 Mar 2001
>> 53 // fixed SetDefinition()
>> 54 //--------------------------------------------------------------
36 #include "G4DynamicParticle.hh" 55 #include "G4DynamicParticle.hh"
37 <<
38 #include "G4DecayProducts.hh" 56 #include "G4DecayProducts.hh"
39 #include "G4IonTable.hh" <<
40 #include "G4LorentzVector.hh" 57 #include "G4LorentzVector.hh"
41 #include "G4ParticleDefinition.hh" 58 #include "G4ParticleDefinition.hh"
42 #include "G4PrimaryParticle.hh" << 59 #include "G4ParticleTable.hh"
>> 60 #include "G4IonTable.hh"
43 61
44 G4Allocator<G4DynamicParticle>*& pDynamicParti << 62 G4Allocator<G4DynamicParticle> aDynamicParticleAllocator;
45 { << 63
46 G4ThreadLocalStatic G4Allocator<G4DynamicPar << 64 static const G4double EnergyMomentumRelationAllowance = keV;
47 return _instance; << 65
>> 66 ////////////////////
>> 67 G4DynamicParticle::G4DynamicParticle():
>> 68 theMomentumDirection(),
>> 69 theParticleDefinition(0),
>> 70 theKineticEnergy(0.0),
>> 71 theProperTime(0.0),
>> 72 thePreAssignedDecayProducts(0),
>> 73 thePreAssignedDecayTime(-1.0),
>> 74 verboseLevel(1)
>> 75 {
>> 76 theDynamicalMass = 0.0;
>> 77 theDynamicalCharge= 0.0;
>> 78 theElectronOccupancy = 0;
>> 79 }
>> 80
>> 81 ////////////////////
>> 82 // -- constructors ----
>> 83 ////////////////////
>> 84 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
>> 85 const G4ThreeVector& aMomentumDirection,
>> 86 G4double aKineticEnergy):
>> 87 theMomentumDirection(aMomentumDirection),
>> 88 theParticleDefinition(aParticleDefinition),
>> 89 theKineticEnergy(aKineticEnergy),
>> 90 theProperTime(0.0),
>> 91 thePreAssignedDecayProducts(0),
>> 92 thePreAssignedDecayTime(-1.0),
>> 93 verboseLevel(1)
>> 94 {
>> 95 // set dynamic charge/mass
>> 96 theDynamicalMass = aParticleDefinition->GetPDGMass();
>> 97 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
>> 98 AllocateElectronOccupancy();
>> 99 }
>> 100
>> 101 ////////////////////
>> 102 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
>> 103 const G4ThreeVector& aParticleMomentum):
>> 104 theParticleDefinition(aParticleDefinition),
>> 105 theProperTime(0.0),
>> 106 thePreAssignedDecayProducts(0),
>> 107 thePreAssignedDecayTime(-1.0),
>> 108 verboseLevel(1)
>> 109 {
>> 110 // set dynamic charge/mass
>> 111 theDynamicalMass = aParticleDefinition->GetPDGMass();
>> 112 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
>> 113 AllocateElectronOccupancy();
>> 114
>> 115 // 3-dim momentum is given
>> 116 G4double pModule2 = aParticleMomentum.mag2();
>> 117 if (pModule2>0.0) {
>> 118 G4double mass = theDynamicalMass;
>> 119 SetKineticEnergy(sqrt(pModule2+mass*mass)-mass);
>> 120 G4double pModule = sqrt(pModule2);
>> 121 SetMomentumDirection(aParticleMomentum.x()/pModule,
>> 122 aParticleMomentum.y()/pModule,
>> 123 aParticleMomentum.z()/pModule);
>> 124 } else {
>> 125 SetMomentumDirection(1.0,0.0,0.0);
>> 126 SetKineticEnergy(0.0);
>> 127 }
48 } 128 }
49 129
50 static const G4double EnergyMomentumRelationAl << 130 ////////////////////
51 static const G4double EnergyMRA2 = << 131 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
52 EnergyMomentumRelationAllowance * EnergyMome << 132 const G4LorentzVector &aParticleMomentum):
53 << 133 theParticleDefinition(aParticleDefinition),
54 G4DynamicParticle::G4DynamicParticle() << 134 theProperTime(0.0),
55 : theMomentumDirection(0.0, 0.0, 1.0), thePo << 135 thePreAssignedDecayProducts(0),
56 {} << 136 thePreAssignedDecayTime(-1.0),
57 << 137 verboseLevel(1)
58 G4DynamicParticle::G4DynamicParticle(const G4P << 138 {
59 const G4T << 139 // set dynamic charge/mass
60 G4double << 140 theDynamicalMass = aParticleDefinition->GetPDGMass();
61 : theMomentumDirection(aMomentumDirection), << 141 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
62 thePolarization(0.0, 0.0, 0.0), << 142 AllocateElectronOccupancy();
63 theParticleDefinition(aParticleDefinition) << 143
64 theKineticEnergy(aKineticEnergy), << 144 // 4-momentum vector (Lorentz vecotr) is given
65 theDynamicalMass(aParticleDefinition->GetP << 145 G4double pModule2 = aParticleMomentum.x()*aParticleMomentum.x()
66 theDynamicalCharge(aParticleDefinition->Ge << 146 + aParticleMomentum.y()*aParticleMomentum.y()
67 theDynamicalSpin(aParticleDefinition->GetP << 147 + aParticleMomentum.z()*aParticleMomentum.z();
68 theDynamicalMagneticMoment(aParticleDefini << 148 if (pModule2>0.0) {
69 {} << 149 G4double pModule = sqrt(pModule2);
70 << 150 SetMomentumDirection(aParticleMomentum.x()/pModule,
71 G4DynamicParticle::G4DynamicParticle(const G4P << 151 aParticleMomentum.y()/pModule,
72 const G4T << 152 aParticleMomentum.z()/pModule);
73 G4double << 153 G4double totalenergy = aParticleMomentum.t();
74 : theMomentumDirection(aMomentumDirection), << 154 if (totalenergy > pModule) {
75 thePolarization(0.0, 0.0, 0.0), << 155 G4double mass = sqrt(totalenergy*totalenergy - pModule2);
76 theParticleDefinition(aParticleDefinition) << 156 theDynamicalMass = mass;
77 theKineticEnergy(aKineticEnergy), << 157 SetKineticEnergy(totalenergy-mass);
78 theDynamicalMass(aParticleDefinition->GetP << 158 } else {
79 theDynamicalCharge(aParticleDefinition->Ge << 159 theDynamicalMass = 0.;
80 theDynamicalSpin(aParticleDefinition->GetP << 160 SetKineticEnergy(totalenergy);
81 theDynamicalMagneticMoment(aParticleDefini << 161 }
82 { << 162 } else {
83 if (std::abs(theDynamicalMass - dynamicalMas << 163 SetMomentumDirection(1.0,0.0,0.0);
84 if (dynamicalMass > EnergyMomentumRelation << 164 SetKineticEnergy(0.0);
85 theDynamicalMass = dynamicalMass; << 165 }
86 else << 166 }
87 theDynamicalMass = 0.0; << 167
88 } << 168 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
89 } << 169 G4double totalEnergy,
90 << 170 const G4ThreeVector &aParticleMomentum):
91 G4DynamicParticle::G4DynamicParticle(const G4P << 171 theParticleDefinition(aParticleDefinition),
92 const G4T << 172 theProperTime(0.0),
93 : thePolarization(0.0, 0.0, 0.0), << 173 thePreAssignedDecayProducts(0),
94 theParticleDefinition(aParticleDefinition) << 174 thePreAssignedDecayTime(-1.0),
95 theDynamicalMass(aParticleDefinition->GetP << 175 verboseLevel(1)
96 theDynamicalCharge(aParticleDefinition->Ge << 176 {
97 theDynamicalSpin(aParticleDefinition->GetP << 177 // set dynamic charge/mass
98 theDynamicalMagneticMoment(aParticleDefini << 178 theDynamicalMass = aParticleDefinition->GetPDGMass();
99 { << 179 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
100 SetMomentum(aParticleMomentum); // 3-dim mo << 180 AllocateElectronOccupancy();
101 } << 181
102 << 182 // total energy and momentum direction are given
103 G4DynamicParticle::G4DynamicParticle(const G4P <<
104 const G4L <<
105 : thePolarization(0.0, 0.0, 0.0), <<
106 theParticleDefinition(aParticleDefinition) <<
107 theDynamicalMass(aParticleDefinition->GetP <<
108 theDynamicalCharge(aParticleDefinition->Ge <<
109 theDynamicalSpin(aParticleDefinition->GetP <<
110 theDynamicalMagneticMoment(aParticleDefini <<
111 { <<
112 Set4Momentum(aParticleMomentum); // 4-momen <<
113 } <<
114 <<
115 G4DynamicParticle::G4DynamicParticle(const G4P <<
116 G4double <<
117 : thePolarization(0.0, 0.0, 0.0), <<
118 theParticleDefinition(aParticleDefinition) <<
119 theDynamicalMass(aParticleDefinition->GetP <<
120 theDynamicalCharge(aParticleDefinition->Ge <<
121 theDynamicalSpin(aParticleDefinition->GetP <<
122 theDynamicalMagneticMoment(aParticleDefini <<
123 { <<
124 // total energy and 3-dim momentum are given <<
125 G4double pModule2 = aParticleMomentum.mag2() 183 G4double pModule2 = aParticleMomentum.mag2();
126 if (pModule2 > 0.0) { << 184 if (pModule2>0.0) {
127 G4double mass2 = totalEnergy * totalEnergy << 185 G4double pModule = sqrt(pModule2);
128 G4double PDGmass2 = (aParticleDefinition-> << 186 SetMomentumDirection(aParticleMomentum.x()/pModule,
129 SetMomentumDirection(aParticleMomentum.uni << 187 aParticleMomentum.y()/pModule,
130 if (mass2 < EnergyMRA2) { << 188 aParticleMomentum.z()/pModule);
>> 189 if (totalEnergy > pModule) {
>> 190 G4double mass = sqrt(totalEnergy*totalEnergy - pModule2);
>> 191 theDynamicalMass = mass;
>> 192 SetKineticEnergy(totalEnergy-mass);
>> 193 } else {
131 theDynamicalMass = 0.; 194 theDynamicalMass = 0.;
132 SetKineticEnergy(totalEnergy); 195 SetKineticEnergy(totalEnergy);
133 } 196 }
134 else { << 197 } else {
135 if (std::abs(PDGmass2 - mass2) > EnergyM << 198 SetMomentumDirection(1.0,0.0,0.0);
136 theDynamicalMass = std::sqrt(mass2); <<
137 SetKineticEnergy(totalEnergy - theDyna <<
138 } <<
139 else { <<
140 SetKineticEnergy(totalEnergy - theDyna <<
141 } <<
142 } <<
143 } <<
144 else { <<
145 SetMomentumDirection(1.0, 0.0, 0.0); <<
146 SetKineticEnergy(0.0); 199 SetKineticEnergy(0.0);
147 } 200 }
148 } 201 }
149 202
150 G4DynamicParticle::G4DynamicParticle(const G4D << 203 ////////////////////
151 : theMomentumDirection(right.theMomentumDire << 204 G4DynamicParticle::G4DynamicParticle(const G4DynamicParticle &right)
152 thePolarization(right.thePolarization), << 205 {
153 theParticleDefinition(right.theParticleDef << 206 theDynamicalMass = right.theDynamicalMass;
154 // Don't copy preassignedDecayProducts << 207 theDynamicalCharge = right.theDynamicalCharge;
155 primaryParticle(right.primaryParticle), << 208 if (right.theElectronOccupancy != 0){
156 theKineticEnergy(right.theKineticEnergy), << 209 theElectronOccupancy =
157 theLogKineticEnergy(right.theLogKineticEne << 210 new G4ElectronOccupancy(*right.theElectronOccupancy);
158 theBeta(right.theBeta), << 211 } else {
159 theProperTime(right.theProperTime), << 212 theElectronOccupancy = 0;
160 theDynamicalMass(right.theDynamicalMass), <<
161 theDynamicalCharge(right.theDynamicalCharg <<
162 theDynamicalSpin(right.theDynamicalSpin), <<
163 theDynamicalMagneticMoment(right.theDynami <<
164 <<
165 verboseLevel(right.verboseLevel), <<
166 thePDGcode(right.thePDGcode) <<
167 { <<
168 if (right.theElectronOccupancy != nullptr) { <<
169 theElectronOccupancy = new G4ElectronOccup <<
170 } 213 }
>> 214
>> 215 theParticleDefinition = right.theParticleDefinition;
>> 216 theMomentumDirection = right.theMomentumDirection;
>> 217 theKineticEnergy = right.theKineticEnergy;
>> 218 thePolarization = right.thePolarization;
>> 219 verboseLevel = right.verboseLevel;
>> 220
>> 221 // proper time is set to zero
>> 222 theProperTime = 0.0;
>> 223
>> 224 // thePreAssignedDecayProducts/Time must not be copied.
>> 225 thePreAssignedDecayProducts = 0;
>> 226 thePreAssignedDecayTime = -1.0;
>> 227
171 } 228 }
172 229
173 G4DynamicParticle::G4DynamicParticle(G4Dynamic << 230 ////////////////////
174 : theMomentumDirection(from.theMomentumDirec << 231 // -- destructor ----
175 thePolarization(from.thePolarization), << 232 ////////////////////
176 theParticleDefinition(from.theParticleDefi << 233 G4DynamicParticle::~G4DynamicParticle() {
177 theElectronOccupancy(from.theElectronOccup <<
178 // Don't move preassignedDecayProducts <<
179 primaryParticle(from.primaryParticle), <<
180 theKineticEnergy(from.theKineticEnergy), <<
181 theLogKineticEnergy(from.theLogKineticEner <<
182 theBeta(from.theBeta), <<
183 theProperTime(from.theProperTime), <<
184 theDynamicalMass(from.theDynamicalMass), <<
185 theDynamicalCharge(from.theDynamicalCharge <<
186 theDynamicalSpin(from.theDynamicalSpin), <<
187 theDynamicalMagneticMoment(from.theDynamic <<
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 234
204 delete theElectronOccupancy; << 235 // delete thePreAssignedDecayProducts
205 theElectronOccupancy = nullptr; << 236 if (thePreAssignedDecayProducts != 0) delete thePreAssignedDecayProducts;
>> 237 thePreAssignedDecayProducts = 0;
>> 238
>> 239 if (theElectronOccupancy != 0) delete theElectronOccupancy;
>> 240 theElectronOccupancy =0;
206 } 241 }
207 242
208 G4DynamicParticle& G4DynamicParticle::operator << 243
>> 244 ////////////////////
>> 245 // -- operators ----
>> 246 ////////////////////
>> 247 G4DynamicParticle & G4DynamicParticle::operator=(const G4DynamicParticle &right)
209 { 248 {
210 if (this != &right) { 249 if (this != &right) {
211 theMomentumDirection = right.theMomentumDi <<
212 theParticleDefinition = right.theParticleD <<
213 thePolarization = right.thePolarization; <<
214 theKineticEnergy = right.theKineticEnergy; <<
215 theProperTime = right.theProperTime; <<
216 <<
217 theDynamicalMass = right.theDynamicalMass; 250 theDynamicalMass = right.theDynamicalMass;
218 theDynamicalCharge = right.theDynamicalCha 251 theDynamicalCharge = right.theDynamicalCharge;
219 theDynamicalSpin = right.theDynamicalSpin; << 252 if (right.theElectronOccupancy != 0){
220 theDynamicalMagneticMoment = right.theDyna << 253 theElectronOccupancy =
221 << 254 new G4ElectronOccupancy(*right.theElectronOccupancy);
222 delete theElectronOccupancy; << 255 } else {
223 if (right.theElectronOccupancy != nullptr) << 256 theElectronOccupancy = 0;
224 theElectronOccupancy = new G4ElectronOcc <<
225 } <<
226 else { <<
227 theElectronOccupancy = nullptr; <<
228 } 257 }
229 258
230 // thePreAssignedDecayProducts must not be << 259 theParticleDefinition = right.theParticleDefinition;
231 thePreAssignedDecayProducts = nullptr; << 260 theMomentumDirection = right.theMomentumDirection;
232 thePreAssignedDecayTime = -1.0; << 261 theKineticEnergy = right.theKineticEnergy;
233 << 262 thePolarization = right.thePolarization;
>> 263 theProperTime = right.theProperTime;
234 verboseLevel = right.verboseLevel; 264 verboseLevel = right.verboseLevel;
235 << 265
236 // Primary particle information must be pr << 266 // thePreAssignedDecayProducts must not be copied.
237 //*** primaryParticle = right.primaryParti << 267 thePreAssignedDecayProducts = 0;
238 <<
239 thePDGcode = right.thePDGcode; <<
240 } <<
241 return *this; <<
242 } <<
243 <<
244 G4DynamicParticle& G4DynamicParticle::operator <<
245 { <<
246 if (this != &from) { <<
247 theMomentumDirection = from.theMomentumDir <<
248 thePolarization = from.thePolarization; <<
249 theKineticEnergy = from.theKineticEnergy; <<
250 theProperTime = from.theProperTime; <<
251 <<
252 theDynamicalMass = from.theDynamicalMass; <<
253 theDynamicalCharge = from.theDynamicalChar <<
254 theDynamicalSpin = from.theDynamicalSpin; <<
255 theDynamicalMagneticMoment = from.theDynam <<
256 <<
257 delete theElectronOccupancy; <<
258 theElectronOccupancy = from.theElectronOcc <<
259 from.theElectronOccupancy = nullptr; <<
260 <<
261 // thePreAssignedDecayProducts must not be <<
262 thePreAssignedDecayProducts = nullptr; <<
263 from.thePreAssignedDecayProducts = nullptr <<
264 thePreAssignedDecayTime = -1.0; 268 thePreAssignedDecayTime = -1.0;
265 269
266 theParticleDefinition = from.theParticleDe <<
267 from.theParticleDefinition = nullptr; <<
268 <<
269 verboseLevel = from.verboseLevel; <<
270 <<
271 primaryParticle = from.primaryParticle; <<
272 from.primaryParticle = nullptr; <<
273 <<
274 thePDGcode = from.thePDGcode; <<
275 } 270 }
276 return *this; 271 return *this;
277 } 272 }
278 273
279 void G4DynamicParticle::SetDefinition(const G4 << 274 ////////////////////
>> 275 void G4DynamicParticle::SetDefinition(G4ParticleDefinition * aParticleDefinition)
280 { 276 {
281 // remove preassigned decay 277 // remove preassigned decay
282 if (thePreAssignedDecayProducts != nullptr) << 278 if (thePreAssignedDecayProducts != 0) {
283 #ifdef G4VERBOSE << 279 G4cout << " G4DynamicParticle::SetDefinition()::";
284 if (verboseLevel > 0) { << 280 G4cout << "!!! Pre-assigned decay products is attached !!!! " << G4endl;
285 G4cout << " G4DynamicParticle::SetDefini << 281 DumpInfo(0);
286 << "!!! Pre-assigned decay produc << 282 G4cout << "!!! New Definition is " << aParticleDefinition->GetParticleName() << " !!! " << G4endl;
287 G4cout << "!!! New Definition is " << aP << 283 G4cout << "!!! Pre-assigned decay products will be deleted !!!! " << G4endl;
288 << G4endl; <<
289 G4cout << "!!! Pre-assigned decay produc <<
290 } <<
291 #endif <<
292 delete thePreAssignedDecayProducts; 284 delete thePreAssignedDecayProducts;
293 } 285 }
294 thePreAssignedDecayProducts = nullptr; << 286 thePreAssignedDecayProducts = 0;
295 287
296 theParticleDefinition = aParticleDefinition; 288 theParticleDefinition = aParticleDefinition;
297 << 289 // set Dynamic mass/chrge
298 // set Dynamic mass/charge << 290 theDynamicalMass = theParticleDefinition->GetPDGMass();
299 SetMass(theParticleDefinition->GetPDGMass()) <<
300 theDynamicalCharge = theParticleDefinition-> 291 theDynamicalCharge = theParticleDefinition->GetPDGCharge();
301 theDynamicalSpin = theParticleDefinition->Ge <<
302 theDynamicalMagneticMoment = theParticleDefi <<
303 292
304 // Set electron orbits 293 // Set electron orbits
305 if (theElectronOccupancy != nullptr) { << 294 if (theElectronOccupancy != 0) delete theElectronOccupancy;
306 delete theElectronOccupancy; << 295 theElectronOccupancy =0;
307 theElectronOccupancy = nullptr; << 296 AllocateElectronOccupancy();
308 } << 297
309 } 298 }
310 299
311 G4bool G4DynamicParticle::operator==(const G4D << 300 ////////////////////
>> 301 G4int G4DynamicParticle::operator==(const G4DynamicParticle &right) const
312 { 302 {
313 return (this == (G4DynamicParticle*)&right); << 303 return (this == (G4DynamicParticle *) &right);
314 } 304 }
315 305
316 G4bool G4DynamicParticle::operator!=(const G4D << 306 ////////////////////
>> 307 G4int G4DynamicParticle::operator!=(const G4DynamicParticle &right) const
317 { 308 {
318 return (this != (G4DynamicParticle*)&right); << 309 return (this != (G4DynamicParticle *) &right);
319 } 310 }
320 311
321 void G4DynamicParticle::AllocateElectronOccupa << 312
>> 313
>> 314 ////////////////////
>> 315 // -- AllocateElectronOccupancy --
>> 316 ////////////////////
>> 317 void G4DynamicParticle::AllocateElectronOccupancy()
322 { 318 {
323 if (G4IonTable::IsIon(theParticleDefinition) << 319 G4ParticleDefinition* particle = GetDefinition();
324 // Only ions can have ElectronOccupancy << 320
>> 321 if (G4IonTable::IsIon(particle)) {
>> 322 // Only ions can have ElectronOccupancy
325 theElectronOccupancy = new G4ElectronOccup 323 theElectronOccupancy = new G4ElectronOccupancy();
326 } << 324
327 else { << 325 } else {
328 theElectronOccupancy = nullptr; << 326 theElectronOccupancy = 0;
>> 327
329 } 328 }
330 } 329 }
331 330
332 void G4DynamicParticle::SetMomentum(const G4Th << 331 ////////////////////
>> 332 // -- methods for setting Energy/Momentum --
>> 333 ////////////////////
>> 334 void G4DynamicParticle::SetMomentum(const G4ThreeVector &momentum)
333 { 335 {
334 G4double pModule2 = momentum.mag2(); 336 G4double pModule2 = momentum.mag2();
335 if (pModule2 > 0.0) { << 337 if (pModule2>0.0) {
336 const G4double mass = theDynamicalMass; << 338 G4double mass = theDynamicalMass;
337 SetMomentumDirection(momentum.unit()); << 339 G4double pModule = sqrt(pModule2);
338 SetKineticEnergy(pModule2 / (std::sqrt(pMo << 340 SetMomentumDirection(momentum.x()/pModule,
339 } << 341 momentum.y()/pModule,
340 else { << 342 momentum.z()/pModule);
341 SetMomentumDirection(1.0, 0.0, 0.0); << 343 SetKineticEnergy(sqrt(pModule2 + mass*mass)-mass);
>> 344 } else {
>> 345 SetMomentumDirection(1.0,0.0,0.0);
342 SetKineticEnergy(0.0); 346 SetKineticEnergy(0.0);
343 } 347 }
344 } 348 }
345 349
346 void G4DynamicParticle::Set4Momentum(const G4L << 350 ////////////////////
>> 351 void G4DynamicParticle::Set4Momentum(const G4LorentzVector &momentum )
347 { 352 {
348 G4double pModule2 = momentum.vect().mag2(); << 353 G4double pModule2 = momentum.x()*momentum.x()
349 if (pModule2 > 0.0) { << 354 + momentum.y()*momentum.y()
350 SetMomentumDirection(momentum.vect().unit( << 355 + momentum.z()*momentum.z();
351 const G4double totalenergy = momentum.t(); << 356 if (pModule2>0.0) {
352 const G4double mass2 = totalenergy * total << 357 G4double pModule = sqrt(pModule2);
353 const G4double PDGmass2 = << 358 SetMomentumDirection(momentum.x()/pModule,
354 (theParticleDefinition->GetPDGMass()) * << 359 momentum.y()/pModule,
355 if (mass2 < EnergyMRA2) { << 360 momentum.z()/pModule);
>> 361 G4double totalenergy = momentum.t();
>> 362 if (totalenergy > pModule) {
>> 363 G4double mass = sqrt(G4std::max(0., totalenergy*totalenergy - pModule2) );
>> 364 theDynamicalMass = mass;
>> 365 SetKineticEnergy(totalenergy-mass);
>> 366 } else {
356 theDynamicalMass = 0.; 367 theDynamicalMass = 0.;
>> 368 SetKineticEnergy(totalenergy);
357 } 369 }
358 else if (std::abs(PDGmass2 - mass2) > Ener << 370 } else {
359 theDynamicalMass = std::sqrt(mass2); << 371 SetMomentumDirection(1.0,0.0,0.0);
360 } <<
361 SetKineticEnergy(totalenergy - theDynamica <<
362 } <<
363 else { <<
364 SetMomentumDirection(1.0, 0.0, 0.0); <<
365 SetKineticEnergy(0.0); 372 SetKineticEnergy(0.0);
366 } 373 }
367 } 374 }
368 375
369 #ifdef G4VERBOSE << 376
>> 377 ////////////////////
>> 378 // --- Dump Information --
>> 379 ////////////////////
370 void G4DynamicParticle::DumpInfo(G4int mode) c 380 void G4DynamicParticle::DumpInfo(G4int mode) const
371 { 381 {
372 if (theParticleDefinition == nullptr) { << 382 if (theParticleDefinition == 0) {
373 G4cout << " G4DynamicParticle::DumpInfo() << 383 G4cout << " G4DynamicParticle::DumpInfo():: !!!Particle type not defined !!!! " << G4endl;
374 } << 384 } else {
375 else { <<
376 G4cout << " Particle type - " << thePartic 385 G4cout << " Particle type - " << theParticleDefinition->GetParticleName() << G4endl
377 << " mass: " << GetMass() << 386 << " mass: " << GetMass()/GeV << "[GeV]" <<G4endl
378 << " charge: " << GetCharge( << 387 << " charge: " << GetCharge()/eplus << "[e]" <<G4endl
379 << " Direction x: " << GetMomentu << 388 << " Direction x: " << GetMomentumDirection().x() << ", y: "
380 << ", y: " << GetMomentumDirection( << 389 << GetMomentumDirection().y() << ", z: "
381 << G4endl << " Total Momentum = " << 390 << GetMomentumDirection().z() << G4endl
382 << G4endl << " Momentum: " << Get << 391 << " Total Momentum = " << GetTotalMomentum() /GeV << "[GeV]" << G4endl
383 << ", y: " << GetMomentum().y() / C << 392 << " Momentum: " << GetMomentum().x() /GeV << "[GeV]" << ", y: "
384 << ", z: " << GetMomentum().z() / C << 393 << GetMomentum().y() /GeV << "[GeV]" << ", z: "
385 << " Total Energy = " << GetTot << 394 << GetMomentum().z() /GeV << "[GeV]" << G4endl
386 << " Kinetic Energy = " << GetKin << 395 << " Total Energy = " << GetTotalEnergy()/GeV << "[GeV]" << G4endl
387 << " MagneticMoment [MeV/T]: " << << 396 << " Kinetic Energy = " << GetKineticEnergy() /GeV << "[GeV]" << G4endl
388 << G4endl << " ProperTime = " << 397 << " ProperTime = " << GetProperTime() /ns << "[ns]" << G4endl;
389 << 398 if (mode>0) {
390 if (mode > 0) { << 399 if( theElectronOccupancy != 0) {
391 if (theElectronOccupancy != nullptr) { << 400 theElectronOccupancy->DumpInfo();
392 theElectronOccupancy->DumpInfo(); <<
393 } 401 }
394 } 402 }
395 } 403 }
396 } 404 }
397 #else <<
398 void G4DynamicParticle::DumpInfo(G4int) const <<
399 { <<
400 return; <<
401 } <<
402 #endif <<
403 405
404 G4double G4DynamicParticle::GetElectronMass() << 406 ////////////////////////
>> 407 G4double G4DynamicParticle::GetElectronMass() const
405 { 408 {
406 return CLHEP::electron_mass_c2; << 409 static G4double electronMass = 0.0;
>> 410
>> 411 // check if electron exits and get the mass
>> 412 if (electronMass<=0.0) {
>> 413 G4ParticleDefinition* electron = G4ParticleTable::GetParticleTable()->FindParticle("e-");
>> 414 if (electron == 0) {
>> 415 G4Exception("G4DynamicParticle: G4Electron is not defined !!");
>> 416 }
>> 417 electronMass = electron->GetPDGMass();
>> 418 }
>> 419
>> 420 return electronMass;
407 } 421 }
408 422