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