| 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.17 2003/08/22 00:03:15 kurasige Exp $
29 // - 2 December 1995, G.Cosmo - first design, << 25 // GEANT4 tag $Name: particles-V05-02-02 $
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 // revised by V.Ivanchenko, 12 June 2003
>> 55 // fixed problem of massless particles
>> 56 // revised by V.Ivanchenko, 18 June 2003
>> 57 // take into account the case of virtual photons
>> 58 //
>> 59 //--------------------------------------------------------------
36 #include "G4DynamicParticle.hh" 60 #include "G4DynamicParticle.hh"
37 <<
38 #include "G4DecayProducts.hh" 61 #include "G4DecayProducts.hh"
39 #include "G4IonTable.hh" <<
40 #include "G4LorentzVector.hh" 62 #include "G4LorentzVector.hh"
41 #include "G4ParticleDefinition.hh" 63 #include "G4ParticleDefinition.hh"
42 #include "G4PrimaryParticle.hh" << 64 #include "G4ParticleTable.hh"
>> 65 #include "G4IonTable.hh"
43 66
44 G4Allocator<G4DynamicParticle>*& pDynamicParti << 67 G4Allocator<G4DynamicParticle> aDynamicParticleAllocator;
45 { << 68
46 G4ThreadLocalStatic G4Allocator<G4DynamicPar << 69 static const G4double EnergyMomentumRelationAllowance = keV;
47 return _instance; << 70
>> 71 ////////////////////
>> 72 G4DynamicParticle::G4DynamicParticle():
>> 73 theMomentumDirection(),
>> 74 theParticleDefinition(0),
>> 75 theKineticEnergy(0.0),
>> 76 theProperTime(0.0),
>> 77 thePreAssignedDecayProducts(0),
>> 78 thePreAssignedDecayTime(-1.0),
>> 79 verboseLevel(1),
>> 80 primaryParticle(0)
>> 81 {
>> 82 theDynamicalMass = 0.0;
>> 83 theDynamicalCharge= 0.0;
>> 84 theElectronOccupancy = 0;
>> 85 }
>> 86
>> 87 ////////////////////
>> 88 // -- constructors ----
>> 89 ////////////////////
>> 90 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
>> 91 const G4ThreeVector& aMomentumDirection,
>> 92 G4double aKineticEnergy):
>> 93 theMomentumDirection(aMomentumDirection),
>> 94 theParticleDefinition(aParticleDefinition),
>> 95 theKineticEnergy(aKineticEnergy),
>> 96 theProperTime(0.0),
>> 97 thePreAssignedDecayProducts(0),
>> 98 thePreAssignedDecayTime(-1.0),
>> 99 verboseLevel(1),
>> 100 primaryParticle(0)
>> 101 {
>> 102 // set dynamic charge/mass
>> 103 theDynamicalMass = aParticleDefinition->GetPDGMass();
>> 104 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
>> 105 AllocateElectronOccupancy();
>> 106 }
>> 107
>> 108 ////////////////////
>> 109 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
>> 110 const G4ThreeVector& aParticleMomentum):
>> 111 theParticleDefinition(aParticleDefinition),
>> 112 theProperTime(0.0),
>> 113 thePreAssignedDecayProducts(0),
>> 114 thePreAssignedDecayTime(-1.0),
>> 115 verboseLevel(1),
>> 116 primaryParticle(0)
>> 117 {
>> 118 // set dynamic charge/mass
>> 119 theDynamicalMass = aParticleDefinition->GetPDGMass();
>> 120 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
>> 121 AllocateElectronOccupancy();
>> 122
>> 123 // 3-dim momentum is given
>> 124 G4double pModule2 = aParticleMomentum.mag2();
>> 125 if (pModule2>0.0) {
>> 126 G4double mass = theDynamicalMass;
>> 127 SetKineticEnergy(sqrt(pModule2+mass*mass)-mass);
>> 128 G4double pModule = sqrt(pModule2);
>> 129 SetMomentumDirection(aParticleMomentum.x()/pModule,
>> 130 aParticleMomentum.y()/pModule,
>> 131 aParticleMomentum.z()/pModule);
>> 132 } else {
>> 133 SetMomentumDirection(1.0,0.0,0.0);
>> 134 SetKineticEnergy(0.0);
>> 135 }
48 } 136 }
49 137
50 static const G4double EnergyMomentumRelationAl << 138 ////////////////////
51 static const G4double EnergyMRA2 = << 139 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
52 EnergyMomentumRelationAllowance * EnergyMome << 140 const G4LorentzVector &aParticleMomentum):
53 << 141 theParticleDefinition(aParticleDefinition),
54 G4DynamicParticle::G4DynamicParticle() << 142 theProperTime(0.0),
55 : theMomentumDirection(0.0, 0.0, 1.0), thePo << 143 thePreAssignedDecayProducts(0),
56 {} << 144 thePreAssignedDecayTime(-1.0),
57 << 145 verboseLevel(1),
58 G4DynamicParticle::G4DynamicParticle(const G4P << 146 primaryParticle(0)
59 const G4T << 147 {
60 G4double << 148 // set dynamic charge/mass
61 : theMomentumDirection(aMomentumDirection), << 149 theDynamicalMass = aParticleDefinition->GetPDGMass();
62 thePolarization(0.0, 0.0, 0.0), << 150 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
63 theParticleDefinition(aParticleDefinition) << 151 AllocateElectronOccupancy();
64 theKineticEnergy(aKineticEnergy), << 152
65 theDynamicalMass(aParticleDefinition->GetP << 153 // 4-momentum vector (Lorentz vecotr) is given
66 theDynamicalCharge(aParticleDefinition->Ge << 154 G4double pModule2 = aParticleMomentum.x()*aParticleMomentum.x()
67 theDynamicalSpin(aParticleDefinition->GetP << 155 + aParticleMomentum.y()*aParticleMomentum.y()
68 theDynamicalMagneticMoment(aParticleDefini << 156 + aParticleMomentum.z()*aParticleMomentum.z();
69 {} << 157 if (pModule2>0.0) {
70 << 158 G4double pModule = sqrt(pModule2);
71 G4DynamicParticle::G4DynamicParticle(const G4P << 159 SetMomentumDirection(aParticleMomentum.x()/pModule,
72 const G4T << 160 aParticleMomentum.y()/pModule,
73 G4double << 161 aParticleMomentum.z()/pModule);
74 : theMomentumDirection(aMomentumDirection), << 162 G4double totalenergy = aParticleMomentum.t();
75 thePolarization(0.0, 0.0, 0.0), << 163 G4double mass2 = totalenergy*totalenergy - pModule2;
76 theParticleDefinition(aParticleDefinition) << 164 if(mass2 < 0.0001*MeV*MeV) {
77 theKineticEnergy(aKineticEnergy), << 165 theDynamicalMass = 0.;
78 theDynamicalMass(aParticleDefinition->GetP << 166 SetKineticEnergy(totalenergy);
79 theDynamicalCharge(aParticleDefinition->Ge << 167 } else {
80 theDynamicalSpin(aParticleDefinition->GetP << 168 theDynamicalMass = sqrt(mass2);
81 theDynamicalMagneticMoment(aParticleDefini << 169 SetKineticEnergy(totalenergy-theDynamicalMass);
82 { << 170 }
83 if (std::abs(theDynamicalMass - dynamicalMas << 171
84 if (dynamicalMass > EnergyMomentumRelation << 172 } else {
85 theDynamicalMass = dynamicalMass; << 173 SetMomentumDirection(1.0,0.0,0.0);
86 else << 174 SetKineticEnergy(0.0);
87 theDynamicalMass = 0.0; << 175 }
88 } << 176 }
89 } << 177
90 << 178 G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
91 G4DynamicParticle::G4DynamicParticle(const G4P << 179 G4double totalEnergy,
92 const G4T << 180 const G4ThreeVector &aParticleMomentum):
93 : thePolarization(0.0, 0.0, 0.0), << 181 theParticleDefinition(aParticleDefinition),
94 theParticleDefinition(aParticleDefinition) << 182 theProperTime(0.0),
95 theDynamicalMass(aParticleDefinition->GetP << 183 thePreAssignedDecayProducts(0),
96 theDynamicalCharge(aParticleDefinition->Ge << 184 thePreAssignedDecayTime(-1.0),
97 theDynamicalSpin(aParticleDefinition->GetP << 185 verboseLevel(1),
98 theDynamicalMagneticMoment(aParticleDefini << 186 primaryParticle(0)
99 { << 187 {
100 SetMomentum(aParticleMomentum); // 3-dim mo << 188 // set dynamic charge/mass
101 } << 189 theDynamicalMass = aParticleDefinition->GetPDGMass();
102 << 190 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
103 G4DynamicParticle::G4DynamicParticle(const G4P << 191 AllocateElectronOccupancy();
104 const G4L << 192
105 : thePolarization(0.0, 0.0, 0.0), << 193 // total energy and momentum direction are given
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() 194 G4double pModule2 = aParticleMomentum.mag2();
126 if (pModule2 > 0.0) { << 195 if (pModule2>0.0) {
127 G4double mass2 = totalEnergy * totalEnergy << 196 G4double pModule = sqrt(pModule2);
128 G4double PDGmass2 = (aParticleDefinition-> << 197 SetMomentumDirection(aParticleMomentum.x()/pModule,
129 SetMomentumDirection(aParticleMomentum.uni << 198 aParticleMomentum.y()/pModule,
130 if (mass2 < EnergyMRA2) { << 199 aParticleMomentum.z()/pModule);
>> 200
>> 201 G4double mass2 = totalEnergy*totalEnergy - pModule2;
>> 202 if(mass2 < 0.0001*MeV*MeV) {
131 theDynamicalMass = 0.; 203 theDynamicalMass = 0.;
132 SetKineticEnergy(totalEnergy); 204 SetKineticEnergy(totalEnergy);
>> 205 } else {
>> 206 theDynamicalMass = sqrt(mass2);
>> 207 SetKineticEnergy(totalEnergy-theDynamicalMass);
133 } 208 }
134 else { << 209 } else {
135 if (std::abs(PDGmass2 - mass2) > EnergyM << 210 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); 211 SetKineticEnergy(0.0);
147 } 212 }
148 } 213 }
149 214
150 G4DynamicParticle::G4DynamicParticle(const G4D << 215 ////////////////////
151 : theMomentumDirection(right.theMomentumDire << 216 G4DynamicParticle::G4DynamicParticle(const G4DynamicParticle &right)
152 thePolarization(right.thePolarization), << 217 {
153 theParticleDefinition(right.theParticleDef << 218 theDynamicalMass = right.theDynamicalMass;
154 // Don't copy preassignedDecayProducts << 219 theDynamicalCharge = right.theDynamicalCharge;
155 primaryParticle(right.primaryParticle), << 220 if (right.theElectronOccupancy != 0){
156 theKineticEnergy(right.theKineticEnergy), << 221 theElectronOccupancy =
157 theLogKineticEnergy(right.theLogKineticEne << 222 new G4ElectronOccupancy(*right.theElectronOccupancy);
158 theBeta(right.theBeta), << 223 } else {
159 theProperTime(right.theProperTime), << 224 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 } 225 }
>> 226
>> 227 theParticleDefinition = right.theParticleDefinition;
>> 228 theMomentumDirection = right.theMomentumDirection;
>> 229 theKineticEnergy = right.theKineticEnergy;
>> 230 thePolarization = right.thePolarization;
>> 231 verboseLevel = right.verboseLevel;
>> 232
>> 233 // proper time is set to zero
>> 234 theProperTime = 0.0;
>> 235
>> 236 // thePreAssignedDecayProducts/Time must not be copied.
>> 237 thePreAssignedDecayProducts = 0;
>> 238 thePreAssignedDecayTime = -1.0;
>> 239
>> 240 primaryParticle = right.primaryParticle;
171 } 241 }
172 242
173 G4DynamicParticle::G4DynamicParticle(G4Dynamic << 243 ////////////////////
174 : theMomentumDirection(from.theMomentumDirec << 244 // -- destructor ----
175 thePolarization(from.thePolarization), << 245 ////////////////////
176 theParticleDefinition(from.theParticleDefi << 246 G4DynamicParticle::~G4DynamicParticle() {
177 theElectronOccupancy(from.theElectronOccup << 247
178 // Don't move preassignedDecayProducts << 248 // delete thePreAssignedDecayProducts
179 primaryParticle(from.primaryParticle), << 249 if (thePreAssignedDecayProducts != 0) delete thePreAssignedDecayProducts;
180 theKineticEnergy(from.theKineticEnergy), << 250 thePreAssignedDecayProducts = 0;
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 251
204 delete theElectronOccupancy; << 252 if (theElectronOccupancy != 0) delete theElectronOccupancy;
205 theElectronOccupancy = nullptr; << 253 theElectronOccupancy =0;
206 } 254 }
207 255
208 G4DynamicParticle& G4DynamicParticle::operator << 256
>> 257 ////////////////////
>> 258 // -- operators ----
>> 259 ////////////////////
>> 260 G4DynamicParticle & G4DynamicParticle::operator=(const G4DynamicParticle &right)
209 { 261 {
210 if (this != &right) { 262 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; 263 theDynamicalMass = right.theDynamicalMass;
218 theDynamicalCharge = right.theDynamicalCha 264 theDynamicalCharge = right.theDynamicalCharge;
219 theDynamicalSpin = right.theDynamicalSpin; << 265 if (theElectronOccupancy != 0) delete theElectronOccupancy;
220 theDynamicalMagneticMoment = right.theDyna << 266 if (right.theElectronOccupancy != 0){
221 << 267 theElectronOccupancy =
222 delete theElectronOccupancy; << 268 new G4ElectronOccupancy(*right.theElectronOccupancy);
223 if (right.theElectronOccupancy != nullptr) << 269 } else {
224 theElectronOccupancy = new G4ElectronOcc << 270 theElectronOccupancy = 0;
225 } <<
226 else { <<
227 theElectronOccupancy = nullptr; <<
228 } 271 }
229 272
230 // thePreAssignedDecayProducts must not be << 273 theParticleDefinition = right.theParticleDefinition;
231 thePreAssignedDecayProducts = nullptr; << 274 theMomentumDirection = right.theMomentumDirection;
232 thePreAssignedDecayTime = -1.0; << 275 theKineticEnergy = right.theKineticEnergy;
233 << 276 thePolarization = right.thePolarization;
>> 277 theProperTime = right.theProperTime;
234 verboseLevel = right.verboseLevel; 278 verboseLevel = right.verboseLevel;
235 279
236 // Primary particle information must be pr << 280 // thePreAssignedDecayProducts must not be copied.
237 //*** primaryParticle = right.primaryParti << 281 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; 282 thePreAssignedDecayTime = -1.0;
265 283
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 } 284 }
276 return *this; 285 return *this;
277 } 286 }
278 287
279 void G4DynamicParticle::SetDefinition(const G4 << 288 ////////////////////
>> 289 void G4DynamicParticle::SetDefinition(G4ParticleDefinition * aParticleDefinition)
280 { 290 {
281 // remove preassigned decay 291 // remove preassigned decay
282 if (thePreAssignedDecayProducts != nullptr) << 292 if (thePreAssignedDecayProducts != 0) {
283 #ifdef G4VERBOSE << 293 G4cout << " G4DynamicParticle::SetDefinition()::";
284 if (verboseLevel > 0) { << 294 G4cout << "!!! Pre-assigned decay products is attached !!!! " << G4endl;
285 G4cout << " G4DynamicParticle::SetDefini << 295 DumpInfo(0);
286 << "!!! Pre-assigned decay produc << 296 G4cout << "!!! New Definition is " << aParticleDefinition->GetParticleName() << " !!! " << G4endl;
287 G4cout << "!!! New Definition is " << aP << 297 G4cout << "!!! Pre-assigned decay products will be deleted !!!! " << G4endl;
288 << G4endl; <<
289 G4cout << "!!! Pre-assigned decay produc <<
290 } <<
291 #endif <<
292 delete thePreAssignedDecayProducts; 298 delete thePreAssignedDecayProducts;
293 } 299 }
294 thePreAssignedDecayProducts = nullptr; << 300 thePreAssignedDecayProducts = 0;
295 301
296 theParticleDefinition = aParticleDefinition; 302 theParticleDefinition = aParticleDefinition;
297 << 303 // set Dynamic mass/chrge
298 // set Dynamic mass/charge << 304 theDynamicalMass = theParticleDefinition->GetPDGMass();
299 SetMass(theParticleDefinition->GetPDGMass()) <<
300 theDynamicalCharge = theParticleDefinition-> 305 theDynamicalCharge = theParticleDefinition->GetPDGCharge();
301 theDynamicalSpin = theParticleDefinition->Ge <<
302 theDynamicalMagneticMoment = theParticleDefi <<
303 306
304 // Set electron orbits 307 // Set electron orbits
305 if (theElectronOccupancy != nullptr) { << 308 if (theElectronOccupancy != 0) delete theElectronOccupancy;
306 delete theElectronOccupancy; << 309 theElectronOccupancy =0;
307 theElectronOccupancy = nullptr; << 310 AllocateElectronOccupancy();
308 } << 311
309 } 312 }
310 313
311 G4bool G4DynamicParticle::operator==(const G4D << 314 ////////////////////
>> 315 G4int G4DynamicParticle::operator==(const G4DynamicParticle &right) const
312 { 316 {
313 return (this == (G4DynamicParticle*)&right); << 317 return (this == (G4DynamicParticle *) &right);
314 } 318 }
315 319
316 G4bool G4DynamicParticle::operator!=(const G4D << 320 ////////////////////
>> 321 G4int G4DynamicParticle::operator!=(const G4DynamicParticle &right) const
317 { 322 {
318 return (this != (G4DynamicParticle*)&right); << 323 return (this != (G4DynamicParticle *) &right);
319 } 324 }
320 325
321 void G4DynamicParticle::AllocateElectronOccupa << 326
>> 327
>> 328 ////////////////////
>> 329 // -- AllocateElectronOccupancy --
>> 330 ////////////////////
>> 331 void G4DynamicParticle::AllocateElectronOccupancy()
322 { 332 {
323 if (G4IonTable::IsIon(theParticleDefinition) << 333 G4ParticleDefinition* particle = GetDefinition();
>> 334
>> 335 if (G4IonTable::IsIon(particle)) {
324 // Only ions can have ElectronOccupancy 336 // Only ions can have ElectronOccupancy
325 theElectronOccupancy = new G4ElectronOccup 337 theElectronOccupancy = new G4ElectronOccupancy();
326 } << 338
327 else { << 339 } else {
328 theElectronOccupancy = nullptr; << 340 theElectronOccupancy = 0;
>> 341
329 } 342 }
330 } 343 }
331 344
332 void G4DynamicParticle::SetMomentum(const G4Th << 345 ////////////////////
>> 346 // -- methods for setting Energy/Momentum --
>> 347 ////////////////////
>> 348 void G4DynamicParticle::SetMomentum(const G4ThreeVector &momentum)
333 { 349 {
334 G4double pModule2 = momentum.mag2(); 350 G4double pModule2 = momentum.mag2();
335 if (pModule2 > 0.0) { << 351 if (pModule2>0.0) {
336 const G4double mass = theDynamicalMass; << 352 G4double mass = theDynamicalMass;
337 SetMomentumDirection(momentum.unit()); << 353 G4double pModule = sqrt(pModule2);
338 SetKineticEnergy(pModule2 / (std::sqrt(pMo << 354 SetMomentumDirection(momentum.x()/pModule,
339 } << 355 momentum.y()/pModule,
340 else { << 356 momentum.z()/pModule);
341 SetMomentumDirection(1.0, 0.0, 0.0); << 357 SetKineticEnergy(sqrt(pModule2 + mass*mass)-mass);
>> 358 } else {
>> 359 SetMomentumDirection(1.0,0.0,0.0);
342 SetKineticEnergy(0.0); 360 SetKineticEnergy(0.0);
343 } 361 }
344 } 362 }
345 363
346 void G4DynamicParticle::Set4Momentum(const G4L << 364 ////////////////////
>> 365 void G4DynamicParticle::Set4Momentum(const G4LorentzVector &momentum )
347 { 366 {
348 G4double pModule2 = momentum.vect().mag2(); << 367 G4double pModule2 = momentum.x()*momentum.x()
349 if (pModule2 > 0.0) { << 368 + momentum.y()*momentum.y()
350 SetMomentumDirection(momentum.vect().unit( << 369 + momentum.z()*momentum.z();
351 const G4double totalenergy = momentum.t(); << 370 if (pModule2>0.0) {
352 const G4double mass2 = totalenergy * total << 371 G4double pModule = sqrt(pModule2);
353 const G4double PDGmass2 = << 372 SetMomentumDirection(momentum.x()/pModule,
354 (theParticleDefinition->GetPDGMass()) * << 373 momentum.y()/pModule,
355 if (mass2 < EnergyMRA2) { << 374 momentum.z()/pModule);
>> 375 G4double totalenergy = momentum.t();
>> 376 G4double mass2 = totalenergy*totalenergy - pModule2;
>> 377 if(mass2 < 0.0001*MeV*MeV) {
356 theDynamicalMass = 0.; 378 theDynamicalMass = 0.;
>> 379 SetKineticEnergy(totalenergy);
>> 380 } else {
>> 381 theDynamicalMass = sqrt(mass2);
>> 382 SetKineticEnergy(totalenergy-theDynamicalMass);
357 } 383 }
358 else if (std::abs(PDGmass2 - mass2) > Ener << 384
359 theDynamicalMass = std::sqrt(mass2); << 385 } else {
360 } << 386 SetMomentumDirection(1.0,0.0,0.0);
361 SetKineticEnergy(totalenergy - theDynamica <<
362 } <<
363 else { <<
364 SetMomentumDirection(1.0, 0.0, 0.0); <<
365 SetKineticEnergy(0.0); 387 SetKineticEnergy(0.0);
366 } 388 }
367 } 389 }
368 390
369 #ifdef G4VERBOSE << 391
>> 392 ////////////////////
>> 393 // --- Dump Information --
>> 394 ////////////////////
370 void G4DynamicParticle::DumpInfo(G4int mode) c 395 void G4DynamicParticle::DumpInfo(G4int mode) const
371 { 396 {
372 if (theParticleDefinition == nullptr) { << 397 if (theParticleDefinition == 0) {
373 G4cout << " G4DynamicParticle::DumpInfo() << 398 G4cout << " G4DynamicParticle::DumpInfo():: !!!Particle type not defined !!!! " << G4endl;
374 } << 399 } else {
375 else { <<
376 G4cout << " Particle type - " << thePartic 400 G4cout << " Particle type - " << theParticleDefinition->GetParticleName() << G4endl
377 << " mass: " << GetMass() << 401 << " mass: " << GetMass()/GeV << "[GeV]" <<G4endl
378 << " charge: " << GetCharge( << 402 << " charge: " << GetCharge()/eplus << "[e]" <<G4endl
379 << " Direction x: " << GetMomentu << 403 << " Direction x: " << GetMomentumDirection().x() << ", y: "
380 << ", y: " << GetMomentumDirection( << 404 << GetMomentumDirection().y() << ", z: "
381 << G4endl << " Total Momentum = " << 405 << GetMomentumDirection().z() << G4endl
382 << G4endl << " Momentum: " << Get << 406 << " Total Momentum = " << GetTotalMomentum() /GeV << "[GeV]" << G4endl
383 << ", y: " << GetMomentum().y() / C << 407 << " Momentum: " << GetMomentum().x() /GeV << "[GeV]" << ", y: "
384 << ", z: " << GetMomentum().z() / C << 408 << GetMomentum().y() /GeV << "[GeV]" << ", z: "
385 << " Total Energy = " << GetTot << 409 << GetMomentum().z() /GeV << "[GeV]" << G4endl
386 << " Kinetic Energy = " << GetKin << 410 << " Total Energy = " << GetTotalEnergy()/GeV << "[GeV]" << G4endl
387 << " MagneticMoment [MeV/T]: " << << 411 << " Kinetic Energy = " << GetKineticEnergy() /GeV << "[GeV]" << G4endl
388 << G4endl << " ProperTime = " << 412 << " ProperTime = " << GetProperTime() /ns << "[ns]" << G4endl;
389 << 413 if (mode>0) {
390 if (mode > 0) { << 414 if( theElectronOccupancy != 0) {
391 if (theElectronOccupancy != nullptr) { << 415 theElectronOccupancy->DumpInfo();
392 theElectronOccupancy->DumpInfo(); <<
393 } 416 }
394 } 417 }
395 } 418 }
396 } 419 }
397 #else <<
398 void G4DynamicParticle::DumpInfo(G4int) const <<
399 { <<
400 return; <<
401 } <<
402 #endif <<
403 420
404 G4double G4DynamicParticle::GetElectronMass() << 421 ////////////////////////
>> 422 G4double G4DynamicParticle::GetElectronMass() const
405 { 423 {
406 return CLHEP::electron_mass_c2; << 424 static G4double electronMass = 0.0;
>> 425
>> 426 // check if electron exits and get the mass
>> 427 if (electronMass<=0.0) {
>> 428 G4ParticleDefinition* electron = G4ParticleTable::GetParticleTable()->FindParticle("e-");
>> 429 if (electron == 0) {
>> 430 G4Exception("G4DynamicParticle: G4Electron is not defined !!");
>> 431 }
>> 432 electronMass = electron->GetPDGMass();
>> 433 }
>> 434
>> 435 return electronMass;
407 } 436 }
408 437