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

Last change on this file since 1202 was 1196, checked in by garnier, 15 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.