| 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 // J.L. Chuma, TRIUMF, 31-Oct-1996 26 // J.L. Chuma, TRIUMF, 31-Oct-1996
27 // last modified: 19-Dec-1996 27 // last modified: 19-Dec-1996
28 // Modified by J.L.Chuma, 05-May-97 28 // Modified by J.L.Chuma, 05-May-97
29 // M. Kelsey 29-Aug-2011 -- Use G4Allocator fo 29 // M. Kelsey 29-Aug-2011 -- Use G4Allocator for better memory management
30 30
31 #include "G4ReactionProduct.hh" 31 #include "G4ReactionProduct.hh"
32 32
33 G4Allocator<G4ReactionProduct>*& aRPAllocator( 33 G4Allocator<G4ReactionProduct>*& aRPAllocator()
34 { 34 {
35 G4ThreadLocalStatic G4Allocator<G4Reaction 35 G4ThreadLocalStatic G4Allocator<G4ReactionProduct>* _instance = nullptr;
36 return _instance; 36 return _instance;
37 } 37 }
38 38
39 G4ReactionProduct::G4ReactionProduct() : 39 G4ReactionProduct::G4ReactionProduct() :
40 theParticleDefinition(NULL), 40 theParticleDefinition(NULL),
41 formationTime(0.0), 41 formationTime(0.0),
42 hasInitialStateParton(false), 42 hasInitialStateParton(false),
43 mass(0.0), 43 mass(0.0),
44 totalEnergy(0.0), 44 totalEnergy(0.0),
45 kineticEnergy(0.0), 45 kineticEnergy(0.0),
46 timeOfFlight(0.0), 46 timeOfFlight(0.0),
47 side(0), 47 side(0),
48 theCreatorModel(-1), 48 theCreatorModel(-1),
49 theParentResonanceDef(nullptr), 49 theParentResonanceDef(nullptr),
50 theParentResonanceID(0), 50 theParentResonanceID(0),
51 NewlyAdded(false), 51 NewlyAdded(false),
52 MayBeKilled(true) 52 MayBeKilled(true)
53 { 53 {
54 SetMomentum( 0.0, 0.0, 0.0 ); 54 SetMomentum( 0.0, 0.0, 0.0 );
55 SetPositionInNucleus( 0.0, 0.0, 0.0 ); 55 SetPositionInNucleus( 0.0, 0.0, 0.0 );
56 } 56 }
57 57
58 G4ReactionProduct::G4ReactionProduct( 58 G4ReactionProduct::G4ReactionProduct(
59 const G4ParticleDefinition *aParticleDefinit 59 const G4ParticleDefinition *aParticleDefinition )
60 { 60 {
61 SetMomentum( 0.0, 0.0, 0.0 ); 61 SetMomentum( 0.0, 0.0, 0.0 );
62 SetPositionInNucleus( 0.0, 0.0, 0.0 ); 62 SetPositionInNucleus( 0.0, 0.0, 0.0 );
63 formationTime = 0.0; 63 formationTime = 0.0;
64 hasInitialStateParton = false; 64 hasInitialStateParton = false;
65 theParticleDefinition = aParticleDefinitio 65 theParticleDefinition = aParticleDefinition;
66 mass = aParticleDefinition->GetPDGMass(); 66 mass = aParticleDefinition->GetPDGMass();
67 totalEnergy = mass; 67 totalEnergy = mass;
68 kineticEnergy = 0.0; 68 kineticEnergy = 0.0;
69 (aParticleDefinition->GetPDGEncoding()<0) 69 (aParticleDefinition->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
70 side = 0; 70 side = 0;
71 theCreatorModel = -1; 71 theCreatorModel = -1;
72 theParentResonanceDef = nullptr; 72 theParentResonanceDef = nullptr;
73 theParentResonanceID = 0; 73 theParentResonanceID = 0;
74 NewlyAdded = false; 74 NewlyAdded = false;
75 MayBeKilled = true; 75 MayBeKilled = true;
76 } 76 }
77 77
78 G4ReactionProduct::G4ReactionProduct( 78 G4ReactionProduct::G4ReactionProduct(
79 const G4ReactionProduct &right ) 79 const G4ReactionProduct &right )
80 { 80 {
81 theParticleDefinition = right.theParticleD 81 theParticleDefinition = right.theParticleDefinition;
82 positionInNucleus = right.positionInNucleu 82 positionInNucleus = right.positionInNucleus;
83 formationTime = right.formationTime; 83 formationTime = right.formationTime;
84 hasInitialStateParton = right.hasInitialSt 84 hasInitialStateParton = right.hasInitialStateParton;
85 momentum = right.momentum; 85 momentum = right.momentum;
86 mass = right.mass; 86 mass = right.mass;
87 totalEnergy = right.totalEnergy; 87 totalEnergy = right.totalEnergy;
88 kineticEnergy = right.kineticEnergy; 88 kineticEnergy = right.kineticEnergy;
89 timeOfFlight = right.timeOfFlight; 89 timeOfFlight = right.timeOfFlight;
90 side = right.side; 90 side = right.side;
91 theCreatorModel = right.theCreatorModel; 91 theCreatorModel = right.theCreatorModel;
92 theParentResonanceDef = right.theParentRes 92 theParentResonanceDef = right.theParentResonanceDef;
93 theParentResonanceID = right.theParentReso 93 theParentResonanceID = right.theParentResonanceID;
94 NewlyAdded = right.NewlyAdded; 94 NewlyAdded = right.NewlyAdded;
95 MayBeKilled = right.MayBeKilled; 95 MayBeKilled = right.MayBeKilled;
96 } 96 }
97 97
98 G4ReactionProduct &G4ReactionProduct::operato 98 G4ReactionProduct &G4ReactionProduct::operator=(
99 const G4ReactionProduct &right ) 99 const G4ReactionProduct &right )
100 { 100 {
101 if( this != &right ) { 101 if( this != &right ) {
102 theParticleDefinition = right.theParticl 102 theParticleDefinition = right.theParticleDefinition;
103 positionInNucleus = right.positionInNucl 103 positionInNucleus = right.positionInNucleus;
104 formationTime = right.formationTime; 104 formationTime = right.formationTime;
105 hasInitialStateParton = right.hasInitial 105 hasInitialStateParton = right.hasInitialStateParton;
106 momentum = right.momentum; 106 momentum = right.momentum;
107 mass = right.mass; 107 mass = right.mass;
108 totalEnergy = right.totalEnergy; 108 totalEnergy = right.totalEnergy;
109 kineticEnergy = right.kineticEnergy; 109 kineticEnergy = right.kineticEnergy;
110 timeOfFlight = right.timeOfFlight; 110 timeOfFlight = right.timeOfFlight;
111 side = right.side; 111 side = right.side;
112 theCreatorModel = right.theCreatorModel; 112 theCreatorModel = right.theCreatorModel;
113 theParentResonanceDef = right.theParentR 113 theParentResonanceDef = right.theParentResonanceDef;
114 theParentResonanceID = right.theParentRe 114 theParentResonanceID = right.theParentResonanceID;
115 NewlyAdded = right.NewlyAdded; 115 NewlyAdded = right.NewlyAdded;
116 MayBeKilled = right.MayBeKilled; 116 MayBeKilled = right.MayBeKilled;
117 } 117 }
118 return *this; 118 return *this;
119 } 119 }
120 120
121 G4ReactionProduct &G4ReactionProduct::operato 121 G4ReactionProduct &G4ReactionProduct::operator=(
122 const G4DynamicParticle &right ) 122 const G4DynamicParticle &right )
123 { 123 {
124 theParticleDefinition = right.GetDefinitio 124 theParticleDefinition = right.GetDefinition();
125 SetPositionInNucleus( 0.0, 0.0, 0.0 ); 125 SetPositionInNucleus( 0.0, 0.0, 0.0 );
126 formationTime = 0.0; 126 formationTime = 0.0;
127 hasInitialStateParton = false; 127 hasInitialStateParton = false;
128 momentum = right.GetMomentum(); 128 momentum = right.GetMomentum();
129 mass = right.GetDefinition()->GetPDGMass() 129 mass = right.GetDefinition()->GetPDGMass();
130 totalEnergy = right.GetTotalEnergy(); 130 totalEnergy = right.GetTotalEnergy();
131 kineticEnergy = right.GetKineticEnergy(); 131 kineticEnergy = right.GetKineticEnergy();
132 (right.GetDefinition()->GetPDGEncoding()<0 132 (right.GetDefinition()->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
133 side = 0; 133 side = 0;
134 theCreatorModel = -1; 134 theCreatorModel = -1;
135 theParentResonanceDef = nullptr; 135 theParentResonanceDef = nullptr;
136 theParentResonanceID = 0; 136 theParentResonanceID = 0;
137 NewlyAdded = false; 137 NewlyAdded = false;
138 MayBeKilled = true; 138 MayBeKilled = true;
139 return *this; 139 return *this;
140 } 140 }
141 141
142 G4ReactionProduct &G4ReactionProduct::operato 142 G4ReactionProduct &G4ReactionProduct::operator=(
143 const G4HadProjectile &right ) 143 const G4HadProjectile &right )
144 { 144 {
145 theParticleDefinition = right.GetDefinitio 145 theParticleDefinition = right.GetDefinition();
146 SetPositionInNucleus( 0.0, 0.0, 0.0 ); 146 SetPositionInNucleus( 0.0, 0.0, 0.0 );
147 formationTime = 0.0; 147 formationTime = 0.0;
148 hasInitialStateParton = false; 148 hasInitialStateParton = false;
149 momentum = right.Get4Momentum().vect(); 149 momentum = right.Get4Momentum().vect();
150 mass = right.GetDefinition()->GetPDGMass() 150 mass = right.GetDefinition()->GetPDGMass();
151 totalEnergy = right.Get4Momentum().e(); 151 totalEnergy = right.Get4Momentum().e();
152 kineticEnergy = right.GetKineticEnergy(); 152 kineticEnergy = right.GetKineticEnergy();
153 (right.GetDefinition()->GetPDGEncoding()<0 153 (right.GetDefinition()->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
154 side = 0; 154 side = 0;
155 theCreatorModel = -1; 155 theCreatorModel = -1;
156 theParentResonanceDef = nullptr; 156 theParentResonanceDef = nullptr;
157 theParentResonanceID = 0; 157 theParentResonanceID = 0;
158 NewlyAdded = false; 158 NewlyAdded = false;
159 MayBeKilled = true; 159 MayBeKilled = true;
160 return *this; 160 return *this;
161 } 161 }
162 162
163 void G4ReactionProduct::SetDefinitionAndUpdat 163 void G4ReactionProduct::SetDefinitionAndUpdateE(
164 const G4ParticleDefinition *aParticleDefinit 164 const G4ParticleDefinition *aParticleDefinition )
165 { G4double aKineticEnergy = GetKineticEne 165 { G4double aKineticEnergy = GetKineticEnergy();
166 G4double pp = GetMomentum().mag(); 166 G4double pp = GetMomentum().mag();
167 G4ThreeVector aMomentum = GetMomentum(); 167 G4ThreeVector aMomentum = GetMomentum();
168 SetDefinition( aParticleDefinition ); 168 SetDefinition( aParticleDefinition );
169 SetKineticEnergy( aKineticEnergy ); 169 SetKineticEnergy( aKineticEnergy );
170 if( pp > DBL_MIN ) 170 if( pp > DBL_MIN )
171 SetMomentum( aMomentum * (std::sqrt(aKin 171 SetMomentum( aMomentum * (std::sqrt(aKineticEnergy*aKineticEnergy +
172 2*aKinetic 172 2*aKineticEnergy*GetMass())/pp) );
173 } 173 }
174 174
175 void G4ReactionProduct::SetDefinition( 175 void G4ReactionProduct::SetDefinition(
176 const G4ParticleDefinition *aParticleDefinit 176 const G4ParticleDefinition *aParticleDefinition )
177 { 177 {
178 theParticleDefinition = aParticleDefinitio 178 theParticleDefinition = aParticleDefinition;
179 mass = aParticleDefinition->GetPDGMass(); 179 mass = aParticleDefinition->GetPDGMass();
180 totalEnergy = mass; 180 totalEnergy = mass;
181 kineticEnergy = 0.0; 181 kineticEnergy = 0.0;
182 (aParticleDefinition->GetPDGEncoding()<0) 182 (aParticleDefinition->GetPDGEncoding()<0) ?
183 timeOfFlight=-1.0 : timeOfFlight=1.0; 183 timeOfFlight=-1.0 : timeOfFlight=1.0;
184 } 184 }
185 185
186 void G4ReactionProduct::SetMomentum( 186 void G4ReactionProduct::SetMomentum(
187 const G4double x, const G4double y, const G4 187 const G4double x, const G4double y, const G4double z )
188 { 188 {
189 momentum.setX( x ); 189 momentum.setX( x );
190 momentum.setY( y ); 190 momentum.setY( y );
191 momentum.setZ( z ); 191 momentum.setZ( z );
192 } 192 }
193 193
194 void G4ReactionProduct::SetMomentum( 194 void G4ReactionProduct::SetMomentum(
195 const G4double x, const G4double y ) 195 const G4double x, const G4double y )
196 { 196 {
197 momentum.setX( x ); 197 momentum.setX( x );
198 momentum.setY( y ); 198 momentum.setY( y );
199 } 199 }
200 200
201 void G4ReactionProduct::SetMomentum( const G4 201 void G4ReactionProduct::SetMomentum( const G4double z )
202 { 202 {
203 momentum.setZ( z ); 203 momentum.setZ( z );
204 } 204 }
205 205
206 void G4ReactionProduct::SetZero() 206 void G4ReactionProduct::SetZero()
207 { 207 {
208 SetMomentum( 0.0, 0.0, 0.0 ); 208 SetMomentum( 0.0, 0.0, 0.0 );
209 totalEnergy = 0.0; 209 totalEnergy = 0.0;
210 kineticEnergy = 0.0; 210 kineticEnergy = 0.0;
211 mass = 0.0; 211 mass = 0.0;
212 timeOfFlight = 0.0; 212 timeOfFlight = 0.0;
213 side = 0; 213 side = 0;
214 theCreatorModel = -1; 214 theCreatorModel = -1;
215 theParentResonanceDef = nullptr; 215 theParentResonanceDef = nullptr;
216 theParentResonanceID = 0; 216 theParentResonanceID = 0;
217 NewlyAdded = false; 217 NewlyAdded = false;
218 SetPositionInNucleus( 0.0, 0.0, 0.0 ); 218 SetPositionInNucleus( 0.0, 0.0, 0.0 );
219 formationTime = 0.0; 219 formationTime = 0.0;
220 hasInitialStateParton = false; 220 hasInitialStateParton = false;
221 } 221 }
222 222
223 void G4ReactionProduct::Lorentz( 223 void G4ReactionProduct::Lorentz(
224 const G4ReactionProduct &p1, const G4Reacti 224 const G4ReactionProduct &p1, const G4ReactionProduct &p2 )
225 { 225 {
226 G4ThreeVector p1M = p1.momentum; 226 G4ThreeVector p1M = p1.momentum;
227 G4ThreeVector p2M = p2.momentum; 227 G4ThreeVector p2M = p2.momentum;
228 G4double p1x = p1M.x(); G4double p1y = p1M 228 G4double p1x = p1M.x(); G4double p1y = p1M.y(); G4double p1z = p1M.z();
229 G4double p2x = p2M.x(); G4double p2y = p2M 229 G4double p2x = p2M.x(); G4double p2y = p2M.y(); G4double p2z = p2M.z();
230 G4double a = ( (p1x*p2x+p1y*p2y+p1z*p2z)/( 230 G4double a = ( (p1x*p2x+p1y*p2y+p1z*p2z)/(p2.totalEnergy+p2.mass) -
231 p1.totalEnergy ) / p2.mass; 231 p1.totalEnergy ) / p2.mass;
232 G4double x = p1x+a*p2x; 232 G4double x = p1x+a*p2x;
233 G4double y = p1y+a*p2y; 233 G4double y = p1y+a*p2y;
234 G4double z = p1z+a*p2z; 234 G4double z = p1z+a*p2z;
235 G4double p = std::sqrt(x*x+y*y+z*z); 235 G4double p = std::sqrt(x*x+y*y+z*z);
236 SetMass( p1.mass ); 236 SetMass( p1.mass );
237 SetTotalEnergy( std::sqrt( (p1.mass+p)*(p1 237 SetTotalEnergy( std::sqrt( (p1.mass+p)*(p1.mass+p) - 2.*p1.mass*p ) );
238 //SetTotalEnergy( std::sqrt( p1.mass*p1.ma 238 //SetTotalEnergy( std::sqrt( p1.mass*p1.mass + x*x + y*y + z*z ) );
239 SetMomentum( x, y, z ); 239 SetMomentum( x, y, z );
240 } 240 }
241 241
242 G4double G4ReactionProduct::Angle( 242 G4double G4ReactionProduct::Angle(
243 const G4ReactionProduct& p ) const 243 const G4ReactionProduct& p ) const
244 { 244 {
245 G4ThreeVector tM = momentum; 245 G4ThreeVector tM = momentum;
246 G4ThreeVector pM = p.momentum; 246 G4ThreeVector pM = p.momentum;
247 G4double tx = tM.x(); G4double ty = tM.y() 247 G4double tx = tM.x(); G4double ty = tM.y(); G4double tz = tM.z();
248 G4double px = pM.x(); G4double py = pM.y() 248 G4double px = pM.x(); G4double py = pM.y(); G4double pz = pM.z();
249 G4double a = std::sqrt( ( px*px + py*py + 249 G4double a = std::sqrt( ( px*px + py*py + pz*pz ) * ( tx*tx + ty*ty + tz*tz ) );
250 if( a == 0.0 ) { 250 if( a == 0.0 ) {
251 return 0.0; 251 return 0.0;
252 } else { 252 } else {
253 a = ( tx*px + ty*py + tz*pz ) / a; 253 a = ( tx*px + ty*py + tz*pz ) / a;
254 if( std::abs(a) > 1.0 ) { a<0.0 ? a=-1.0 254 if( std::abs(a) > 1.0 ) { a<0.0 ? a=-1.0 : a=1.0; }
255 return std::acos( a ); 255 return std::acos( a );
256 } 256 }
257 } 257 }
258 258
259 G4ReactionProduct operator+( 259 G4ReactionProduct operator+(
260 const G4ReactionProduct& p1, const G4Reactio 260 const G4ReactionProduct& p1, const G4ReactionProduct& p2 )
261 { 261 {
262 G4double totEnergy = p1.totalEnergy + p2.t 262 G4double totEnergy = p1.totalEnergy + p2.totalEnergy;
263 G4double x = p1.momentum.x() + p2.momentum 263 G4double x = p1.momentum.x() + p2.momentum.x();
264 G4double y = p1.momentum.y() + p2.momentum 264 G4double y = p1.momentum.y() + p2.momentum.y();
265 G4double z = p1.momentum.z() + p2.momentum 265 G4double z = p1.momentum.z() + p2.momentum.z();
266 G4double newMass = totEnergy*totEnergy - ( 266 G4double newMass = totEnergy*totEnergy - ( x*x + y*y + z*z );
267 if( newMass < 0.0 ) 267 if( newMass < 0.0 )
268 newMass = -1. * std::sqrt( -newMass ); 268 newMass = -1. * std::sqrt( -newMass );
269 else 269 else
270 newMass = std::sqrt( newMass ); 270 newMass = std::sqrt( newMass );
271 G4ReactionProduct result; 271 G4ReactionProduct result;
272 result.SetMass( newMass ); 272 result.SetMass( newMass );
273 result.SetMomentum( x, y, z ); 273 result.SetMomentum( x, y, z );
274 result.SetTotalEnergy( totEnergy ); 274 result.SetTotalEnergy( totEnergy );
275 result.SetPositionInNucleus( 0.0, 0.0, 0.0 275 result.SetPositionInNucleus( 0.0, 0.0, 0.0 );
276 result.SetFormationTime(0.0); 276 result.SetFormationTime(0.0);
277 result.HasInitialStateParton(false); 277 result.HasInitialStateParton(false);
278 return result; 278 return result;
279 } 279 }
280 280
281 G4ReactionProduct operator-( 281 G4ReactionProduct operator-(
282 const G4ReactionProduct& p1, const G4Reactio 282 const G4ReactionProduct& p1, const G4ReactionProduct& p2 )
283 { 283 {
284 G4double totEnergy = p1.totalEnergy - p2.t 284 G4double totEnergy = p1.totalEnergy - p2.totalEnergy;
285 G4double x = p1.momentum.x() - p2.momentum 285 G4double x = p1.momentum.x() - p2.momentum.x();
286 G4double y = p1.momentum.y() - p2.momentum 286 G4double y = p1.momentum.y() - p2.momentum.y();
287 G4double z = p1.momentum.z() - p2.momentum 287 G4double z = p1.momentum.z() - p2.momentum.z();
288 G4double newMass = totEnergy*totEnergy - ( 288 G4double newMass = totEnergy*totEnergy - ( x*x + y*y + z*z );
289 if( newMass < 0.0 ) 289 if( newMass < 0.0 )
290 newMass = -1. * std::sqrt( -newMass ); 290 newMass = -1. * std::sqrt( -newMass );
291 else 291 else
292 newMass = std::sqrt( newMass ); 292 newMass = std::sqrt( newMass );
293 G4ReactionProduct result; 293 G4ReactionProduct result;
294 result.SetMass( newMass ); 294 result.SetMass( newMass );
295 result.SetMomentum( x, y, z ); 295 result.SetMomentum( x, y, z );
296 result.SetTotalEnergy( totEnergy ); 296 result.SetTotalEnergy( totEnergy );
297 result.SetPositionInNucleus( 0.0, 0.0, 0.0 297 result.SetPositionInNucleus( 0.0, 0.0, 0.0 );
298 result.SetFormationTime(0.0); 298 result.SetFormationTime(0.0);
299 result.HasInitialStateParton(false); 299 result.HasInitialStateParton(false);
300 return result; 300 return result;
301 } 301 }
302 /* end of code */ 302 /* end of code */
303 303
304 304