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

Last change on this file since 1039 was 992, checked in by garnier, 17 years ago

fichiers oublies
File size: 16.0 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.25 2007/03/11 07:17:35 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-02 $
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 < 0.0001*MeV*MeV) {
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 < 0.0001*MeV*MeV) {
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 G4cout << " G4DynamicParticle::SetDefinition()::";
312 G4cout << "!!! Pre-assigned decay products is attached !!!! " << G4endl;
313 DumpInfo(0);
314 G4cout << "!!! New Definition is " << aParticleDefinition->GetParticleName() << " !!! " << G4endl;
315 G4cout << "!!! Pre-assigned decay products will be deleted !!!! " << G4endl;
316 delete thePreAssignedDecayProducts;
317 }
318 thePreAssignedDecayProducts = 0;
319
320 theParticleDefinition = aParticleDefinition;
321 // set Dynamic mass/chrge
322 theDynamicalMass = theParticleDefinition->GetPDGMass();
323 theDynamicalCharge = theParticleDefinition->GetPDGCharge();
324 theDynamicalMagneticMoment = theParticleDefinition->GetPDGMagneticMoment();
325
326 // Set electron orbits
327 if (theElectronOccupancy != 0) delete theElectronOccupancy;
328 theElectronOccupancy =0;
329 AllocateElectronOccupancy();
330
331}
332
333////////////////////
334G4int G4DynamicParticle::operator==(const G4DynamicParticle &right) const
335{
336 return (this == (G4DynamicParticle *) &right);
337}
338
339////////////////////
340G4int G4DynamicParticle::operator!=(const G4DynamicParticle &right) const
341{
342 return (this != (G4DynamicParticle *) &right);
343}
344
345
346
347////////////////////
348// -- AllocateElectronOccupancy --
349////////////////////
350void G4DynamicParticle::AllocateElectronOccupancy()
351{
352 G4ParticleDefinition* particle = GetDefinition();
353
354 if (G4IonTable::IsIon(particle)) {
355 // Only ions can have ElectronOccupancy
356 theElectronOccupancy = new G4ElectronOccupancy();
357
358 } else {
359 theElectronOccupancy = 0;
360
361 }
362}
363
364////////////////////
365// -- methods for setting Energy/Momentum --
366////////////////////
367void G4DynamicParticle::SetMomentum(const G4ThreeVector &momentum)
368{
369 G4double pModule2 = momentum.mag2();
370 if (pModule2>0.0) {
371 G4double mass = theDynamicalMass;
372 G4double pModule = std::sqrt(pModule2);
373 SetMomentumDirection(momentum.x()/pModule,
374 momentum.y()/pModule,
375 momentum.z()/pModule);
376 SetKineticEnergy(std::sqrt(pModule2 + mass*mass)-mass);
377 } else {
378 SetMomentumDirection(1.0,0.0,0.0);
379 SetKineticEnergy(0.0);
380 }
381}
382
383////////////////////
384void G4DynamicParticle::Set4Momentum(const G4LorentzVector &momentum )
385{
386 G4double pModule2 = momentum.x()*momentum.x()
387 + momentum.y()*momentum.y()
388 + momentum.z()*momentum.z();
389 if (pModule2>0.0) {
390 G4double pModule = std::sqrt(pModule2);
391 SetMomentumDirection(momentum.x()/pModule,
392 momentum.y()/pModule,
393 momentum.z()/pModule);
394 G4double totalenergy = momentum.t();
395 G4double mass2 = totalenergy*totalenergy - pModule2;
396 if(mass2 < 0.0001*MeV*MeV) {
397 theDynamicalMass = 0.;
398 SetKineticEnergy(totalenergy);
399 } else {
400 theDynamicalMass = std::sqrt(mass2);
401 SetKineticEnergy(totalenergy-theDynamicalMass);
402 }
403
404 } else {
405 SetMomentumDirection(1.0,0.0,0.0);
406 SetKineticEnergy(0.0);
407 }
408}
409
410
411////////////////////
412// --- Dump Information --
413////////////////////
414void G4DynamicParticle::DumpInfo(G4int mode) const
415{
416 if (theParticleDefinition == 0) {
417 G4cout << " G4DynamicParticle::DumpInfo():: !!!Particle type not defined !!!! " << G4endl;
418 } else {
419 G4cout << " Particle type - " << theParticleDefinition->GetParticleName() << G4endl
420 << " mass: " << GetMass()/GeV << "[GeV]" <<G4endl
421 << " charge: " << GetCharge()/eplus << "[e]" <<G4endl
422 << " Direction x: " << GetMomentumDirection().x() << ", y: "
423 << GetMomentumDirection().y() << ", z: "
424 << GetMomentumDirection().z() << G4endl
425 << " Total Momentum = " << GetTotalMomentum() /GeV << "[GeV]" << G4endl
426 << " Momentum: " << GetMomentum().x() /GeV << "[GeV]" << ", y: "
427 << GetMomentum().y() /GeV << "[GeV]" << ", z: "
428 << GetMomentum().z() /GeV << "[GeV]" << G4endl
429 << " Total Energy = " << GetTotalEnergy()/GeV << "[GeV]" << G4endl
430 << " Kinetic Energy = " << GetKineticEnergy() /GeV << "[GeV]" << G4endl
431 << " MagneticMoment [MeV/T]: " << GetMagneticMoment()/MeV*tesla << G4endl
432 << " ProperTime = " << GetProperTime() /ns << "[ns]" << G4endl;
433
434 if (mode>0) {
435 if( theElectronOccupancy != 0) {
436 theElectronOccupancy->DumpInfo();
437 }
438 }
439 }
440}
441
442////////////////////////
443G4double G4DynamicParticle::GetElectronMass() const
444{
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;
457}
Note: See TracBrowser for help on using the repository browser.