Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4GHEKinematicsVector.hh
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 //
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 "G4ios.hh"
47 
48 #include <CLHEP/Units/PhysicalConstants.h>
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::fabs(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::fabs(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::fabs((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 L)
675  {
676  G4cout << "G4GHEKinematicsVector: "
677  << L << " " << 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
Generated on Sat May 25 2013 14:34:09 for Geant4 by   doxygen 1.8.4