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