| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // G4DecayProducts class implementation 26 // G4DecayProducts class implementation
27 // 27 //
28 // Author: H.Kurashige, 12 July 1996 28 // Author: H.Kurashige, 12 July 1996
29 // ------------------------------------------- 29 // ------------------------------------------------------------
30 30
>> 31 #include "G4ios.hh"
>> 32 #include "globals.hh"
>> 33 #include "G4PhysicalConstants.hh"
>> 34 #include "G4SystemOfUnits.hh"
31 #include "G4DecayProducts.hh" 35 #include "G4DecayProducts.hh"
32 36
33 #include "G4LorentzRotation.hh" <<
34 #include "G4LorentzVector.hh" 37 #include "G4LorentzVector.hh"
35 #include "G4PhysicalConstants.hh" << 38 #include "G4LorentzRotation.hh"
36 #include "G4SystemOfUnits.hh" <<
37 #include "G4ios.hh" <<
38 #include "globals.hh" <<
39 39
40 G4DecayProducts::G4DecayProducts() 40 G4DecayProducts::G4DecayProducts()
41 { << 41 {
42 theProductVector = new G4DecayProductVector( 42 theProductVector = new G4DecayProductVector();
43 } 43 }
44 44
45 G4DecayProducts::G4DecayProducts(const G4Dynam 45 G4DecayProducts::G4DecayProducts(const G4DynamicParticle& aParticle)
46 { 46 {
47 theParentParticle = new G4DynamicParticle(aP 47 theParentParticle = new G4DynamicParticle(aParticle);
48 theProductVector = new G4DecayProductVector( 48 theProductVector = new G4DecayProductVector();
49 } 49 }
50 50
51 G4DecayProducts::G4DecayProducts(const G4Decay << 51 G4DecayProducts::G4DecayProducts(const G4DecayProducts& right)
52 { 52 {
53 theProductVector = new G4DecayProductVector( 53 theProductVector = new G4DecayProductVector();
54 54
55 // copy parent (Deep Copy) 55 // copy parent (Deep Copy)
56 theParentParticle = new G4DynamicParticle(*r 56 theParentParticle = new G4DynamicParticle(*right.theParentParticle);
57 57
58 // copy daughters (Deep Copy) << 58 //copy daughters (Deep Copy)
59 for (G4int index = 0; index < right.numberOf << 59 for (G4int index=0; index<right.numberOfProducts; ++index)
>> 60 {
60 G4DynamicParticle* daughter = right.thePro 61 G4DynamicParticle* daughter = right.theProductVector->at(index);
61 auto pDaughter = new G4DynamicParticle(*da << 62 G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
62 63
63 G4double properTime = daughter->GetPreAssi 64 G4double properTime = daughter->GetPreAssignedDecayProperTime();
64 if (properTime > 0.0) pDaughter->SetPreAss << 65 if(properTime>0.0) pDaughter->SetPreAssignedDecayProperTime(properTime);
65 66
66 const G4DecayProducts* pPreAssigned = daug << 67 const G4DecayProducts* pPreAssigned
67 if (pPreAssigned != nullptr) { << 68 = daughter->GetPreAssignedDecayProducts();
68 auto pPA = new G4DecayProducts(*pPreAssi << 69 if (pPreAssigned != nullptr)
>> 70 {
>> 71 G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
69 pDaughter->SetPreAssignedDecayProducts(p 72 pDaughter->SetPreAssignedDecayProducts(pPA);
70 } 73 }
71 theProductVector->push_back(pDaughter); << 74 theProductVector->push_back( pDaughter );
72 } 75 }
73 numberOfProducts = right.numberOfProducts; 76 numberOfProducts = right.numberOfProducts;
74 } 77 }
75 78
76 G4DecayProducts& G4DecayProducts::operator=(co 79 G4DecayProducts& G4DecayProducts::operator=(const G4DecayProducts& right)
77 { 80 {
78 G4int index; 81 G4int index;
79 82
80 if (this != &right) { << 83 if (this != &right)
>> 84 {
81 // recreate parent 85 // recreate parent
82 delete theParentParticle; << 86 if (theParentParticle != nullptr) delete theParentParticle;
83 theParentParticle = new G4DynamicParticle( 87 theParentParticle = new G4DynamicParticle(*right.theParentParticle);
84 88
85 // delete G4DynamicParticle objects 89 // delete G4DynamicParticle objects
86 for (index = 0; index < numberOfProducts; << 90 for (index=0; index<numberOfProducts; ++index)
>> 91 {
87 delete theProductVector->at(index); 92 delete theProductVector->at(index);
88 } 93 }
89 theProductVector->clear(); 94 theProductVector->clear();
90 95
91 // copy daughters (Deep Copy) << 96 //copy daughters (Deep Copy)
92 for (index = 0; index < right.numberOfProd << 97 for (index=0; index<right.numberOfProducts; ++index)
>> 98 {
93 G4DynamicParticle* daughter = right.theP 99 G4DynamicParticle* daughter = right.theProductVector->at(index);
94 auto pDaughter = new G4DynamicParticle(* << 100 G4DynamicParticle* pDaughter = new G4DynamicParticle(*daughter);
95 101
96 G4double properTime = daughter->GetPreAs 102 G4double properTime = daughter->GetPreAssignedDecayProperTime();
97 if (properTime > 0.0) pDaughter->SetPreA << 103 if(properTime>0.0) pDaughter->SetPreAssignedDecayProperTime(properTime);
98 << 104
99 const G4DecayProducts* pPreAssigned = da << 105 const G4DecayProducts* pPreAssigned
100 if (pPreAssigned != nullptr) { << 106 = daughter->GetPreAssignedDecayProducts();
101 auto pPA = new G4DecayProducts(*pPreAs << 107 if (pPreAssigned != nullptr)
>> 108 {
>> 109 G4DecayProducts* pPA = new G4DecayProducts(*pPreAssigned);
102 pDaughter->SetPreAssignedDecayProducts 110 pDaughter->SetPreAssignedDecayProducts(pPA);
103 } 111 }
104 theProductVector->push_back(pDaughter); << 112 theProductVector->push_back( pDaughter );
105 } << 113 }
106 numberOfProducts = right.numberOfProducts; 114 numberOfProducts = right.numberOfProducts;
107 } 115 }
108 return *this; 116 return *this;
109 } 117 }
110 118
111 G4DecayProducts::~G4DecayProducts() 119 G4DecayProducts::~G4DecayProducts()
112 { 120 {
113 // delete parent << 121 //delete parent
114 delete theParentParticle; << 122 if (theParentParticle != nullptr) delete theParentParticle;
115 theParentParticle = nullptr; 123 theParentParticle = nullptr;
116 124
117 // delete G4DynamicParticle object 125 // delete G4DynamicParticle object
118 for (G4int index = 0; index < numberOfProduc << 126 for (G4int index=0; index<numberOfProducts; ++index)
>> 127 {
119 delete theProductVector->at(index); 128 delete theProductVector->at(index);
120 } 129 }
121 theProductVector->clear(); 130 theProductVector->clear();
122 numberOfProducts = 0; << 131 numberOfProducts = 0;
123 delete theProductVector; 132 delete theProductVector;
124 theProductVector = nullptr; 133 theProductVector = nullptr;
125 } 134 }
126 135
127 G4DynamicParticle* G4DecayProducts::PopProduct 136 G4DynamicParticle* G4DecayProducts::PopProducts()
128 { 137 {
129 if (numberOfProducts > 0) { << 138 if ( numberOfProducts >0 )
130 numberOfProducts -= 1; << 139 {
>> 140 numberOfProducts -= 1;
131 G4DynamicParticle* part = theProductVector 141 G4DynamicParticle* part = theProductVector->back();
132 theProductVector->pop_back(); 142 theProductVector->pop_back();
133 return part; 143 return part;
134 } 144 }
135 << 145 else
136 return nullptr; << 146 {
>> 147 return nullptr;
>> 148 }
137 } 149 }
138 150
139 G4int G4DecayProducts::PushProducts(G4DynamicP 151 G4int G4DecayProducts::PushProducts(G4DynamicParticle* aParticle)
140 { 152 {
141 theProductVector->push_back(aParticle); 153 theProductVector->push_back(aParticle);
142 numberOfProducts += 1; << 154 numberOfProducts += 1;
143 return numberOfProducts; 155 return numberOfProducts;
144 } 156 }
145 157
146 G4DynamicParticle* G4DecayProducts::operator[] 158 G4DynamicParticle* G4DecayProducts::operator[](G4int anIndex) const
147 { 159 {
148 if ((numberOfProducts > anIndex) && (anIndex << 160 if ((numberOfProducts > anIndex) && (anIndex >=0) )
149 return theProductVector->at(anIndex); << 161 {
>> 162 return theProductVector->at(anIndex);
>> 163 }
>> 164 else
>> 165 {
>> 166 return nullptr;
150 } 167 }
151 <<
152 return nullptr; <<
153 } 168 }
154 169
155 void G4DecayProducts::SetParentParticle(const 170 void G4DecayProducts::SetParentParticle(const G4DynamicParticle& aParticle)
156 { 171 {
157 delete theParentParticle; << 172 if (theParentParticle != nullptr) delete theParentParticle;
158 theParentParticle = new G4DynamicParticle(aP 173 theParentParticle = new G4DynamicParticle(aParticle);
159 } 174 }
160 175
161 void G4DecayProducts::Boost(G4double totalEner << 176
>> 177 void G4DecayProducts::Boost(G4double totalEnergy,
>> 178 const G4ThreeVector& momentumDirection)
162 { 179 {
163 // calculate new beta 180 // calculate new beta
164 G4double mass = theParentParticle->GetMass() 181 G4double mass = theParentParticle->GetMass();
165 G4double totalMomentum(0); 182 G4double totalMomentum(0);
166 if (totalEnergy > mass) { << 183 if ( totalEnergy > mass )
167 totalMomentum = std::sqrt((totalEnergy - m << 184 {
168 } << 185 totalMomentum = std::sqrt( (totalEnergy - mass)*(totalEnergy + mass) );
169 G4double betax = momentumDirection.x() * tot << 186 }
170 G4double betay = momentumDirection.y() * tot << 187 G4double betax = momentumDirection.x()*totalMomentum/totalEnergy;
171 G4double betaz = momentumDirection.z() * tot << 188 G4double betay = momentumDirection.y()*totalMomentum/totalEnergy;
>> 189 G4double betaz = momentumDirection.z()*totalMomentum/totalEnergy;
172 Boost(betax, betay, betaz); 190 Boost(betax, betay, betaz);
173 } 191 }
174 192
175 void G4DecayProducts::Boost(G4double newbetax, << 193 void G4DecayProducts::Boost(G4double newbetax,
176 { << 194 G4double newbetay,
>> 195 G4double newbetaz)
>> 196 {
177 G4double mass = theParentParticle->GetMass() 197 G4double mass = theParentParticle->GetMass();
178 G4double energy = theParentParticle->GetTota << 198 G4double energy = theParentParticle->GetTotalEnergy();
179 G4double momentum = 0.0; << 199 G4double momentum = 0.0;
180 200
181 G4ThreeVector direction(0.0, 0.0, 1.0); << 201 G4ThreeVector direction(0.0,0.0,1.0);
182 G4LorentzVector p4; 202 G4LorentzVector p4;
183 203
184 if (energy - mass > DBL_MIN) { << 204 if (energy - mass > DBL_MIN)
>> 205 {
185 // calcurate beta of initial state 206 // calcurate beta of initial state
186 momentum = theParentParticle->GetTotalMome << 207 momentum = theParentParticle->GetTotalMomentum();
187 direction = theParentParticle->GetMomentum 208 direction = theParentParticle->GetMomentumDirection();
188 G4double betax = -1.0 * direction.x() * mo << 209 G4double betax = -1.0*direction.x()*momentum/energy;
189 G4double betay = -1.0 * direction.y() * mo << 210 G4double betay = -1.0*direction.y()*momentum/energy;
190 G4double betaz = -1.0 * direction.z() * mo << 211 G4double betaz = -1.0*direction.z()*momentum/energy;
191 << 212
192 for (G4int index = 0; index < numberOfProd << 213 for (G4int index=0; index<numberOfProducts; ++index)
>> 214 {
193 // make G4LorentzVector for secondaries 215 // make G4LorentzVector for secondaries
194 p4 = (theProductVector->at(index))->Get4 216 p4 = (theProductVector->at(index))->Get4Momentum();
195 217
196 // boost secondaries to theParentParticl << 218 // boost secondaries to theParentParticle's rest frame
197 p4.boost(betax, betay, betaz); 219 p4.boost(betax, betay, betaz);
198 220
199 // boost secondaries to new frame << 221 // boost secondaries to new frame
200 p4.boost(newbetax, newbetay, newbetaz); 222 p4.boost(newbetax, newbetay, newbetaz);
201 223
202 // change energy/momentum 224 // change energy/momentum
203 (theProductVector->at(index))->Set4Momen 225 (theProductVector->at(index))->Set4Momentum(p4);
204 } 226 }
205 } 227 }
206 else { << 228 else
207 for (G4int index = 0; index < numberOfProd << 229 {
>> 230 for (G4int index=0; index<numberOfProducts; ++index)
>> 231 {
208 // make G4LorentzVector for secondaries 232 // make G4LorentzVector for secondaries
209 p4 = (theProductVector->at(index))->Get4 233 p4 = (theProductVector->at(index))->Get4Momentum();
210 234
211 // boost secondaries to new frame << 235 // boost secondaries to new frame
212 p4.boost(newbetax, newbetay, newbetaz); 236 p4.boost(newbetax, newbetay, newbetaz);
213 237
214 // change energy/momentum 238 // change energy/momentum
215 (theProductVector->at(index))->Set4Momen 239 (theProductVector->at(index))->Set4Momentum(p4);
216 } 240 }
217 } 241 }
218 242
219 // make G4LorentzVector for parent in its re 243 // make G4LorentzVector for parent in its rest frame
220 mass = theParentParticle->GetMass(); 244 mass = theParentParticle->GetMass();
221 G4LorentzVector parent4(0.0, 0.0, 0.0, mass) << 245 G4LorentzVector parent4( 0.0, 0.0, 0.0, mass);
222 246
223 // boost parent to new frame << 247 // boost parent to new frame
224 parent4.boost(newbetax, newbetay, newbetaz); 248 parent4.boost(newbetax, newbetay, newbetaz);
225 249
226 // change energy/momentum 250 // change energy/momentum
227 theParentParticle->Set4Momentum(parent4); 251 theParentParticle->Set4Momentum(parent4);
228 } 252 }
229 253
230 G4bool G4DecayProducts::IsChecked() const 254 G4bool G4DecayProducts::IsChecked() const
231 { 255 {
232 G4bool returnValue = true; 256 G4bool returnValue = true;
233 257
234 // check parent << 258 // check parent
235 // energy/momentum 259 // energy/momentum
236 G4double parent_energy = theParentParticle-> << 260 G4double parent_energy = theParentParticle->GetTotalEnergy();
237 G4ThreeVector direction = theParentParticle- 261 G4ThreeVector direction = theParentParticle->GetMomentumDirection();
238 G4ThreeVector parent_momentum = direction * << 262 G4ThreeVector parent_momentum
>> 263 = direction*(theParentParticle->GetTotalMomentum());
239 264
240 // check momentum direction is a unit vector 265 // check momentum direction is a unit vector
241 if ((parent_momentum.mag() > 0.0) && (std::f << 266 if ((parent_momentum.mag() >0.0) && (std::fabs(direction.mag()-1.0) >1.0e-6))
>> 267 {
242 #ifdef G4VERBOSE 268 #ifdef G4VERBOSE
243 G4cout << "G4DecayProducts::IsChecked():: 269 G4cout << "G4DecayProducts::IsChecked():: "
244 << " Momentum Direction Vector of P 270 << " Momentum Direction Vector of Parent is not normalized "
245 << " (=" << direction.mag() << ")" 271 << " (=" << direction.mag() << ")" << G4endl;
246 #endif 272 #endif
247 returnValue = false; 273 returnValue = false;
248 parent_momentum = parent_momentum * (1. / << 274 parent_momentum = parent_momentum * (1./direction.mag());
249 } 275 }
250 276
251 // daughters 277 // daughters
252 G4double mass, energy; 278 G4double mass, energy;
253 G4ThreeVector momentum; 279 G4ThreeVector momentum;
254 G4double total_energy = parent_energy; 280 G4double total_energy = parent_energy;
255 G4ThreeVector total_momentum = parent_moment 281 G4ThreeVector total_momentum = parent_momentum;
256 282
257 for (G4int index = 0; index < numberOfProduc << 283 for (G4int index=0; index < numberOfProducts; ++index)
>> 284 {
258 G4DynamicParticle* part = theProductVector 285 G4DynamicParticle* part = theProductVector->at(index);
259 mass = part->GetMass(); 286 mass = part->GetMass();
260 energy = part->GetTotalEnergy(); 287 energy = part->GetTotalEnergy();
261 direction = part->GetMomentumDirection(); 288 direction = part->GetMomentumDirection();
262 momentum = direction * (part->GetTotalMome << 289 momentum = direction*(part->GetTotalMomentum());
263 290
264 // check momentum direction is a unit vect 291 // check momentum direction is a unit vector
265 if ((momentum.mag() > 0.0) && (std::fabs(d << 292 if ( (momentum.mag()>0.0) && (std::fabs(direction.mag()-1.0) > 1.0e-6))
>> 293 {
266 #ifdef G4VERBOSE 294 #ifdef G4VERBOSE
267 G4cout << "G4DecayProducts::IsChecked(): << 295 G4cout << "G4DecayProducts::IsChecked():: "
268 << " Momentum Direction Vector of 296 << " Momentum Direction Vector of Daughter [" << index
269 << "] is not normalized (=" << d 297 << "] is not normalized (=" << direction.mag() << ")" << G4endl;
270 #endif 298 #endif
271 returnValue = false; 299 returnValue = false;
272 momentum = momentum * (1. / direction.ma << 300 momentum = momentum * (1./direction.mag());
273 } 301 }
274 // whether daughter stops or not 302 // whether daughter stops or not
275 if (energy - mass < DBL_MIN) { << 303 if (energy - mass < DBL_MIN )
>> 304 {
276 #ifdef G4VERBOSE 305 #ifdef G4VERBOSE
277 G4cout << "G4DecayProducts::IsChecked(): << 306 G4cout << "G4DecayProducts::IsChecked():: "
278 << " Daughter [" << index << "] << 307 << " Daughter [" << index << "] has no kinetic energy "<< G4endl;
279 #endif 308 #endif
280 returnValue = false; 309 returnValue = false;
281 } 310 }
282 total_energy -= energy; << 311 total_energy -= energy;
283 total_momentum -= momentum; 312 total_momentum -= momentum;
284 } 313 }
285 // check energy/momentum conservation 314 // check energy/momentum conservation
286 if ((std::fabs(total_energy) > 1.0e-9 * MeV) << 315 if ( (std::fabs(total_energy) >1.0e-9*MeV)
>> 316 || (total_momentum.mag() >1.0e-9*MeV ) )
>> 317 {
287 #ifdef G4VERBOSE 318 #ifdef G4VERBOSE
288 G4cout << "G4DecayProducts::IsChecked():: << 319 G4cout << "G4DecayProducts::IsChecked():: "
289 << " Energy/Momentum is not conserv << 320 << " Energy/Momentum is not conserved "<< G4endl;
290 G4cout << " difference between parent ener << 321 G4cout << " difference between parent energy & sum of daughters energy: "
291 << "[MeV] " << G4endl; << 322 << total_energy /MeV << "[MeV] " << G4endl;
292 G4cout << " difference between parent mome 323 G4cout << " difference between parent momentum & sum of daughters momentum: "
293 << " x:" << total_momentum.getX() / << 324 << " x:" << total_momentum.getX()/MeV
294 << " z:" << total_momentum.getZ() / << 325 << " y:" << total_momentum.getY()/MeV
>> 326 << " z:" << total_momentum.getZ()/MeV
>> 327 << G4endl;
295 #endif 328 #endif
296 returnValue = false; 329 returnValue = false;
297 } 330 }
298 return returnValue; 331 return returnValue;
299 } 332 }
300 333
301 void G4DecayProducts::DumpInfo() const 334 void G4DecayProducts::DumpInfo() const
302 { 335 {
303 G4cout << " ----- List of DecayProducts --- << 336 G4cout << " ----- List of DecayProducts -----" << G4endl;
304 G4cout << " ------ Parent Particle --------- << 337 G4cout << " ------ Parent Particle ----------" << G4endl;
305 if (theParentParticle != nullptr) theParentP << 338 if (theParentParticle != 0) theParentParticle->DumpInfo();
306 G4cout << " ------ Daughter Particles ----- << 339 G4cout << " ------ Daughter Particles ------" << G4endl;
307 for (G4int index = 0; index < numberOfProduc << 340 for (G4int index=0; index<numberOfProducts; ++index)
308 G4cout << " ----------" << index + 1 << " << 341 {
309 (theProductVector->at(index))->DumpInfo(); << 342 G4cout << " ----------" << index+1 << " -------------" << G4endl;
310 } << 343 (theProductVector->at(index))-> DumpInfo();
311 G4cout << " ----- End List of DecayProducts << 344 }
312 G4cout << G4endl; << 345 G4cout << " ----- End List of DecayProducts -----" << G4endl;
313 } << 346 G4cout << G4endl;
>> 347 }
314 348