ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GHEKinematicsVector.hh
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4GHEKinematicsVector.hh
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 // ------------------------------------------------------------
29 // GEANT 4 class header file --- Copyright CERN 1998
30 // CERN Geneva Switzerland
31 //
32 // History: first implementation, based on object model of
33 // 2nd December 1995, G.Cosmo
34 // ------------ G4GHEKinematicsVector utility class ------
35 // by Larry Felawka (TRIUMF), March 1997
36 // E-mail: felawka@alph04.triumf.ca
37 // ************************************************************
38 //-----------------------------------------------------------------------------
39 
40 // Store, Retrieve and manipulate particle data.
41 // Based on "G4GHEVector" class of H. Fesefeldt.
42 
43 #ifndef G4GHEKinematicsVector_h
44 #define G4GHEKinematicsVector_h 1
45 
46 #include <CLHEP/Units/PhysicalConstants.h>
47 #include "G4ios.hh"
48 #include "G4ParticleMomentum.hh"
49 
50 class G4GHEKinematicsVector
51  {
52  public:
53  inline
54  G4GHEKinematicsVector()
55  {
56  momentum.setX( 0.0 );
57  momentum.setY( 0.0 );
58  momentum.setZ( 0.0 );
59  energy = 0.0;
60  kineticEnergy = 0.0;
61  mass = 0.0;
62  charge = 0.0;
63  timeOfFlight = 0.0;
64  side = 0;
65  flag = false;
66  code = 0;
67  particleDef = NULL;
68  }
69 
70  ~G4GHEKinematicsVector() {}
71 
72  inline
73  G4GHEKinematicsVector( const G4GHEKinematicsVector & p )
74  {
75  momentum.setX( p.momentum.x() );
76  momentum.setY( p.momentum.y() );
77  momentum.setZ( p.momentum.z() );
78  energy = p.energy;
79  kineticEnergy = p.kineticEnergy;
80  mass = p.mass;
81  charge = p.charge;
82  timeOfFlight = p.timeOfFlight;
83  side = p.side;
84  flag = p.flag;
85  code = p.code;
86  particleDef = p.particleDef;
87  }
88 
89  inline
90  G4GHEKinematicsVector & operator = ( const G4GHEKinematicsVector & p )
91  {
92  if (this != &p)
93  {
94  momentum.setX( p.momentum.x() );
95  momentum.setY( p.momentum.y() );
96  momentum.setZ( p.momentum.z() );
97  energy = p.energy;
98  kineticEnergy = p.kineticEnergy;
99  mass = p.mass;
100  charge = p.charge;
101  timeOfFlight = p.timeOfFlight;
102  side = p.side;
103  flag = p.flag;
104  code = p.code;
105  particleDef = p.particleDef;
106  }
107  return *this;
108  }
109 
110  inline
111  void SetMomentum( G4ParticleMomentum mom ) { momentum = mom; return; };
112 
113  inline
114  void SetMomentumAndUpdate( G4ParticleMomentum mom )
115  {
116  momentum = mom;
117  energy = std::sqrt(mass*mass + momentum.mag2());
118  kineticEnergy = std::max(0.,energy - mass);
119  return;
120  }
121 
122  inline const
123  G4ParticleMomentum GetMomentum() const { return momentum; }
124 
125  inline
126  void SetMomentum( G4double x, G4double y, G4double z)
127  {
128  momentum.setX( x );
129  momentum.setY( y );
130  momentum.setZ( z );
131  return;
132  }
133 
134  inline
135  void SetMomentumAndUpdate( G4double x, G4double y, G4double z )
136  {
137  momentum.setX( x );
138  momentum.setY( y );
139  momentum.setZ( z );
140  energy = std::sqrt(mass*mass + momentum.mag2());
141  kineticEnergy = std::max(0.,energy-mass);
142  return;
143  }
144 
145  inline
146  void SetMomentum( G4double x, G4double y )
147  {
148  momentum.setX( x );
149  momentum.setY( y );
150  return;
151  }
152 
153  inline
154  void SetMomentumAndUpdate( G4double x, G4double y )
155  {
156  momentum.setX( x );
157  momentum.setY( y );
158  energy = std::sqrt(mass*mass + momentum.mag2());
159  kineticEnergy = std::max(0.,energy-mass);
160  return;
161  }
162 
163  inline
164  void SetMomentum( G4double z )
165  {
166  momentum.setZ( z );
167  return;
168  }
169 
170  inline
171  void SetMomentumAndUpdate( G4double z )
172  {
173  momentum.setZ( z );
174  energy = std::sqrt(mass*mass + momentum.mag2());
175  kineticEnergy = std::max(0.,energy-mass);
176  return;
177  }
178 
179  inline
180  void SetEnergy( G4double e ) { energy = e; return; }
181 
182  inline
183  void SetEnergyAndUpdate( G4double e )
184  {
185  if (e <= mass)
186  {
187  energy = mass;
188  kineticEnergy = 0.;
189  momentum.setX( 0.);
190  momentum.setY( 0.);
191  momentum.setZ( 0.);
192  }
193  else
194  {
195  energy = e;
196  kineticEnergy = energy - mass;
197  G4double momold = momentum.mag();
198  G4double momnew = std::sqrt(energy*energy - mass*mass);
199  if (momold == 0.)
200  {
201  G4double cost = 1.0- 2.0*G4UniformRand();
202  G4double sint = std::sqrt(1. - cost*cost);
203  G4double phi = CLHEP::twopi* G4UniformRand();
204  momentum.setX( momnew * sint * std::cos(phi));
205  momentum.setY( momnew * sint * std::sin(phi));
206  momentum.setZ( momnew * cost);
207  }
208  else
209  {
210  momnew /= momold;
211  momentum.setX(momentum.x()*momnew);
212  momentum.setY(momentum.y()*momnew);
213  momentum.setZ(momentum.z()*momnew);
214  }
215  }
216  return;
217  }
218 
219  inline
220  void SetKineticEnergy( G4double ekin ) { kineticEnergy = ekin; return; }
221 
222  inline
223  void SetKineticEnergyAndUpdate(G4double ekin)
224  {
225  if (ekin <= 0.)
226  {
227  energy = mass;
228  kineticEnergy = 0.;
229  momentum.setX( 0.);
230  momentum.setY( 0.);
231  momentum.setZ( 0.);
232  }
233  else
234  {
235  energy = ekin + mass;
236  kineticEnergy = ekin;
237  G4double momold = momentum.mag();
238  G4double momnew = std::sqrt(energy*energy - mass*mass);
239  if (momold == 0.)
240  {
241  G4double cost = 1.0-2.0*G4UniformRand();
242  G4double sint = std::sqrt(1. - cost*cost);
243  G4double phi = CLHEP::twopi* G4UniformRand();
244  momentum.setX( momnew * sint * std::cos(phi));
245  momentum.setY( momnew * sint * std::sin(phi));
246  momentum.setZ( momnew * cost);
247  }
248  else
249  {
250  momnew /= momold;
251  momentum.setX(momentum.x()*momnew);
252  momentum.setY(momentum.y()*momnew);
253  momentum.setZ(momentum.z()*momnew);
254  }
255  }
256  return;
257  }
258 
259  inline
260  G4double GetEnergy() {return energy;}
261 
262  inline
263  G4double GetKineticEnergy() {return kineticEnergy;}
264 
265  inline
266  void SetMass( G4double mas ) { mass = mas; return; }
267 
268  inline
269  void SetMassAndUpdate( G4double mas )
270  {
271  kineticEnergy = std::max(0., energy - mas);
272  mass = mas;
273  energy = kineticEnergy + mass;
274  G4double momnew = std::sqrt(std::max(0., energy*energy - mass*mass));
275  if ( momnew == 0.0)
276  {
277  momentum.setX( 0.0 );
278  momentum.setY( 0.0 );
279  momentum.setZ( 0.0 );
280  }
281  else
282  {
283  G4double momold = momentum.mag();
284  if (momold == 0.)
285  {
286  G4double cost = 1.-2.*G4UniformRand();
287  G4double sint = std::sqrt(1.-cost*cost);
288  G4double phi = CLHEP::twopi*G4UniformRand();
289  momentum.setX( momnew*sint*std::cos(phi));
290  momentum.setY( momnew*sint*std::sin(phi));
291  momentum.setZ( momnew*cost);
292  }
293  else
294  {
295  momnew /= momold;
296  momentum.setX( momentum.x()*momnew );
297  momentum.setY( momentum.y()*momnew );
298  momentum.setZ( momentum.z()*momnew );
299  }
300  }
301  return;
302  }
303 
304  inline
305  G4double GetMass() { return mass; }
306 
307  inline
308  void SetCharge( G4double c ) { charge = c; return; }
309 
310  inline
311  G4double GetCharge() {return charge; }
312 
313  inline
314  void SetTOF( G4double t ) { timeOfFlight = t; return; }
315 
316  inline
317  G4double GetTOF() { return timeOfFlight; }
318 
319  inline
320  void SetSide( G4int sid ) { side = sid; return; }
321 
322  inline
323  G4int GetSide() { return side; }
324 
325  inline
326  void setFlag( G4bool f ) { flag = f; return; }
327 
328  inline
329  G4bool getFlag() { return flag; }
330 
331  inline
332  void SetCode( G4int c ) { code = c; return; }
333 
334  inline
335  void SetParticleDef( G4ParticleDefinition * c ) { particleDef = c; return; }
336 
337  inline
338  G4int GetCode() { return code; }
339 
340  inline
341  G4ParticleDefinition * GetParticleDef() { return particleDef; }
342 
343  inline
344  void SetZero()
345  {
346  momentum.setX( 0.0 );
347  momentum.setY( 0.0 );
348  momentum.setZ( 0.0 );
349  energy = 0.0;
350  kineticEnergy = 0.0;
351  mass = 0.0;
352  charge = 0.0;
353  timeOfFlight = 0.0;
354  side = 0;
355  flag = false;
356  code = 0;
357  particleDef = NULL;
358  }
359 
360  inline
361  void Add( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
362  {
363  momentum = p1.momentum + p2.momentum;
364  energy = p1.energy + p2.energy;
365  G4double b = energy*energy - momentum.mag2();
366  if( b < 0 )
367  mass = -1. * std::sqrt( -b );
368  else
369  mass = std::sqrt( b );
370  kineticEnergy = std::max(0.,energy - mass);
371  charge = p1.charge + p2.charge;
372  code = p1.code + p2.code;
373  particleDef = p1.particleDef;
374  }
375 
376  inline
377  void Sub( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
378  {
379  momentum = p1.momentum - p2.momentum;
380  energy = p1.energy - p2.energy;
381  G4double b = energy*energy - momentum.mag2();
382  if( b < 0 )
383  mass = -1. * std::sqrt( -b );
384  else
385  mass = std::sqrt( b );
386  kineticEnergy = std::max(0.,energy - mass);
387  charge = p1.charge - p2.charge;
388  code = p1.code - p2.code;
389  particleDef = p1.particleDef;
390  }
391 
392  inline
393  void Lor( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
394  {
395  G4double a;
396  a = ( p1.momentum.dot(p2.momentum)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
397  momentum.setX( p1.momentum.x()+a*p2.momentum.x() );
398  momentum.setY( p1.momentum.y()+a*p2.momentum.y() );
399  momentum.setZ( p1.momentum.z()+a*p2.momentum.z() );
400  energy = std::sqrt( sqr(p1.mass) + momentum.mag2() );
401  mass = p1.mass;
402  kineticEnergy = std::max(0.,energy - mass);
403  timeOfFlight = p1.timeOfFlight;
404  side = p1.side;
405  flag = p1.flag;
406  code = p1.code;
407  particleDef = p1.particleDef;
408  }
409 
410  inline
411  G4double CosAng( const G4GHEKinematicsVector & p )
412  {
413  G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
414  if( a != 0.0 )
415  {
416  a = (momentum.x()*p.momentum.x() +
417  momentum.y()*p.momentum.y() +
418  momentum.z()*p.momentum.z()) / a;
419  if( std::abs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
420  }
421  return a;
422  }
423  inline
424  G4double Ang(const G4GHEKinematicsVector & p )
425  {
426  G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
427  if( a != 0.0 )
428  {
429  a = (momentum.x()*p.momentum.x() +
430  momentum.y()*p.momentum.y() +
431  momentum.z()*p.momentum.z()) / a;
432  if( std::abs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
433  }
434  return std::acos(a);
435  }
436 
437  inline
438  G4double Dot4( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
439  {
440  return ( p1.energy * p2.energy
441  - p1.momentum.x() * p2.momentum.x()
442  - p1.momentum.y() * p2.momentum.y()
443  - p1.momentum.z() * p2.momentum.z() );
444  }
445 
446  inline
447  G4double Impu( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
448  {
449  return ( - sqr( p1.energy - p2.energy)
450  + sqr(p1.momentum.x() - p2.momentum.x())
451  + sqr(p1.momentum.y() - p2.momentum.y())
452  + sqr(p1.momentum.z() - p2.momentum.z()) );
453  }
454 
455  inline
456  void Add3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
457  {
458  momentum.setX( p1.momentum.x() + p2.momentum.x());
459  momentum.setY( p1.momentum.y() + p2.momentum.y());
460  momentum.setZ( p1.momentum.z() + p2.momentum.z());
461  return;
462  }
463 
464  inline
465  void Sub3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
466  {
467  momentum.setX( p1.momentum.x() - p2.momentum.x());
468  momentum.setY( p1.momentum.y() - p2.momentum.y());
469  momentum.setZ( p1.momentum.z() - p2.momentum.z());
470  return;
471  }
472 
473  inline
474  void Cross( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
475  {
476  G4double px, py, pz;
477  px = p1.momentum.y() * p2.momentum.z() - p1.momentum.z() * p2.momentum.y();
478  py = p1.momentum.z() * p2.momentum.x() - p1.momentum.x() * p2.momentum.z();
479  pz = p1.momentum.x() * p2.momentum.y() - p1.momentum.y() * p2.momentum.x();
480  momentum.setX( px );
481  momentum.setY( py );
482  momentum.setZ( pz );
483  return;
484  }
485 
486  inline
487  G4double Dot( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
488  {
489  return ( p1.momentum.x() * p2.momentum.x()
490  + p1.momentum.y() * p2.momentum.y()
491  + p1.momentum.z() * p2.momentum.z() );
492  }
493 
494  inline
495  void Smul( const G4GHEKinematicsVector & p, G4double h)
496  {
497  momentum.setX( h * p.momentum.x());
498  momentum.setY( h * p.momentum.y());
499  momentum.setZ( h * p.momentum.z());
500  return;
501  }
502 
503  inline
504  void SmulAndUpdate( const G4GHEKinematicsVector & p, G4double h)
505  {
506  momentum.setX( h * p.momentum.x());
507  momentum.setY( h * p.momentum.y());
508  momentum.setZ( h * p.momentum.z());
509  mass = p.mass;
510  energy = std::sqrt(momentum.mag2() + mass*mass);
511  kineticEnergy = energy - mass;
512  charge = p.charge;
513  timeOfFlight = p.timeOfFlight;
514  side = p.side;
515  flag = p.flag;
516  code = p.code;
517  particleDef = p.particleDef;
518  return;
519  }
520 
521  inline
522  void Norz( const G4GHEKinematicsVector & p )
523  {
524  G4double a = p.momentum.mag2();
525  if (a > 0.0) a = 1./std::sqrt(a);
526  momentum.setX( a * p.momentum.x() );
527  momentum.setY( a * p.momentum.y() );
528  momentum.setZ( a * p.momentum.z() );
529  mass = p.mass;
530  energy = std::sqrt(momentum.mag2() + mass*mass);
531  kineticEnergy = energy - mass;
532  charge = p.charge;
533  timeOfFlight = p.timeOfFlight;
534  side = p.side;
535  flag = p.flag;
536  code = p.code;
537  particleDef = p.particleDef;
538  return;
539  }
540 
541  inline
542  G4double Length()
543  {
544  return momentum.mag() ;
545  }
546 
547  inline
548  void Exch( G4GHEKinematicsVector & p1)
549  {
550  G4GHEKinematicsVector mx = *this;
551 // mx.momentum.SetX( momentum.x());
552 // mx.momentum.SetY( momentum.y());
553 // mx.momentum.SetZ( momentum.z());
554 // mx.energy = energy;
555 // mx.kineticEnergy = kineticEnergy;
556 // mx.mass = mass;
557 // mx.charge = charge;
558 // mx.timeOfFlight = timeOfFlight;
559 // mx.side = side;
560 // mx.flag = flag;
561 // mx.code = code;
562 // momentum.setX( p1.momentum.x());
563 // momentum.setY( p1.momentum.y());
564 // momentum.setZ( p1.momentum.z());
565 // energy = p1.energy;
566 // kineticEnergy = p1.kineticEnergy;
567 // mass = p1.mass;
568 // charge = p1.charge;
569 // timeOfFlight = p1.timeOfFlight;
570 // side = p1.side
571 // flag = p1.flag;
572 // code = p1.code;
573  *this = p1;
574  p1 = mx;
575  return;
576  }
577 
578  inline
579  void Defs1( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
580  {
581  G4double pt2 = sqr(p1.momentum.x()) + sqr(p1.momentum.y());
582  if (pt2 > 0.0)
583  {
584  G4double ph, px, py, pz;
585  G4double cost = p2.momentum.z()/p2.momentum.mag();
586  G4double sint = 0.5 * ( std::sqrt(std::abs((1.-cost)*(1.+cost)))
587  + std::sqrt(pt2)/p2.momentum.mag());
588  (p2.momentum.y() < 0.) ? ph = 1.5*CLHEP::pi : ph = CLHEP::halfpi;
589  if( p2.momentum.x() != 0.0)
590  ph = std::atan2(p2.momentum.y(),p2.momentum.x());
591  px = cost*std::cos(ph)*p1.momentum.x() - std::sin(ph)*p1.momentum.y()
592  + sint*std::cos(ph)*p1.momentum.z();
593  py = cost*std::sin(ph)*p1.momentum.x() + std::cos(ph)*p1.momentum.y()
594  + sint*std::sin(ph)*p1.momentum.z();
595  pz = - sint *p1.momentum.x()
596  + cost *p1.momentum.z();
597  momentum.setX( px );
598  momentum.setY( py );
599  momentum.setZ( pz );
600  }
601  else
602  {
603  momentum = p1.momentum;
604  }
605  }
606 
607  inline
608  void Defs( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2,
609  G4GHEKinematicsVector & my, G4GHEKinematicsVector & mz )
610  {
611  my = p1;
612  mz = p2;
613  momentum.setX( my.momentum.y()*mz.momentum.z()
614  - my.momentum.z()*mz.momentum.y());
615  momentum.setY( my.momentum.z()*mz.momentum.x()
616  - my.momentum.x()*mz.momentum.z());
617  momentum.setZ( my.momentum.x()*mz.momentum.y()
618  - my.momentum.y()*mz.momentum.x());
619  my.momentum.setX( mz.momentum.y()*momentum.z()
620  - mz.momentum.z()*momentum.y());
621  my.momentum.setY( mz.momentum.z()*momentum.x()
622  - mz.momentum.x()*momentum.z());
623  my.momentum.setZ( mz.momentum.x()*momentum.y()
624  - mz.momentum.y()*momentum.x());
625  G4double pp;
626  pp = momentum.mag();
627  if (pp > 0.)
628  {
629  pp = 1./pp;
630  momentum.setX( momentum.x()*pp );
631  momentum.setY( momentum.y()*pp );
632  momentum.setZ( momentum.z()*pp );
633  }
634  pp = my.momentum.mag();
635  if (pp > 0.)
636  {
637  pp = 1./pp;
638  my.momentum.setX( my.momentum.x()*pp );
639  my.momentum.setY( my.momentum.y()*pp );
640  my.momentum.setZ( my.momentum.z()*pp );
641  }
642  pp = mz.momentum.mag();
643  if (pp > 0.)
644  {
645  pp = 1./pp;
646  mz.momentum.setX( mz.momentum.x()*pp );
647  mz.momentum.setY( mz.momentum.y()*pp );
648  mz.momentum.setZ( mz.momentum.z()*pp );
649  }
650  return;
651  }
652 
653  inline
654  void Trac( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & mx,
655  const G4GHEKinematicsVector & my, const G4GHEKinematicsVector & mz)
656  {
657  double px, py, pz;
658  px = mx.momentum.x()*p1.momentum.x()
659  + mx.momentum.y()*p1.momentum.y()
660  + mx.momentum.z()*p1.momentum.z();
661  py = my.momentum.x()*p1.momentum.x()
662  + my.momentum.y()*p1.momentum.y()
663  + my.momentum.z()*p1.momentum.z();
664  pz = mz.momentum.x()*p1.momentum.x()
665  + mz.momentum.y()*p1.momentum.y()
666  + mz.momentum.z()*p1.momentum.z();
667  momentum.setX( px );
668  momentum.setY( py );
669  momentum.setZ( pz );
670  return;
671  }
672 
673  inline
674  void Print( G4int LLL)
675  {
676  G4cout << "G4GHEKinematicsVector: "
677  << LLL << " " << momentum.x() << " " << momentum.y() << " " << momentum.z() << " "
678  << energy << " " << kineticEnergy << " " << mass << " " << charge << " "
679  << timeOfFlight << " " << side << " " << flag << " " << code << particleDef << G4endl;
680  return;
681  }
682 
683  G4ParticleMomentum momentum;
684  G4double energy;
685  G4double kineticEnergy;
686  G4double mass;
687  G4double charge;
688  G4double timeOfFlight;
689  G4int side;
690  G4bool flag;
691  G4int code;
692  G4ParticleDefinition * particleDef;
693 };
694 
695 #endif