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