| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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 G4Dynam
46 {
47 theParentParticle = new G4DynamicParticle(aP
48 theProductVector = new G4DecayProductVector(
49 }
50
51 G4DecayProducts::G4DecayProducts(const G4Decay
52 {
53 theProductVector = new G4DecayProductVector(
54
55 // copy parent (Deep Copy)
56 theParentParticle = new G4DynamicParticle(*r
57
58 // copy daughters (Deep Copy)
59 for (G4int index = 0; index < right.numberOf
60 G4DynamicParticle* daughter = right.thePro
61 auto pDaughter = new G4DynamicParticle(*da
62
63 G4double properTime = daughter->GetPreAssi
64 if (properTime > 0.0) pDaughter->SetPreAss
65
66 const G4DecayProducts* pPreAssigned = daug
67 if (pPreAssigned != nullptr) {
68 auto pPA = new G4DecayProducts(*pPreAssi
69 pDaughter->SetPreAssignedDecayProducts(p
70 }
71 theProductVector->push_back(pDaughter);
72 }
73 numberOfProducts = right.numberOfProducts;
74 }
75
76 G4DecayProducts& G4DecayProducts::operator=(co
77 {
78 G4int index;
79
80 if (this != &right) {
81 // recreate parent
82 delete theParentParticle;
83 theParentParticle = new G4DynamicParticle(
84
85 // delete G4DynamicParticle objects
86 for (index = 0; index < numberOfProducts;
87 delete theProductVector->at(index);
88 }
89 theProductVector->clear();
90
91 // copy daughters (Deep Copy)
92 for (index = 0; index < right.numberOfProd
93 G4DynamicParticle* daughter = right.theP
94 auto pDaughter = new G4DynamicParticle(*
95
96 G4double properTime = daughter->GetPreAs
97 if (properTime > 0.0) pDaughter->SetPreA
98
99 const G4DecayProducts* pPreAssigned = da
100 if (pPreAssigned != nullptr) {
101 auto pPA = new G4DecayProducts(*pPreAs
102 pDaughter->SetPreAssignedDecayProducts
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 < numberOfProduc
119 delete theProductVector->at(index);
120 }
121 theProductVector->clear();
122 numberOfProducts = 0;
123 delete theProductVector;
124 theProductVector = nullptr;
125 }
126
127 G4DynamicParticle* G4DecayProducts::PopProduct
128 {
129 if (numberOfProducts > 0) {
130 numberOfProducts -= 1;
131 G4DynamicParticle* part = theProductVector
132 theProductVector->pop_back();
133 return part;
134 }
135
136 return nullptr;
137 }
138
139 G4int G4DecayProducts::PushProducts(G4DynamicP
140 {
141 theProductVector->push_back(aParticle);
142 numberOfProducts += 1;
143 return numberOfProducts;
144 }
145
146 G4DynamicParticle* G4DecayProducts::operator[]
147 {
148 if ((numberOfProducts > anIndex) && (anIndex
149 return theProductVector->at(anIndex);
150 }
151
152 return nullptr;
153 }
154
155 void G4DecayProducts::SetParentParticle(const
156 {
157 delete theParentParticle;
158 theParentParticle = new G4DynamicParticle(aP
159 }
160
161 void G4DecayProducts::Boost(G4double totalEner
162 {
163 // calculate new beta
164 G4double mass = theParentParticle->GetMass()
165 G4double totalMomentum(0);
166 if (totalEnergy > mass) {
167 totalMomentum = std::sqrt((totalEnergy - m
168 }
169 G4double betax = momentumDirection.x() * tot
170 G4double betay = momentumDirection.y() * tot
171 G4double betaz = momentumDirection.z() * tot
172 Boost(betax, betay, betaz);
173 }
174
175 void G4DecayProducts::Boost(G4double newbetax,
176 {
177 G4double mass = theParentParticle->GetMass()
178 G4double energy = theParentParticle->GetTota
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->GetTotalMome
187 direction = theParentParticle->GetMomentum
188 G4double betax = -1.0 * direction.x() * mo
189 G4double betay = -1.0 * direction.y() * mo
190 G4double betaz = -1.0 * direction.z() * mo
191
192 for (G4int index = 0; index < numberOfProd
193 // make G4LorentzVector for secondaries
194 p4 = (theProductVector->at(index))->Get4
195
196 // boost secondaries to theParentParticl
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))->Set4Momen
204 }
205 }
206 else {
207 for (G4int index = 0; index < numberOfProd
208 // make G4LorentzVector for secondaries
209 p4 = (theProductVector->at(index))->Get4
210
211 // boost secondaries to new frame
212 p4.boost(newbetax, newbetay, newbetaz);
213
214 // change energy/momentum
215 (theProductVector->at(index))->Set4Momen
216 }
217 }
218
219 // make G4LorentzVector for parent in its re
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->
237 G4ThreeVector direction = theParentParticle-
238 G4ThreeVector parent_momentum = direction *
239
240 // check momentum direction is a unit vector
241 if ((parent_momentum.mag() > 0.0) && (std::f
242 #ifdef G4VERBOSE
243 G4cout << "G4DecayProducts::IsChecked()::
244 << " Momentum Direction Vector of P
245 << " (=" << direction.mag() << ")"
246 #endif
247 returnValue = false;
248 parent_momentum = parent_momentum * (1. /
249 }
250
251 // daughters
252 G4double mass, energy;
253 G4ThreeVector momentum;
254 G4double total_energy = parent_energy;
255 G4ThreeVector total_momentum = parent_moment
256
257 for (G4int index = 0; index < numberOfProduc
258 G4DynamicParticle* part = theProductVector
259 mass = part->GetMass();
260 energy = part->GetTotalEnergy();
261 direction = part->GetMomentumDirection();
262 momentum = direction * (part->GetTotalMome
263
264 // check momentum direction is a unit vect
265 if ((momentum.mag() > 0.0) && (std::fabs(d
266 #ifdef G4VERBOSE
267 G4cout << "G4DecayProducts::IsChecked():
268 << " Momentum Direction Vector of
269 << "] is not normalized (=" << d
270 #endif
271 returnValue = false;
272 momentum = momentum * (1. / direction.ma
273 }
274 // whether daughter stops or not
275 if (energy - mass < DBL_MIN) {
276 #ifdef G4VERBOSE
277 G4cout << "G4DecayProducts::IsChecked():
278 << " Daughter [" << index << "]
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)
287 #ifdef G4VERBOSE
288 G4cout << "G4DecayProducts::IsChecked()::
289 << " Energy/Momentum is not conserv
290 G4cout << " difference between parent ener
291 << "[MeV] " << G4endl;
292 G4cout << " difference between parent mome
293 << " x:" << total_momentum.getX() /
294 << " z:" << total_momentum.getZ() /
295 #endif
296 returnValue = false;
297 }
298 return returnValue;
299 }
300
301 void G4DecayProducts::DumpInfo() const
302 {
303 G4cout << " ----- List of DecayProducts ---
304 G4cout << " ------ Parent Particle ---------
305 if (theParentParticle != nullptr) theParentP
306 G4cout << " ------ Daughter Particles -----
307 for (G4int index = 0; index < numberOfProduc
308 G4cout << " ----------" << index + 1 << "
309 (theProductVector->at(index))->DumpInfo();
310 }
311 G4cout << " ----- End List of DecayProducts
312 G4cout << G4endl;
313 }
314