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