G4GHEKinematicsVector.hh

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
//
// ------------------------------------------------------------
//      GEANT 4 class header file --- Copyright CERN 1998
//      CERN Geneva Switzerland
//
//      History: first implementation, based on object model of
//      2nd December 1995, G.Cosmo
//      ------------ G4GHEKinematicsVector utility class ------
//                   by Larry Felawka (TRIUMF), March 1997
//                     E-mail: felawka@alph04.triumf.ca
// ************************************************************
//-----------------------------------------------------------------------------

// Store, Retrieve and manipulate particle data.
// Based on "G4GHEVector" class of H. Fesefeldt.

#ifndef G4GHEKinematicsVector_h
#define G4GHEKinematicsVector_h 1

#include "G4ios.hh"

#include <CLHEP/Units/PhysicalConstants.h>

class G4GHEKinematicsVector
 {
 public:
  inline
   G4GHEKinematicsVector()
   {
     momentum.setX(  0.0 );
     momentum.setY(  0.0 );
     momentum.setZ(  0.0 );
     energy        = 0.0;
     kineticEnergy = 0.0;
     mass          = 0.0;
     charge        = 0.0;
     timeOfFlight  = 0.0;
     side          = 0;
     flag          = false;
     code          = 0;
     particleDef   = NULL;
   }

  ~G4GHEKinematicsVector() {}

  inline
   G4GHEKinematicsVector( const G4GHEKinematicsVector & p )
   {
     momentum.setX( p.momentum.x() );
     momentum.setY( p.momentum.y() );
     momentum.setZ( p.momentum.z() );
     energy        = p.energy;
     kineticEnergy = p.kineticEnergy;
     mass          = p.mass;
     charge        = p.charge;
     timeOfFlight  = p.timeOfFlight;
     side          = p.side;
     flag          = p.flag;
     code          = p.code;
     particleDef   = p.particleDef;
   }

  inline
   G4GHEKinematicsVector & operator = ( const G4GHEKinematicsVector & p )
   {
     if (this != &p)
     {
        momentum.setX( p.momentum.x() );
        momentum.setY( p.momentum.y() );
        momentum.setZ( p.momentum.z() );
        energy        = p.energy;
        kineticEnergy = p.kineticEnergy;
        mass          = p.mass;
        charge        = p.charge;
        timeOfFlight  = p.timeOfFlight;
        side          = p.side;
        flag          = p.flag;
        code          = p.code;
        particleDef   = p.particleDef;
     }
    return *this;
   }

  inline
   void SetMomentum( G4ParticleMomentum mom ) { momentum = mom; return; };

  inline
   void SetMomentumAndUpdate( G4ParticleMomentum mom )
   {
     momentum      = mom;
     energy        = std::sqrt(mass*mass + momentum.mag2());
     kineticEnergy = std::max(0.,energy - mass);
     return;
   }

  inline const
  G4ParticleMomentum GetMomentum() const { return momentum; }

  inline
   void SetMomentum( G4double x, G4double y, G4double z)
   { 
     momentum.setX( x );
     momentum.setY( y );
     momentum.setZ( z );
     return;
   } 

  inline
   void SetMomentumAndUpdate( G4double x, G4double y, G4double z )
   {
     momentum.setX( x );
     momentum.setY( y );
     momentum.setZ( z );
     energy        = std::sqrt(mass*mass + momentum.mag2());
     kineticEnergy = std::max(0.,energy-mass);
     return;
   }

  inline
   void SetMomentum( G4double x, G4double y )
   {
     momentum.setX( x );
     momentum.setY( y );
     return;
   }

  inline
   void SetMomentumAndUpdate( G4double x, G4double y )
   {
     momentum.setX( x );
     momentum.setY( y );
     energy = std::sqrt(mass*mass + momentum.mag2());
     kineticEnergy = std::max(0.,energy-mass);
     return;
   }

  inline
   void SetMomentum( G4double z )
   {
     momentum.setZ( z );
     return;
   }

  inline
   void SetMomentumAndUpdate( G4double z )
   {
     momentum.setZ( z );
     energy = std::sqrt(mass*mass + momentum.mag2());
     kineticEnergy = std::max(0.,energy-mass);
     return;
   }

  inline 
   void SetEnergy( G4double e ) { energy = e; return; }

  inline
   void SetEnergyAndUpdate( G4double e )
   {
     if (e <= mass)
       { 
         energy        = mass;
         kineticEnergy = 0.;
         momentum.setX(  0.);
         momentum.setY(  0.);
         momentum.setZ(  0.);
       }
     else
       {
         energy = e;
         kineticEnergy   = energy - mass;
         G4double momold = momentum.mag();
         G4double momnew = std::sqrt(energy*energy - mass*mass);
         if (momold == 0.)
           {
             G4double cost = 1.0- 2.0*G4UniformRand();
             G4double sint = std::sqrt(1. - cost*cost);
             G4double phi  = CLHEP::twopi* G4UniformRand();
             momentum.setX( momnew * sint * std::cos(phi));
             momentum.setY( momnew * sint * std::sin(phi));
             momentum.setZ( momnew * cost);
           }
         else
           {
             momnew /= momold;
             momentum.setX(momentum.x()*momnew);
             momentum.setY(momentum.y()*momnew);
             momentum.setZ(momentum.z()*momnew);
           }
       }    
     return;
   }

  inline 
   void SetKineticEnergy( G4double ekin ) { kineticEnergy = ekin; return; }

  inline 
   void SetKineticEnergyAndUpdate(G4double ekin) 
   {
     if (ekin <= 0.)
       { 
         energy        = mass;
         kineticEnergy = 0.;
         momentum.setX(  0.);
         momentum.setY(  0.);
         momentum.setZ(  0.);
       }
     else
       {
         energy = ekin + mass;
         kineticEnergy   = ekin;
         G4double momold = momentum.mag();
         G4double momnew = std::sqrt(energy*energy - mass*mass);
         if (momold == 0.)
           {
             G4double cost = 1.0-2.0*G4UniformRand();
             G4double sint = std::sqrt(1. - cost*cost);
             G4double phi  = CLHEP::twopi* G4UniformRand();
             momentum.setX( momnew * sint * std::cos(phi));
             momentum.setY( momnew * sint * std::sin(phi));
             momentum.setZ( momnew * cost);
           }
         else
           {
             momnew /= momold;
             momentum.setX(momentum.x()*momnew);
             momentum.setY(momentum.y()*momnew);
             momentum.setZ(momentum.z()*momnew);
           }
       }    
     return;
   }

  inline
  G4double GetEnergy() {return energy;}

  inline
  G4double GetKineticEnergy() {return kineticEnergy;}

  inline
  void SetMass( G4double mas ) { mass = mas; return; }

  inline
  void SetMassAndUpdate( G4double mas )
  {
    kineticEnergy = std::max(0., energy - mas);
    mass = mas;
    energy = kineticEnergy + mass;
    G4double momnew = std::sqrt(std::max(0., energy*energy - mass*mass));
    if ( momnew == 0.0) 
       {
         momentum.setX( 0.0 );
         momentum.setY( 0.0 );
         momentum.setZ( 0.0 );
       }
    else
       {
         G4double momold = momentum.mag();
         if (momold == 0.)
            { 
              G4double cost = 1.-2.*G4UniformRand();
              G4double sint = std::sqrt(1.-cost*cost);
              G4double phi  = CLHEP::twopi*G4UniformRand();
              momentum.setX( momnew*sint*std::cos(phi));
              momentum.setY( momnew*sint*std::sin(phi));
              momentum.setZ( momnew*cost);
            }
         else
            {
              momnew /= momold;
              momentum.setX( momentum.x()*momnew );
              momentum.setY( momentum.y()*momnew );
              momentum.setZ( momentum.z()*momnew );
            }
       }     
    return;
  }    

  inline
  G4double GetMass() { return mass; }

  inline
  void SetCharge( G4double c ) { charge = c; return; }

  inline
  G4double GetCharge() {return charge; }

  inline
  void SetTOF( G4double t ) { timeOfFlight = t; return; }

  inline
  G4double GetTOF() { return timeOfFlight; }

  inline
  void SetSide( G4int sid ) { side = sid; return; }

  inline
  G4int GetSide() { return side; }

  inline
  void setFlag( G4bool f ) { flag = f; return; }

  inline
  G4bool getFlag() { return flag; }

  inline
  void SetCode( G4int c ) { code = c; return; }

  inline
  void SetParticleDef( G4ParticleDefinition * c ) { particleDef = c; return; }

  inline
  G4int GetCode() { return code; } 

  inline
  G4ParticleDefinition * GetParticleDef() { return particleDef; } 

  inline
  void SetZero()
   {
     momentum.setX(  0.0 );
     momentum.setY(  0.0 );
     momentum.setZ(  0.0 );
     energy        = 0.0;
     kineticEnergy = 0.0;
     mass          = 0.0;
     charge        = 0.0;
     timeOfFlight  = 0.0;
     side          = 0;
     flag          = false;
     code          = 0;
     particleDef   = NULL;
   }

  inline
   void Add( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
   {
     momentum = p1.momentum + p2.momentum;
     energy = p1.energy + p2.energy;
     G4double b = energy*energy - momentum.mag2();
     if( b < 0 )
       mass = -1. * std::sqrt( -b );
     else
       mass = std::sqrt( b );
     kineticEnergy = std::max(0.,energy - mass);
     charge        = p1.charge + p2.charge;
     code          = p1.code   + p2.code;
     particleDef   = p1.particleDef;
   }

  inline
   void Sub( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
   {
     momentum = p1.momentum - p2.momentum;
     energy = p1.energy - p2.energy;
     G4double b = energy*energy - momentum.mag2();
     if( b < 0 )
       mass = -1. * std::sqrt( -b );
     else
       mass = std::sqrt( b );
     kineticEnergy = std::max(0.,energy - mass);
     charge        = p1.charge - p2.charge;
     code          = p1.code   - p2.code;
     particleDef   = p1.particleDef;
   }

  inline
   void Lor( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2 )
   {
     G4double a;
     a = ( p1.momentum.dot(p2.momentum)/(p2.energy+p2.mass) - p1.energy ) / p2.mass;
     momentum.setX( p1.momentum.x()+a*p2.momentum.x() );
     momentum.setY( p1.momentum.y()+a*p2.momentum.y() );
     momentum.setZ( p1.momentum.z()+a*p2.momentum.z() );
     energy = std::sqrt( sqr(p1.mass) + momentum.mag2() );
     mass = p1.mass;
     kineticEnergy = std::max(0.,energy - mass);
     timeOfFlight  = p1.timeOfFlight;
     side          = p1.side;
     flag          = p1.flag;
     code          = p1.code;
     particleDef   = p1.particleDef;
   }

  inline
   G4double CosAng( const G4GHEKinematicsVector & p )
   {
     G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
     if( a != 0.0 ) 
       {
         a = (momentum.x()*p.momentum.x() +
              momentum.y()*p.momentum.y() +
              momentum.z()*p.momentum.z()) / a;
         if( std::fabs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
       }
     return a;
   }
  inline
   G4double Ang(const G4GHEKinematicsVector & p )
   {
     G4double a = std::sqrt( momentum.mag2() * p.momentum.mag2() );
     if( a != 0.0 ) 
       {
         a = (momentum.x()*p.momentum.x() +
              momentum.y()*p.momentum.y() +
              momentum.z()*p.momentum.z()) / a;
         if( std::fabs(a) > 1.0 ) a<0.0 ? a=-1.0 : a=1.0;
       }
     return std::acos(a);
   }      

  inline
   G4double Dot4( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     return (   p1.energy       * p2.energy
              - p1.momentum.x() * p2.momentum.x()
              - p1.momentum.y() * p2.momentum.y()
              - p1.momentum.z() * p2.momentum.z() );
   } 
  
  inline 
   G4double Impu( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     return ( - sqr( p1.energy      - p2.energy)  
              + sqr(p1.momentum.x() - p2.momentum.x())
              + sqr(p1.momentum.y() - p2.momentum.y()) 
              + sqr(p1.momentum.z() - p2.momentum.z()) );
   }

  inline 
   void Add3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     momentum.setX( p1.momentum.x() + p2.momentum.x());
     momentum.setY( p1.momentum.y() + p2.momentum.y());
     momentum.setZ( p1.momentum.z() + p2.momentum.z());    
     return;
   }

  inline
   void Sub3( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     momentum.setX( p1.momentum.x() - p2.momentum.x());
     momentum.setY( p1.momentum.y() - p2.momentum.y());
     momentum.setZ( p1.momentum.z() - p2.momentum.z());
     return;
   }

  inline
   void Cross( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     G4double px, py, pz;
     px = p1.momentum.y() * p2.momentum.z() - p1.momentum.z() * p2.momentum.y();
     py = p1.momentum.z() * p2.momentum.x() - p1.momentum.x() * p2.momentum.z();
     pz = p1.momentum.x() * p2.momentum.y() - p1.momentum.y() * p2.momentum.x();
     momentum.setX( px );
     momentum.setY( py );
     momentum.setZ( pz ); 
     return;
   }  

  inline 
   G4double Dot( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     return (   p1.momentum.x() * p2.momentum.x()
              + p1.momentum.y() * p2.momentum.y()
              + p1.momentum.z() * p2.momentum.z() );
   }

  inline
   void Smul( const G4GHEKinematicsVector & p, G4double h)
   {
     momentum.setX( h * p.momentum.x());
     momentum.setY( h * p.momentum.y());
     momentum.setZ( h * p.momentum.z());
     return;
   }   

  inline
   void SmulAndUpdate( const G4GHEKinematicsVector & p, G4double h)
   {
     momentum.setX( h * p.momentum.x());
     momentum.setY( h * p.momentum.y());
     momentum.setZ( h * p.momentum.z());
     mass          = p.mass;
     energy        = std::sqrt(momentum.mag2() + mass*mass);
     kineticEnergy = energy - mass;
     charge        = p.charge;
     timeOfFlight  = p.timeOfFlight;
     side          = p.side;
     flag          = p.flag;
     code          = p.code;
     particleDef   = p.particleDef;
     return;
   }
 
  inline
   void Norz( const G4GHEKinematicsVector & p )
   {
     G4double a =   p.momentum.mag2();
     if (a > 0.0) a = 1./std::sqrt(a);
     momentum.setX( a * p.momentum.x() );
     momentum.setY( a * p.momentum.y() );
     momentum.setZ( a * p.momentum.z() );
     mass          = p.mass;
     energy        = std::sqrt(momentum.mag2() + mass*mass);
     kineticEnergy = energy - mass;
     charge        = p.charge;
     timeOfFlight  = p.timeOfFlight;
     side          = p.side;
     flag          = p.flag;
     code          = p.code; 
     particleDef   = p.particleDef; 
     return;
   }

  inline 
   G4double Length()
   {
     return  momentum.mag() ;
   }

  inline
   void Exch( G4GHEKinematicsVector & p1)
   {
     G4GHEKinematicsVector mx = *this;
//     mx.momentum.SetX( momentum.x());
//     mx.momentum.SetY( momentum.y());
//     mx.momentum.SetZ( momentum.z());
//     mx.energy        = energy;
//     mx.kineticEnergy = kineticEnergy;
//     mx.mass          = mass;
//     mx.charge        = charge;
//     mx.timeOfFlight  = timeOfFlight;
//     mx.side          = side;
//     mx.flag          = flag;
//     mx.code          = code; 
//     momentum.setX( p1.momentum.x());
//     momentum.setY( p1.momentum.y());
//     momentum.setZ( p1.momentum.z());
//     energy        = p1.energy;
//     kineticEnergy = p1.kineticEnergy;
//     mass          = p1.mass;
//     charge        = p1.charge;
//     timeOfFlight  = p1.timeOfFlight;
//     side          = p1.side
//     flag          = p1.flag;
//     code          = p1.code;
     *this = p1;
     p1 = mx;
     return; 
   }      

  inline
   void Defs1( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2)
   {
     G4double pt2 = sqr(p1.momentum.x()) + sqr(p1.momentum.y());
     if (pt2 > 0.0)
        {
          G4double ph, px, py, pz;
          G4double cost  = p2.momentum.z()/p2.momentum.mag(); 
          G4double sint  = 0.5 * (  std::sqrt(std::fabs((1.-cost)*(1.+cost))) 
                                  + std::sqrt(pt2)/p2.momentum.mag());
          (p2.momentum.y() < 0.) ? ph = 1.5*CLHEP::pi : ph = CLHEP::halfpi;
          if( p2.momentum.x() != 0.0) 
             ph = std::atan2(p2.momentum.y(),p2.momentum.x());             
          px =   cost*std::cos(ph)*p1.momentum.x() - std::sin(ph)*p1.momentum.y()
               + sint*std::cos(ph)*p1.momentum.z();
          py =   cost*std::sin(ph)*p1.momentum.x() + std::cos(ph)*p1.momentum.y()
               + sint*std::sin(ph)*p1.momentum.z();
          pz = - sint        *p1.momentum.x() 
               + cost        *p1.momentum.z();
          momentum.setX( px ); 
          momentum.setY( py );
          momentum.setZ( pz );     
        }
     else
        {
          momentum = p1.momentum;             
    } 
   }

  inline 
   void Defs( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & p2,
                    G4GHEKinematicsVector & my,       G4GHEKinematicsVector & mz )
   {
     my = p1;
     mz = p2;
     momentum.setX(   my.momentum.y()*mz.momentum.z()
                    - my.momentum.z()*mz.momentum.y());
     momentum.setY(   my.momentum.z()*mz.momentum.x()
                    - my.momentum.x()*mz.momentum.z());
     momentum.setZ(   my.momentum.x()*mz.momentum.y()
                    - my.momentum.y()*mz.momentum.x());
     my.momentum.setX(   mz.momentum.y()*momentum.z()
                       - mz.momentum.z()*momentum.y());
     my.momentum.setY(   mz.momentum.z()*momentum.x()
                       - mz.momentum.x()*momentum.z());
     my.momentum.setZ(   mz.momentum.x()*momentum.y()
                       - mz.momentum.y()*momentum.x());
     G4double pp;
     pp = momentum.mag();
     if (pp > 0.)
        {
          pp = 1./pp; 
          momentum.setX( momentum.x()*pp );
          momentum.setY( momentum.y()*pp );
          momentum.setZ( momentum.z()*pp );
        }
     pp = my.momentum.mag();
     if (pp > 0.)
        {
          pp = 1./pp;
          my.momentum.setX( my.momentum.x()*pp );
          my.momentum.setY( my.momentum.y()*pp );
          my.momentum.setZ( my.momentum.z()*pp );
        }
     pp = mz.momentum.mag();
     if (pp > 0.)
        {
          pp = 1./pp;
          mz.momentum.setX( mz.momentum.x()*pp );
          mz.momentum.setY( mz.momentum.y()*pp );
          mz.momentum.setZ( mz.momentum.z()*pp );
        }
     return; 
   }  
             
  inline
   void Trac( const G4GHEKinematicsVector & p1, const G4GHEKinematicsVector & mx,
              const G4GHEKinematicsVector & my, const G4GHEKinematicsVector & mz)
   {
     double px, py, pz;
     px =   mx.momentum.x()*p1.momentum.x()
          + mx.momentum.y()*p1.momentum.y()
          + mx.momentum.z()*p1.momentum.z();
     py =   my.momentum.x()*p1.momentum.x()
          + my.momentum.y()*p1.momentum.y()
          + my.momentum.z()*p1.momentum.z();
     pz =   mz.momentum.x()*p1.momentum.x()
          + mz.momentum.y()*p1.momentum.y()
          + mz.momentum.z()*p1.momentum.z();
     momentum.setX( px );
     momentum.setY( py );
     momentum.setZ( pz );
     return;
   }                 

  inline
   void Print( G4int L)
   {
     G4cout << "G4GHEKinematicsVector: " 
          << L << " " << momentum.x() << " " <<  momentum.y() << " " <<  momentum.z() << " "
          << energy << " " << kineticEnergy << " " << mass << " " << charge << " " 
          << timeOfFlight << " " << side << " " << flag << " " << code << particleDef << G4endl;
     return;                         
   }

  G4ParticleMomentum momentum;
  G4double energy;
  G4double kineticEnergy;
  G4double mass;
  G4double charge;
  G4double timeOfFlight;
  G4int side;
  G4bool flag;
  G4int code;
  G4ParticleDefinition * particleDef;
};

#endif