Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DecayProducts.cc
Go to the documentation of this file.
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: G4DecayProducts.cc 69015 2013-04-15 09:46:48Z gcosmo $
28 //
29 //
30 // ------------------------------------------------------------
31 // GEANT 4 class implementation file
32 //
33 // History: first implementation, based on object model of
34 // 10 July 1996 H.Kurashige
35 // 21 Oct 1996 H.Kurashige
36 // 12 Dec 1997 H.Kurashige
37 // 4 Apr. 2012 H.Kurashige use std::vector
38 // ------------------------------------------------------------
39 
40 #include "G4ios.hh"
41 #include "globals.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4DecayProducts.hh"
45 
46 #include "G4LorentzVector.hh"
47 #include "G4LorentzRotation.hh"
48 
49 
51  :numberOfProducts(0),theParentParticle(0)
52 {
53  theProductVector = new G4DecayProductVector();
54 }
55 
57  :numberOfProducts(0),theParentParticle(0)
58 {
59  theParentParticle = new G4DynamicParticle(aParticle);
60  theProductVector = new G4DecayProductVector();
61 }
62 
64  :numberOfProducts(0)
65 {
66  theProductVector = new G4DecayProductVector();
67 
68  // copy parent (Deep Copy)
69  theParentParticle = new G4DynamicParticle(*right.theParentParticle);
70 
71  //copy daughters (Deep Copy)
72  for (G4int index=0; index < right.numberOfProducts; index++) {
73  G4DynamicParticle* daughter = right.theProductVector->at(index);
74  G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
75 
76  G4double properTime = daughter->GetPreAssignedDecayProperTime();
77  if(properTime>0.0)pDaughter->SetPreAssignedDecayProperTime(properTime);
78 
79  const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts();
80  if (pPreAssigned) {
81  G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
82  pDaughter->SetPreAssignedDecayProducts(pPA);
83  }
84 
85  theProductVector->push_back( pDaughter );
86  }
87  numberOfProducts = right.numberOfProducts;
88 }
89 
91 {
92  G4int index;
93 
94  if (this != &right)
95  {
96  // recreate parent
97  if (theParentParticle != 0) delete theParentParticle;
98  theParentParticle = new G4DynamicParticle(*right.theParentParticle);
99 
100  // delete G4DynamicParticle objects
101  for (index=0; index < numberOfProducts; index++) {
102  delete theProductVector->at(index);
103  }
104  theProductVector->clear();
105 
106  //copy daughters (Deep Copy)
107  for (index=0; index < right.numberOfProducts; index++) {
108  G4DynamicParticle* daughter = right.theProductVector->at(index);
109  G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
110 
111  G4double properTime = daughter->GetPreAssignedDecayProperTime();
112  if(properTime>0.0) pDaughter->SetPreAssignedDecayProperTime(properTime);
113 
114  const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts();
115  if (pPreAssigned) {
116  G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
117  pDaughter->SetPreAssignedDecayProducts(pPA);
118  }
119  theProductVector->push_back( pDaughter );
120  }
121  numberOfProducts = right.numberOfProducts;
122 
123  }
124  return *this;
125 }
126 
128 {
129  //delete parent
130  if (theParentParticle != 0) delete theParentParticle;
131 
132  // delete G4DynamicParticle object
133  for (G4int index=0; index < numberOfProducts; index++) {
134  delete theProductVector->at(index);
135  }
136  theProductVector->clear();
137  numberOfProducts = 0;
138  delete theProductVector;
139 }
140 
142 {
143  if ( numberOfProducts >0 ) {
144  numberOfProducts -= 1;
145  G4DynamicParticle* part = theProductVector->back();
146  theProductVector->pop_back();
147  return part;
148  } else {
149  return 0;
150  }
151 }
152 
154 {
155  theProductVector->push_back(aParticle);
156  numberOfProducts += 1;
157  return numberOfProducts;
158 }
159 
161 {
162  if ((numberOfProducts > anIndex) && (anIndex >=0) ) {
163  return theProductVector->at(anIndex);
164  } else {
165  return 0;
166  }
167 }
168 
170 {
171  if (theParentParticle != 0) delete theParentParticle;
172  theParentParticle = new G4DynamicParticle(aParticle);
173 }
174 
175 
176 void G4DecayProducts::Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
177 {
178  // calcurate new beta
179  G4double mass = theParentParticle->GetMass();
180  G4double totalMomentum(0);
181  if (totalEnergy > mass ) totalMomentum = std::sqrt( (totalEnergy - mass)*(totalEnergy + mass) );
182  G4double betax = momentumDirection.x()*totalMomentum/totalEnergy;
183  G4double betay = momentumDirection.y()*totalMomentum/totalEnergy;
184  G4double betaz = momentumDirection.z()*totalMomentum/totalEnergy;
185  this->Boost(betax, betay, betaz);
186 }
187 
188 void G4DecayProducts::Boost(G4double newbetax, G4double newbetay, G4double newbetaz)
189 {
190  G4double mass = theParentParticle->GetMass();
191  G4double energy = theParentParticle->GetTotalEnergy();
192  G4double momentum = 0.0;
193 
194  G4ThreeVector direction(0.0,0.0,1.0);
195  G4LorentzVector p4;
196 
197  if (energy - mass > DBL_MIN) {
198  // calcurate beta of initial state
199  momentum = theParentParticle->GetTotalMomentum();
200  direction = theParentParticle->GetMomentumDirection();
201  G4double betax = -1.0*direction.x()*momentum/energy;
202  G4double betay = -1.0*direction.y()*momentum/energy;
203  G4double betaz = -1.0*direction.z()*momentum/energy;
204 
205  for (G4int index=0; index < numberOfProducts; index++) {
206  // make G4LorentzVector for secondaries
207  p4 = (theProductVector->at(index))->Get4Momentum();
208 
209  // boost secondaries to theParentParticle's rest frame
210  p4.boost(betax, betay, betaz);
211 
212  // boost secondaries to new frame
213  p4.boost(newbetax, newbetay, newbetaz);
214 
215  // change energy/momentum
216  (theProductVector->at(index))->Set4Momentum(p4);
217  }
218  } else {
219  for (G4int index=0; index < numberOfProducts; index++) {
220  // make G4LorentzVector for secondaries
221  p4 = (theProductVector->at(index))->Get4Momentum();
222 
223  // boost secondaries to new frame
224  p4.boost(newbetax, newbetay, newbetaz);
225 
226  // change energy/momentum
227  (theProductVector->at(index))->Set4Momentum(p4);
228  }
229  }
230  // make G4LorentzVector for parent in its rest frame
231  mass = theParentParticle->GetMass();
232  G4LorentzVector parent4( 0.0, 0.0, 0.0, mass);
233 
234  // boost parent to new frame
235  parent4.boost(newbetax, newbetay, newbetaz);
236 
237  // change energy/momentum
238  theParentParticle->Set4Momentum(parent4);
239 }
240 
242 {
243  G4bool returnValue = true;
244  // check parent
245  // energy/momentum
246  G4double parent_energy = theParentParticle->GetTotalEnergy();
247  G4ThreeVector direction = theParentParticle->GetMomentumDirection();
248  G4ThreeVector parent_momentum = direction*(theParentParticle->GetTotalMomentum());
249  // check momentum dirction is a unit vector
250  if ( (parent_momentum.mag() >0.0) && (std::fabs(direction.mag()-1.0) >1.0e-6 ) ) {
251 #ifdef G4VERBOSE
252  G4cout << "G4DecayProducts::IsChecked():: "
253  << " Momentum Direction Vector of Parent is not normalized "
254  << " (=" << direction.mag() << ")" << G4endl;
255 #endif
256  returnValue = false;
257  parent_momentum = parent_momentum * (1./direction.mag());
258  }
259 
260  //daughters
261  G4double mass, energy;
262  G4ThreeVector momentum;
263  G4double total_energy = parent_energy;
264  G4ThreeVector total_momentum = parent_momentum;
265  for (G4int index=0; index < numberOfProducts; index++)
266  {
267  G4DynamicParticle* part = theProductVector->at(index);
268  mass = part->GetMass();
269  energy = part->GetTotalEnergy();
270  direction = part->GetMomentumDirection();
271  momentum = direction*(part->GetTotalMomentum());
272  // check momentum dirction is a unit vector
273  if ( (momentum.mag()>0.0) && (std::fabs(direction.mag()-1.0) > 1.0e-6)) {
274 #ifdef G4VERBOSE
275  G4cout << "G4DecayProducts::IsChecked():: "
276  << " Momentum Direction Vector of Daughter [" << index
277  << "] is not normalized (=" << direction.mag() << ")" << G4endl;
278 #endif
279  returnValue = false;
280  momentum = momentum * (1./direction.mag());
281  }
282  // whether daughter stops or not
283  if (energy - mass < DBL_MIN ) {
284 #ifdef G4VERBOSE
285  G4cout << "G4DecayProducts::IsChecked():: "
286  << " Daughter [" << index << "] has no kinetic energy "<< G4endl;
287 #endif
288  returnValue = false;
289  }
290  total_energy -= energy;
291  total_momentum -= momentum;
292  }
293  // check energy/momentum conservation
294  if ( (std::fabs(total_energy) >1.0e-9*MeV) || (total_momentum.mag() >1.0e-9*MeV ) ){
295 #ifdef G4VERBOSE
296  G4cout << "G4DecayProducts::IsChecked():: "
297  << " Energy/Momentum is not conserved "<< G4endl;
298  G4cout << " difference between parent energy and sum of dughters' energy : "
299  << total_energy /MeV << "[MeV] " << G4endl;
300  G4cout << " difference between parent momentum and sum of dughters' momentum : "
301  << " x:" << total_momentum.getX()/MeV
302  << " y:" << total_momentum.getY()/MeV
303  << " z:" << total_momentum.getZ()/MeV
304  << G4endl;
305 #endif
306  returnValue = false;
307  }
308  return returnValue;
309 }
310 
312 {
313  G4cout << " ----- List of DecayProducts -----" << G4endl;
314  G4cout << " ------ Parent Particle ----------" << G4endl;
315  if (theParentParticle != 0) theParentParticle->DumpInfo();
316  G4cout << " ------ Daughter Particles ------" << G4endl;
317  for (G4int index=0; index < numberOfProducts; index++)
318  {
319  G4cout << " ----------" << index+1 << " -------------" << G4endl;
320  (theProductVector->at(index))-> DumpInfo();
321  }
322  G4cout << " ----- End List of DecayProducts -----" << G4endl;
323  G4cout << G4endl;
324 }
void SetPreAssignedDecayProperTime(G4double)
G4double GetTotalEnergy() const
double x() const
G4int PushProducts(G4DynamicParticle *aParticle)
void DumpInfo(G4int mode=0) const
const G4DecayProducts * GetPreAssignedDecayProducts() const
double getY() const
G4DecayProducts & operator=(const G4DecayProducts &right)
void Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
int G4int
Definition: G4Types.hh:78
double z() const
G4DynamicParticle * operator[](G4int anIndex) const
G4double GetTotalMomentum() const
double getX() const
G4GLOB_DLL std::ostream G4cout
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
void DumpInfo() const
const G4ThreeVector & GetMomentumDirection() const
HepLorentzVector & boost(double, double, double)
std::vector< G4DynamicParticle * > G4DecayProductVector
void Set4Momentum(const G4LorentzVector &momentum)
void SetParentParticle(const G4DynamicParticle &aParticle)
G4double energy(const ThreeVector &p, const G4double m)
double getZ() const
double y() const
G4DynamicParticle * PopProducts()
#define DBL_MIN
Definition: templates.hh:75
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4bool IsChecked() const
double G4double
Definition: G4Types.hh:76
void SetPreAssignedDecayProducts(G4DecayProducts *aDecayProducts)
G4double GetPreAssignedDecayProperTime() const
double mag() const