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