source: trunk/source/particles/management/src/G4DynamicParticle.cc@ 1243

Last change on this file since 1243 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 16.2 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4DynamicParticle.cc,v 1.26 2009/08/17 14:52:19 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30//
31// --------------------------------------------------------------
32// GEANT 4 class implementation file
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//--------------------------------------------------------------
63
64#include "G4DynamicParticle.hh"
65#include "G4DecayProducts.hh"
66#include "G4LorentzVector.hh"
67#include "G4ParticleDefinition.hh"
68#include "G4ParticleTable.hh"
69#include "G4IonTable.hh"
70
71G4Allocator<G4DynamicParticle> aDynamicParticleAllocator;
72
73static const G4double EnergyMomentumRelationAllowance = keV;
74
75////////////////////
76G4DynamicParticle::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////////////////////
95G4DynamicParticle::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////////////////////
116G4DynamicParticle::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 }
145}
146
147////////////////////
148G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
149 const G4LorentzVector &aParticleMomentum):
150 theParticleDefinition(aParticleDefinition),
151 theProperTime(0.0),
152 thePreAssignedDecayProducts(0),
153 thePreAssignedDecayTime(-1.0),
154 verboseLevel(1),
155 primaryParticle(0),
156 thePDGcode(0)
157{
158 // set dynamic charge/mass
159 theDynamicalMass = aParticleDefinition->GetPDGMass();
160 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
161 theDynamicalMagneticMoment = aParticleDefinition->GetPDGMagneticMoment();
162 AllocateElectronOccupancy();
163
164 // 4-momentum vector (Lorentz vecotr) is given
165 G4double pModule2 = aParticleMomentum.x()*aParticleMomentum.x()
166 + aParticleMomentum.y()*aParticleMomentum.y()
167 + aParticleMomentum.z()*aParticleMomentum.z();
168 if (pModule2>0.0) {
169 G4double pModule = std::sqrt(pModule2);
170 SetMomentumDirection(aParticleMomentum.x()/pModule,
171 aParticleMomentum.y()/pModule,
172 aParticleMomentum.z()/pModule);
173 G4double totalenergy = aParticleMomentum.t();
174 G4double mass2 = totalenergy*totalenergy - pModule2;
175 if(mass2 < EnergyMomentumRelationAllowance*EnergyMomentumRelationAllowance) {
176 theDynamicalMass = 0.;
177 SetKineticEnergy(totalenergy);
178 } else {
179 theDynamicalMass = std::sqrt(mass2);
180 SetKineticEnergy(totalenergy-theDynamicalMass);
181 }
182
183 } else {
184 SetMomentumDirection(1.0,0.0,0.0);
185 SetKineticEnergy(0.0);
186 }
187}
188
189G4DynamicParticle::G4DynamicParticle(G4ParticleDefinition * aParticleDefinition,
190 G4double totalEnergy,
191 const G4ThreeVector &aParticleMomentum):
192 theParticleDefinition(aParticleDefinition),
193 theProperTime(0.0),
194 thePreAssignedDecayProducts(0),
195 thePreAssignedDecayTime(-1.0),
196 verboseLevel(1),
197 primaryParticle(0),
198 thePDGcode(0)
199{
200 // set dynamic charge/mass
201 theDynamicalMass = aParticleDefinition->GetPDGMass();
202 theDynamicalCharge = aParticleDefinition->GetPDGCharge();
203 theDynamicalMagneticMoment = aParticleDefinition->GetPDGMagneticMoment();
204 AllocateElectronOccupancy();
205
206 // total energy and momentum direction are given
207 G4double pModule2 = aParticleMomentum.mag2();
208 if (pModule2>0.0) {
209 G4double pModule = std::sqrt(pModule2);
210 SetMomentumDirection(aParticleMomentum.x()/pModule,
211 aParticleMomentum.y()/pModule,
212 aParticleMomentum.z()/pModule);
213
214 G4double mass2 = totalEnergy*totalEnergy - pModule2;
215 if(mass2 < EnergyMomentumRelationAllowance*EnergyMomentumRelationAllowance) {
216 theDynamicalMass = 0.;
217 SetKineticEnergy(totalEnergy);
218 } else {
219 theDynamicalMass = std::sqrt(mass2);
220 SetKineticEnergy(totalEnergy-theDynamicalMass);
221 }
222 } else {
223 SetMomentumDirection(1.0,0.0,0.0);
224 SetKineticEnergy(0.0);
225 }
226}
227
228////////////////////
229G4DynamicParticle::G4DynamicParticle(const G4DynamicParticle &right)
230{
231 theDynamicalMass = right.theDynamicalMass;
232 theDynamicalCharge = right.theDynamicalCharge;
233 theDynamicalMagneticMoment = right.theDynamicalMagneticMoment;
234
235 if (right.theElectronOccupancy != 0){
236 theElectronOccupancy =
237 new G4ElectronOccupancy(*right.theElectronOccupancy);
238 } else {
239 theElectronOccupancy = 0;
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;
257}
258
259////////////////////
260// -- destructor ----
261////////////////////
262G4DynamicParticle::~G4DynamicParticle() {
263
264 // delete thePreAssignedDecayProducts
265 if (thePreAssignedDecayProducts != 0) delete thePreAssignedDecayProducts;
266 thePreAssignedDecayProducts = 0;
267
268 if (theElectronOccupancy != 0) delete theElectronOccupancy;
269 theElectronOccupancy =0;
270}
271
272
273////////////////////
274// -- operators ----
275////////////////////
276G4DynamicParticle & G4DynamicParticle::operator=(const G4DynamicParticle &right)
277{
278 if (this != &right) {
279 theDynamicalMass = right.theDynamicalMass;
280 theDynamicalCharge = right.theDynamicalCharge;
281 theDynamicalMagneticMoment = right.theDynamicalMagneticMoment;
282
283 if (theElectronOccupancy != 0) delete theElectronOccupancy;
284 if (right.theElectronOccupancy != 0){
285 theElectronOccupancy =
286 new G4ElectronOccupancy(*right.theElectronOccupancy);
287 } else {
288 theElectronOccupancy = 0;
289 }
290
291 theParticleDefinition = right.theParticleDefinition;
292 theMomentumDirection = right.theMomentumDirection;
293 theKineticEnergy = right.theKineticEnergy;
294 thePolarization = right.thePolarization;
295 theProperTime = right.theProperTime;
296 verboseLevel = right.verboseLevel;
297
298 // thePreAssignedDecayProducts must not be copied.
299 thePreAssignedDecayProducts = 0;
300 thePreAssignedDecayTime = -1.0;
301
302 }
303 return *this;
304}
305
306////////////////////
307void G4DynamicParticle::SetDefinition(G4ParticleDefinition * aParticleDefinition)
308{
309 // remove preassigned decay
310 if (thePreAssignedDecayProducts != 0) {
311#ifdef G4VERBOSE
312 if (verboseLevel>0) {
313 G4cout << " G4DynamicParticle::SetDefinition()::";
314 G4cout << "!!! Pre-assigned decay products is attached !!!! " << G4endl;
315 DumpInfo(0);
316 G4cout << "!!! New Definition is " << aParticleDefinition->GetParticleName()
317 << " !!! " << G4endl;
318 G4cout << "!!! Pre-assigned decay products will be deleted !!!! " << G4endl;
319 }
320#endif
321 delete thePreAssignedDecayProducts;
322 }
323 thePreAssignedDecayProducts = 0;
324
325 theParticleDefinition = aParticleDefinition;
326 // set Dynamic mass/chrge
327 theDynamicalMass = theParticleDefinition->GetPDGMass();
328 theDynamicalCharge = theParticleDefinition->GetPDGCharge();
329 theDynamicalMagneticMoment = theParticleDefinition->GetPDGMagneticMoment();
330
331 // Set electron orbits
332 if (theElectronOccupancy != 0) delete theElectronOccupancy;
333 theElectronOccupancy =0;
334 AllocateElectronOccupancy();
335
336}
337
338////////////////////
339G4int G4DynamicParticle::operator==(const G4DynamicParticle &right) const
340{
341 return (this == (G4DynamicParticle *) &right);
342}
343
344////////////////////
345G4int G4DynamicParticle::operator!=(const G4DynamicParticle &right) const
346{
347 return (this != (G4DynamicParticle *) &right);
348}
349
350
351
352////////////////////
353// -- AllocateElectronOccupancy --
354////////////////////
355void G4DynamicParticle::AllocateElectronOccupancy()
356{
357 G4ParticleDefinition* particle = GetDefinition();
358
359 if (G4IonTable::IsIon(particle)) {
360 // Only ions can have ElectronOccupancy
361 theElectronOccupancy = new G4ElectronOccupancy();
362
363 } else {
364 theElectronOccupancy = 0;
365
366 }
367}
368
369////////////////////
370// -- methods for setting Energy/Momentum --
371////////////////////
372void G4DynamicParticle::SetMomentum(const G4ThreeVector &momentum)
373{
374 G4double pModule2 = momentum.mag2();
375 if (pModule2>0.0) {
376 G4double mass = theDynamicalMass;
377 G4double pModule = std::sqrt(pModule2);
378 SetMomentumDirection(momentum.x()/pModule,
379 momentum.y()/pModule,
380 momentum.z()/pModule);
381 SetKineticEnergy(std::sqrt(pModule2 + mass*mass)-mass);
382 } else {
383 SetMomentumDirection(1.0,0.0,0.0);
384 SetKineticEnergy(0.0);
385 }
386}
387
388////////////////////
389void G4DynamicParticle::Set4Momentum(const G4LorentzVector &momentum )
390{
391 G4double pModule2 = momentum.x()*momentum.x()
392 + momentum.y()*momentum.y()
393 + momentum.z()*momentum.z();
394 if (pModule2>0.0) {
395 G4double pModule = std::sqrt(pModule2);
396 SetMomentumDirection(momentum.x()/pModule,
397 momentum.y()/pModule,
398 momentum.z()/pModule);
399 G4double totalenergy = momentum.t();
400 G4double mass2 = totalenergy*totalenergy - pModule2;
401 if(mass2 < 0.0001*MeV*MeV) {
402 theDynamicalMass = 0.;
403 SetKineticEnergy(totalenergy);
404 } else {
405 theDynamicalMass = std::sqrt(mass2);
406 SetKineticEnergy(totalenergy-theDynamicalMass);
407 }
408
409 } else {
410 SetMomentumDirection(1.0,0.0,0.0);
411 SetKineticEnergy(0.0);
412 }
413}
414
415
416////////////////////
417// --- Dump Information --
418////////////////////
419void G4DynamicParticle::DumpInfo(G4int mode) const
420{
421#ifdef G4VERBOSE
422 if (theParticleDefinition == 0) {
423 G4cout << " G4DynamicParticle::DumpInfo():: !!!Particle type not defined !!!! " << G4endl;
424 } else {
425 G4cout << " Particle type - " << theParticleDefinition->GetParticleName() << G4endl
426 << " mass: " << GetMass()/GeV << "[GeV]" <<G4endl
427 << " charge: " << GetCharge()/eplus << "[e]" <<G4endl
428 << " Direction x: " << GetMomentumDirection().x() << ", y: "
429 << GetMomentumDirection().y() << ", z: "
430 << GetMomentumDirection().z() << G4endl
431 << " Total Momentum = " << GetTotalMomentum() /GeV << "[GeV]" << G4endl
432 << " Momentum: " << GetMomentum().x() /GeV << "[GeV]" << ", y: "
433 << GetMomentum().y() /GeV << "[GeV]" << ", z: "
434 << GetMomentum().z() /GeV << "[GeV]" << G4endl
435 << " Total Energy = " << GetTotalEnergy()/GeV << "[GeV]" << G4endl
436 << " Kinetic Energy = " << GetKineticEnergy() /GeV << "[GeV]" << G4endl
437 << " MagneticMoment [MeV/T]: " << GetMagneticMoment()/MeV*tesla << G4endl
438 << " ProperTime = " << GetProperTime() /ns << "[ns]" << G4endl;
439
440 if (mode>0) {
441 if( theElectronOccupancy != 0) {
442 theElectronOccupancy->DumpInfo();
443 }
444 }
445 }
446#endif
447}
448
449////////////////////////
450G4double G4DynamicParticle::GetElectronMass() const
451{
452 static G4double electronMass = 0.0;
453
454 // check if electron exits and get the mass
455 if (electronMass<=0.0) {
456 G4ParticleDefinition* electron = G4ParticleTable::GetParticleTable()->FindParticle("e-");
457 if (electron == 0) {
458 G4Exception("G4DynamicParticle: G4Electron is not defined !!");
459 }
460 electronMass = electron->GetPDGMass();
461 }
462
463 return electronMass;
464}
Note: See TracBrowser for help on using the repository browser.