source: trunk/source/processes/hadronic/util/include/G4GHEKinematicsVector.hh @ 1326

Last change on this file since 1326 was 819, checked in by garnier, 16 years ago

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