ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DecayProducts.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DecayProducts.cc
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 //
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 class implementation file
31 //
32 // History: first implementation, based on object model of
33 // 10 July 1996 H.Kurashige
34 // 21 Oct 1996 H.Kurashige
35 // 12 Dec 1997 H.Kurashige
36 // 4 Apr. 2012 H.Kurashige use std::vector
37 // ------------------------------------------------------------
38 
39 #include "G4ios.hh"
40 #include "globals.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4DecayProducts.hh"
44 
45 #include "G4LorentzVector.hh"
46 #include "G4LorentzRotation.hh"
47 
48 
50  :numberOfProducts(0),theParentParticle(nullptr)
51 {
53 }
54 
56  :numberOfProducts(0),theParentParticle(nullptr)
57 {
58  theParentParticle = new G4DynamicParticle(aParticle);
60 }
61 
63  :numberOfProducts(0)
64 {
66 
67  // copy parent (Deep Copy)
69 
70  //copy daughters (Deep Copy)
71  for (G4int index=0; index < right.numberOfProducts; index++) {
72  G4DynamicParticle* daughter = right.theProductVector->at(index);
73  G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
74 
75  G4double properTime = daughter->GetPreAssignedDecayProperTime();
76  if(properTime>0.0)pDaughter->SetPreAssignedDecayProperTime(properTime);
77 
78  const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts();
79  if (pPreAssigned) {
80  G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
81  pDaughter->SetPreAssignedDecayProducts(pPA);
82  }
83 
84  theProductVector->push_back( pDaughter );
85  }
87 }
88 
90 {
91  G4int index;
92 
93  if (this != &right){
94  // recreate parent
95  if (theParentParticle != nullptr) delete theParentParticle;
97 
98  // delete G4DynamicParticle objects
99  for (index=0; index < numberOfProducts; index++) {
100  delete theProductVector->at(index);
101  }
102  theProductVector->clear();
103 
104  //copy daughters (Deep Copy)
105  for (index=0; index < right.numberOfProducts; index++) {
106  G4DynamicParticle* daughter = right.theProductVector->at(index);
107  G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
108 
109  G4double properTime = daughter->GetPreAssignedDecayProperTime();
110  if(properTime>0.0) pDaughter->SetPreAssignedDecayProperTime(properTime);
111 
112  const G4DecayProducts* pPreAssigned = daughter->GetPreAssignedDecayProducts();
113  if (pPreAssigned) {
114  G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
115  pDaughter->SetPreAssignedDecayProducts(pPA);
116  }
117  theProductVector->push_back( pDaughter );
118  }
119  numberOfProducts = right.numberOfProducts;
120 
121  }
122  return *this;
123 }
124 
126 {
127  //delete parent
128  if (theParentParticle != nullptr) delete theParentParticle;
129  theParentParticle = nullptr;
130 
131  // delete G4DynamicParticle object
132  for (G4int index=0; index < numberOfProducts; index++) {
133  delete theProductVector->at(index);
134  }
135  theProductVector->clear();
136  numberOfProducts = 0;
137  delete theProductVector;
138  theProductVector = nullptr;
139 }
140 
142 {
143  if ( numberOfProducts >0 ) {
144  numberOfProducts -= 1;
146  theProductVector->pop_back();
147  return part;
148  } else {
149  return nullptr;
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 nullptr;
166  }
167 }
168 
170 {
171  if (theParentParticle != nullptr) delete theParentParticle;
172  theParentParticle = new G4DynamicParticle(aParticle);
173 }
174 
175 
176 void G4DecayProducts::Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
177 {
178  // calcurate new beta
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 {
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();
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
239 }
240 
242 {
243  G4bool returnValue = true;
244  // check parent
245  // energy/momentum
246  G4double parent_energy = theParentParticle->GetTotalEnergy();
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
263  G4double total_energy = parent_energy;
264  G4ThreeVector total_momentum = parent_momentum;
265 
266  for (G4int index=0; index < numberOfProducts; 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;
316  G4cout << " ------ Daughter Particles ------" << G4endl;
317  for (G4int index=0; index < numberOfProducts; index++) {
318  G4cout << " ----------" << index+1 << " -------------" << G4endl;
319  (theProductVector->at(index))-> DumpInfo();
320  }
321  G4cout << " ----- End List of DecayProducts -----" << G4endl;
322  G4cout << G4endl;
323 }