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