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