00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "G4ReactionProduct.hh"
00032
00033 G4Allocator<G4ReactionProduct> aRPAllocator;
00034
00035
00036 G4ReactionProduct::G4ReactionProduct() :
00037 theParticleDefinition(NULL),
00038 formationTime(0.0),
00039 hasInitialStateParton(false),
00040 mass(0.0),
00041 totalEnergy(0.0),
00042 kineticEnergy(0.0),
00043 timeOfFlight(0.0),
00044 side(0),
00045 NewlyAdded(false),
00046 MayBeKilled(true)
00047 {
00048 SetMomentum( 0.0, 0.0, 0.0 );
00049 SetPositionInNucleus( 0.0, 0.0, 0.0 );
00050 }
00051
00052 G4ReactionProduct::G4ReactionProduct(
00053 G4ParticleDefinition *aParticleDefinition )
00054 {
00055 SetMomentum( 0.0, 0.0, 0.0 );
00056 SetPositionInNucleus( 0.0, 0.0, 0.0 );
00057 formationTime = 0.0;
00058 hasInitialStateParton = false;
00059 theParticleDefinition = aParticleDefinition;
00060 mass = aParticleDefinition->GetPDGMass();
00061 totalEnergy = mass;
00062 kineticEnergy = 0.0;
00063 (aParticleDefinition->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
00064 side = 0;
00065 NewlyAdded = false;
00066 MayBeKilled = true;
00067 }
00068
00069 G4ReactionProduct::G4ReactionProduct(
00070 const G4ReactionProduct &right )
00071 {
00072 theParticleDefinition = right.theParticleDefinition;
00073 positionInNucleus = right.positionInNucleus;
00074 formationTime = right.formationTime;
00075 hasInitialStateParton = right.hasInitialStateParton;
00076 momentum = right.momentum;
00077 mass = right.mass;
00078 totalEnergy = right.totalEnergy;
00079 kineticEnergy = right.kineticEnergy;
00080 timeOfFlight = right.timeOfFlight;
00081 side = right.side;
00082 NewlyAdded = right.NewlyAdded;
00083 MayBeKilled = right.MayBeKilled;
00084 }
00085
00086 G4ReactionProduct &G4ReactionProduct::operator=(
00087 const G4ReactionProduct &right )
00088 {
00089 if( this != &right ) {
00090 theParticleDefinition = right.theParticleDefinition;
00091 positionInNucleus = right.positionInNucleus;
00092 formationTime = right.formationTime;
00093 hasInitialStateParton = right.hasInitialStateParton;
00094 momentum = right.momentum;
00095 mass = right.mass;
00096 totalEnergy = right.totalEnergy;
00097 kineticEnergy = right.kineticEnergy;
00098 timeOfFlight = right.timeOfFlight;
00099 side = right.side;
00100 NewlyAdded = right.NewlyAdded;
00101 MayBeKilled = right.MayBeKilled;
00102 }
00103 return *this;
00104 }
00105
00106 G4ReactionProduct &G4ReactionProduct::operator=(
00107 const G4DynamicParticle &right )
00108 {
00109 theParticleDefinition = right.GetDefinition();
00110 SetPositionInNucleus( 0.0, 0.0, 0.0 );
00111 formationTime = 0.0;
00112 hasInitialStateParton = false;
00113 momentum = right.GetMomentum();
00114 mass = right.GetDefinition()->GetPDGMass();
00115 totalEnergy = right.GetTotalEnergy();
00116 kineticEnergy = right.GetKineticEnergy();
00117 (right.GetDefinition()->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
00118 side = 0;
00119 NewlyAdded = false;
00120 MayBeKilled = true;
00121 return *this;
00122 }
00123
00124 G4ReactionProduct &G4ReactionProduct::operator=(
00125 const G4HadProjectile &right )
00126 {
00127 theParticleDefinition = const_cast<G4ParticleDefinition *>(right.GetDefinition());
00128 SetPositionInNucleus( 0.0, 0.0, 0.0 );
00129 formationTime = 0.0;
00130 hasInitialStateParton = false;
00131 momentum = right.Get4Momentum().vect();
00132 mass = right.GetDefinition()->GetPDGMass();
00133 totalEnergy = right.Get4Momentum().e();
00134 kineticEnergy = right.GetKineticEnergy();
00135 (right.GetDefinition()->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
00136 side = 0;
00137 NewlyAdded = false;
00138 MayBeKilled = true;
00139 return *this;
00140 }
00141
00142 void G4ReactionProduct::SetDefinitionAndUpdateE(
00143 G4ParticleDefinition *aParticleDefinition )
00144 {
00145 G4double aKineticEnergy = GetKineticEnergy();
00146 G4double pp = GetMomentum().mag();
00147 G4ThreeVector aMomentum = GetMomentum();
00148 SetDefinition( aParticleDefinition );
00149 SetKineticEnergy( aKineticEnergy );
00150 if( pp > DBL_MIN )
00151 SetMomentum( aMomentum * (std::sqrt(aKineticEnergy*aKineticEnergy +
00152 2*aKineticEnergy*GetMass())/pp) );
00153 }
00154
00155 void G4ReactionProduct::SetDefinition(
00156 G4ParticleDefinition *aParticleDefinition )
00157 {
00158 theParticleDefinition = aParticleDefinition;
00159 mass = aParticleDefinition->GetPDGMass();
00160 totalEnergy = mass;
00161 kineticEnergy = 0.0;
00162 (aParticleDefinition->GetPDGEncoding()<0) ?
00163 timeOfFlight=-1.0 : timeOfFlight=1.0;
00164 }
00165
00166 void G4ReactionProduct::SetMomentum(
00167 const G4double x, const G4double y, const G4double z )
00168 {
00169 momentum.setX( x );
00170 momentum.setY( y );
00171 momentum.setZ( z );
00172 }
00173
00174 void G4ReactionProduct::SetMomentum(
00175 const G4double x, const G4double y )
00176 {
00177 momentum.setX( x );
00178 momentum.setY( y );
00179 }
00180
00181 void G4ReactionProduct::SetMomentum( const G4double z )
00182 {
00183 momentum.setZ( z );
00184 }
00185
00186 void G4ReactionProduct::SetZero()
00187 {
00188 SetMomentum( 0.0, 0.0, 0.0 );
00189 totalEnergy = 0.0;
00190 kineticEnergy = 0.0;
00191 mass = 0.0;
00192 timeOfFlight = 0.0;
00193 side = 0;
00194 NewlyAdded = false;
00195 SetPositionInNucleus( 0.0, 0.0, 0.0 );
00196 formationTime = 0.0;
00197 hasInitialStateParton = false;
00198 }
00199
00200 void G4ReactionProduct::Lorentz(
00201 const G4ReactionProduct &p1, const G4ReactionProduct &p2 )
00202 {
00203 G4ThreeVector p1M = p1.momentum;
00204 G4ThreeVector p2M = p2.momentum;
00205 G4double p1x = p1M.x(); G4double p1y = p1M.y(); G4double p1z = p1M.z();
00206 G4double p2x = p2M.x(); G4double p2y = p2M.y(); G4double p2z = p2M.z();
00207 G4double a = ( (p1x*p2x+p1y*p2y+p1z*p2z)/(p2.totalEnergy+p2.mass) -
00208 p1.totalEnergy ) / p2.mass;
00209 G4double x = p1x+a*p2x;
00210 G4double y = p1y+a*p2y;
00211 G4double z = p1z+a*p2z;
00212 G4double p = std::sqrt(x*x+y*y+z*z);
00213 SetMass( p1.mass );
00214 SetTotalEnergy( std::sqrt( (p1.mass+p)*(p1.mass+p) - 2.*p1.mass*p ) );
00215
00216 SetMomentum( x, y, z );
00217 }
00218
00219 G4double G4ReactionProduct::Angle(
00220 const G4ReactionProduct& p ) const
00221 {
00222 G4ThreeVector tM = momentum;
00223 G4ThreeVector pM = p.momentum;
00224 G4double tx = tM.x(); G4double ty = tM.y(); G4double tz = tM.z();
00225 G4double px = pM.x(); G4double py = pM.y(); G4double pz = pM.z();
00226 G4double a = std::sqrt( ( px*px + py*py + pz*pz ) * ( tx*tx + ty*ty + tz*tz ) );
00227 if( a == 0.0 ) {
00228 return 0.0;
00229 } else {
00230 a = ( tx*px + ty*py + tz*pz ) / a;
00231 if( std::fabs(a) > 1.0 ) { a<0.0 ? a=-1.0 : a=1.0; }
00232 return std::acos( a );
00233 }
00234 }
00235
00236 G4ReactionProduct operator+(
00237 const G4ReactionProduct& p1, const G4ReactionProduct& p2 )
00238 {
00239 G4double totEnergy = p1.totalEnergy + p2.totalEnergy;
00240 G4double x = p1.momentum.x() + p2.momentum.x();
00241 G4double y = p1.momentum.y() + p2.momentum.y();
00242 G4double z = p1.momentum.z() + p2.momentum.z();
00243 G4double newMass = totEnergy*totEnergy - ( x*x + y*y + z*z );
00244 if( newMass < 0.0 )
00245 newMass = -1. * std::sqrt( -newMass );
00246 else
00247 newMass = std::sqrt( newMass );
00248 G4ReactionProduct result;
00249 result.SetMass( newMass );
00250 result.SetMomentum( x, y, z );
00251 result.SetTotalEnergy( totEnergy );
00252 result.SetPositionInNucleus( 0.0, 0.0, 0.0 );
00253 result.SetFormationTime(0.0);
00254 result.HasInitialStateParton(false);
00255 return result;
00256 }
00257
00258 G4ReactionProduct operator-(
00259 const G4ReactionProduct& p1, const G4ReactionProduct& p2 )
00260 {
00261 G4double totEnergy = p1.totalEnergy - p2.totalEnergy;
00262 G4double x = p1.momentum.x() - p2.momentum.x();
00263 G4double y = p1.momentum.y() - p2.momentum.y();
00264 G4double z = p1.momentum.z() - p2.momentum.z();
00265 G4double newMass = totEnergy*totEnergy - ( x*x + y*y + z*z );
00266 if( newMass < 0.0 )
00267 newMass = -1. * std::sqrt( -newMass );
00268 else
00269 newMass = std::sqrt( newMass );
00270 G4ReactionProduct result;
00271 result.SetMass( newMass );
00272 result.SetMomentum( x, y, z );
00273 result.SetTotalEnergy( totEnergy );
00274 result.SetPositionInNucleus( 0.0, 0.0, 0.0 );
00275 result.SetFormationTime(0.0);
00276 result.HasInitialStateParton(false);
00277 return result;
00278 }
00279
00280