source: trunk/source/processes/hadronic/util/src/G4ReactionProduct.cc@ 1199

Last change on this file since 1199 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 9.3 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27//
28 // J.L. Chuma, TRIUMF, 31-Oct-1996
29 // last modified: 19-Dec-1996
30 // Modified by J.L.Chuma, 05-May-97
31
32#include "G4ReactionProduct.hh"
33
34 G4ReactionProduct::G4ReactionProduct() :
35 theParticleDefinition(NULL),
36 formationTime(0.0),
37 hasInitialStateParton(false),
38 mass(0.0),
39 totalEnergy(0.0),
40 kineticEnergy(0.0),
41 timeOfFlight(0.0),
42 side(0),
43 NewlyAdded(false),
44 MayBeKilled(true)
45 {
46 SetMomentum( 0.0, 0.0, 0.0 );
47 SetPositionInNucleus( 0.0, 0.0, 0.0 );
48 }
49
50 G4ReactionProduct::G4ReactionProduct(
51 G4ParticleDefinition *aParticleDefinition )
52 {
53 SetMomentum( 0.0, 0.0, 0.0 );
54 SetPositionInNucleus( 0.0, 0.0, 0.0 );
55 formationTime = 0.0;
56 hasInitialStateParton = false;
57 theParticleDefinition = aParticleDefinition;
58 mass = aParticleDefinition->GetPDGMass();
59 totalEnergy = mass;
60 kineticEnergy = 0.0;
61 (aParticleDefinition->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
62 side = 0;
63 NewlyAdded = false;
64 MayBeKilled = true;
65 }
66
67 G4ReactionProduct::G4ReactionProduct(
68 const G4ReactionProduct &right )
69 {
70 theParticleDefinition = right.theParticleDefinition;
71 positionInNucleus = right.positionInNucleus;
72 formationTime = right.formationTime;
73 hasInitialStateParton = right.hasInitialStateParton;
74 momentum = right.momentum;
75 mass = right.mass;
76 totalEnergy = right.totalEnergy;
77 kineticEnergy = right.kineticEnergy;
78 timeOfFlight = right.timeOfFlight;
79 side = right.side;
80 NewlyAdded = right.NewlyAdded;
81 MayBeKilled = right.MayBeKilled;
82 }
83
84 G4ReactionProduct &G4ReactionProduct::operator=(
85 const G4ReactionProduct &right )
86 {
87 if( this != &right ) {
88 theParticleDefinition = right.theParticleDefinition;
89 positionInNucleus = right.positionInNucleus;
90 formationTime = right.formationTime;
91 hasInitialStateParton = right.hasInitialStateParton;
92 momentum = right.momentum;
93 mass = right.mass;
94 totalEnergy = right.totalEnergy;
95 kineticEnergy = right.kineticEnergy;
96 timeOfFlight = right.timeOfFlight;
97 side = right.side;
98 NewlyAdded = right.NewlyAdded;
99 MayBeKilled = right.MayBeKilled;
100 }
101 return *this;
102 }
103
104 G4ReactionProduct &G4ReactionProduct::operator=(
105 const G4DynamicParticle &right )
106 {
107 theParticleDefinition = right.GetDefinition();
108 SetPositionInNucleus( 0.0, 0.0, 0.0 );
109 formationTime = 0.0;
110 hasInitialStateParton = false;
111 momentum = right.GetMomentum();
112 mass = right.GetDefinition()->GetPDGMass();
113 totalEnergy = right.GetTotalEnergy();
114 kineticEnergy = right.GetKineticEnergy();
115 (right.GetDefinition()->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
116 side = 0;
117 NewlyAdded = false;
118 MayBeKilled = true;
119 return *this;
120 }
121
122 G4ReactionProduct &G4ReactionProduct::operator=(
123 const G4HadProjectile &right )
124 {
125 theParticleDefinition = const_cast<G4ParticleDefinition *>(right.GetDefinition());
126 SetPositionInNucleus( 0.0, 0.0, 0.0 );
127 formationTime = 0.0;
128 hasInitialStateParton = false;
129 momentum = right.Get4Momentum().vect();
130 mass = right.GetDefinition()->GetPDGMass();
131 totalEnergy = right.Get4Momentum().e();
132 kineticEnergy = right.GetKineticEnergy();
133 (right.GetDefinition()->GetPDGEncoding()<0) ? timeOfFlight=-1.0 : timeOfFlight=1.0;
134 side = 0;
135 NewlyAdded = false;
136 MayBeKilled = true;
137 return *this;
138 }
139
140 void G4ReactionProduct::SetDefinitionAndUpdateE(
141 G4ParticleDefinition *aParticleDefinition )
142 {
143 G4double aKineticEnergy = GetKineticEnergy();
144 G4double pp = GetMomentum().mag();
145 G4ThreeVector aMomentum = GetMomentum();
146 SetDefinition( aParticleDefinition );
147 SetKineticEnergy( aKineticEnergy );
148 if( pp > DBL_MIN )
149 SetMomentum( aMomentum * (std::sqrt(aKineticEnergy*aKineticEnergy +
150 2*aKineticEnergy*GetMass())/pp) );
151 }
152
153 void G4ReactionProduct::SetDefinition(
154 G4ParticleDefinition *aParticleDefinition )
155 {
156 theParticleDefinition = aParticleDefinition;
157 mass = aParticleDefinition->GetPDGMass();
158 totalEnergy = mass;
159 kineticEnergy = 0.0;
160 (aParticleDefinition->GetPDGEncoding()<0) ?
161 timeOfFlight=-1.0 : timeOfFlight=1.0;
162 }
163
164 void G4ReactionProduct::SetMomentum(
165 const G4double x, const G4double y, const G4double z )
166 {
167 momentum.setX( x );
168 momentum.setY( y );
169 momentum.setZ( z );
170 }
171
172 void G4ReactionProduct::SetMomentum(
173 const G4double x, const G4double y )
174 {
175 momentum.setX( x );
176 momentum.setY( y );
177 }
178
179 void G4ReactionProduct::SetMomentum( const G4double z )
180 {
181 momentum.setZ( z );
182 }
183
184 void G4ReactionProduct::SetZero()
185 {
186 SetMomentum( 0.0, 0.0, 0.0 );
187 totalEnergy = 0.0;
188 kineticEnergy = 0.0;
189 mass = 0.0;
190 timeOfFlight = 0.0;
191 side = 0;
192 NewlyAdded = false;
193 SetPositionInNucleus( 0.0, 0.0, 0.0 );
194 formationTime = 0.0;
195 hasInitialStateParton = false;
196 }
197
198 void G4ReactionProduct::Lorentz(
199 const G4ReactionProduct &p1, const G4ReactionProduct &p2 )
200 {
201 G4ThreeVector p1M = p1.momentum;
202 G4ThreeVector p2M = p2.momentum;
203 G4double p1x = p1M.x(); G4double p1y = p1M.y(); G4double p1z = p1M.z();
204 G4double p2x = p2M.x(); G4double p2y = p2M.y(); G4double p2z = p2M.z();
205 G4double a = ( (p1x*p2x+p1y*p2y+p1z*p2z)/(p2.totalEnergy+p2.mass) -
206 p1.totalEnergy ) / p2.mass;
207 G4double x = p1x+a*p2x;
208 G4double y = p1y+a*p2y;
209 G4double z = p1z+a*p2z;
210 G4double p = std::sqrt(x*x+y*y+z*z);
211 SetMass( p1.mass );
212 SetTotalEnergy( std::sqrt( (p1.mass+p)*(p1.mass+p) - 2.*p1.mass*p ) );
213 //SetTotalEnergy( std::sqrt( p1.mass*p1.mass + x*x + y*y + z*z ) );
214 SetMomentum( x, y, z );
215 }
216
217 G4double G4ReactionProduct::Angle(
218 const G4ReactionProduct& p ) const
219 {
220 G4ThreeVector tM = momentum;
221 G4ThreeVector pM = p.momentum;
222 G4double tx = tM.x(); G4double ty = tM.y(); G4double tz = tM.z();
223 G4double px = pM.x(); G4double py = pM.y(); G4double pz = pM.z();
224 G4double a = std::sqrt( ( px*px + py*py + pz*pz ) * ( tx*tx + ty*ty + tz*tz ) );
225 if( a == 0.0 ) {
226 return 0.0;
227 } else {
228 a = ( tx*px + ty*py + tz*pz ) / a;
229 if( std::fabs(a) > 1.0 ) { a<0.0 ? a=-1.0 : a=1.0; }
230 return std::acos( a );
231 }
232 }
233
234 G4ReactionProduct operator+(
235 const G4ReactionProduct& p1, const G4ReactionProduct& p2 )
236 {
237 G4double totEnergy = p1.totalEnergy + p2.totalEnergy;
238 G4double x = p1.momentum.x() + p2.momentum.x();
239 G4double y = p1.momentum.y() + p2.momentum.y();
240 G4double z = p1.momentum.z() + p2.momentum.z();
241 G4double newMass = totEnergy*totEnergy - ( x*x + y*y + z*z );
242 if( newMass < 0.0 )
243 newMass = -1. * std::sqrt( -newMass );
244 else
245 newMass = std::sqrt( newMass );
246 G4ReactionProduct result;
247 result.SetMass( newMass );
248 result.SetMomentum( x, y, z );
249 result.SetTotalEnergy( totEnergy );
250 result.SetPositionInNucleus( 0.0, 0.0, 0.0 );
251 result.SetFormationTime(0.0);
252 result.HasInitialStateParton(false);
253 return result;
254 }
255
256 G4ReactionProduct operator-(
257 const G4ReactionProduct& p1, const G4ReactionProduct& p2 )
258 {
259 G4double totEnergy = p1.totalEnergy - p2.totalEnergy;
260 G4double x = p1.momentum.x() - p2.momentum.x();
261 G4double y = p1.momentum.y() - p2.momentum.y();
262 G4double z = p1.momentum.z() - p2.momentum.z();
263 G4double newMass = totEnergy*totEnergy - ( x*x + y*y + z*z );
264 if( newMass < 0.0 )
265 newMass = -1. * std::sqrt( -newMass );
266 else
267 newMass = std::sqrt( newMass );
268 G4ReactionProduct result;
269 result.SetMass( newMass );
270 result.SetMomentum( x, y, z );
271 result.SetTotalEnergy( totEnergy );
272 result.SetPositionInNucleus( 0.0, 0.0, 0.0 );
273 result.SetFormationTime(0.0);
274 result.HasInitialStateParton(false);
275 return result;
276 }
277 /* end of code */
278
Note: See TracBrowser for help on using the repository browser.