Ex03dParticle.h

Go to the documentation of this file.

#ifndef EX03_PARTICLE_H
#define EX03_PARTICLE_H
 
//------------------------------------------------
// The Virtual Monte Carlo examples
// Copyright (C) 2014 - 2018 Ivana Hrivnacova
// All rights reserved.
//
// For the licensing terms see geant4_vmc/LICENSE.
// Contact: root-vmc@cern.ch
//-------------------------------------------------
 
/// \file  Ex03dParticle.h
/// \brief Definition of the Ex03dParticle class
///
/// Simplified version of ROOT's TParticle to store in RNTuple
///
/// \author Radoslaw Karabowicz; GSI
 
#include "TLorentzVector.h"
 
/// \ingroup E03
/// \brief Replacement of the TParticle
///
/// A variant of the TParticle class
/// that can be stored in the RNTuple.
///
 
/// \date 07/07/2026
/// \author Radoslaw Karabowicz, GSI
 
class Ex03dParticle
{
 public:
  // ****** constructors and destructor
  Ex03dParticle();
 
  Ex03dParticle(Int_t pdg, Int_t status, Int_t mother1, Int_t mother2,
    Int_t daughter1, Int_t daughter2, Double_t px, Double_t py, Double_t pz,
    Double_t etot, Double_t vx, Double_t vy, Double_t vz, Double_t time);
 
  Ex03dParticle(Int_t pdg, Int_t status, Int_t mother1, Int_t mother2,
    Int_t daughter1, Int_t daughter2, const TLorentzVector& p,
    const TLorentzVector& v);
 
  Ex03dParticle(const Ex03dParticle& part);
 
  ~Ex03dParticle();
 
  // methods
  Ex03dParticle& operator=(const Ex03dParticle&);
 
  Double_t Ek() const { return fE - fCalcMass; }
  Int_t GetStatusCode() const { return fStatusCode; }
  Int_t GetPdgCode() const { return fPdgCode; }
  Int_t GetFirstMother() const { return fMother[0]; }
  Int_t GetMother(Int_t i) const { return fMother[i]; }
  Int_t GetSecondMother() const { return fMother[1]; }
  Bool_t IsPrimary() const
  {
    return fMother[0] > -1 ? kFALSE : kTRUE;
  } // Is this particle primary one?
  Int_t GetFirstDaughter() const { return fDaughter[0]; }
  Int_t GetDaughter(Int_t i) const { return fDaughter[i]; }
  Int_t GetLastDaughter() const { return fDaughter[1]; }
  Double_t GetCalcMass() const { return fCalcMass; }
  Double_t GetMass() const;
  Int_t GetNDaughters() const
  {
    return fDaughter[1] > 0 ? fDaughter[1] - fDaughter[0] + 1 : 0;
  }
  Float_t GetWeight() const { return fWeight; }
  void GetPolarisation(TVector3& v) const;
  //  TParticlePDG*  GetPDG          (Int_t mode = 0) const;
  Int_t Beauty() const;
  Int_t Charm() const;
  Int_t Strangeness() const;
  Double_t PhiX() const
  {
    return TMath::Pi() + TMath::ATan2(-fPz, -fPy);
  } // note that PhiX() returns an angle between 0 and 2pi
  Double_t PhiY() const
  {
    return TMath::Pi() + TMath::ATan2(-fPx, -fPz);
  } // note that PhiY() returns an angle between 0 and 2pi
  Double_t PhiZ() const
  {
    return TMath::Pi() + TMath::ATan2(-fPy, -fPx);
  } // note that PhiZ() returns an angle between 0 and 2pi
  Double_t ThetaX() const
  {
    return (fPx == 0) ? TMath::PiOver2() : TMath::ACos(fPx / P());
  }
  Double_t ThetaY() const
  {
    return (fPy == 0) ? TMath::PiOver2() : TMath::ACos(fPy / P());
  }
  Double_t ThetaZ() const
  {
    return (fPz == 0) ? TMath::PiOver2() : TMath::ACos(fPz / P());
  }
  Double_t GetPolarTheta() const { return fPolarTheta; }
  Double_t GetPolarPhi() const { return fPolarPhi; }
  void GetPolarisation(Double_t& theta, Double_t& phi) const
  {
    theta = fPolarTheta;
    phi = fPolarPhi;
  }
  void SetPolarTheta(Double_t theta) { fPolarTheta = theta; }
  void SetPolarPhi(Double_t phi) { fPolarPhi = phi; }
  void SetPolarisation(Double_t theta, Double_t phi)
  {
    fPolarTheta = theta;
    fPolarPhi = phi;
  }
  void Momentum(TLorentzVector& v) const { v.SetPxPyPzE(fPx, fPy, fPz, fE); }
  void ProductionVertex(TLorentzVector& v) const
  {
    v.SetXYZT(fVx, fVy, fVz, fVt);
  }
 
  Double_t Theta(
    const Ex03dParticle& p) // Returns the angle between momenta of particles
  {
    Double_t v = P() * p.P();
    if (v == 0)
      v = 1;
    else
      v = (fPx * p.Px() + fPy * p.Py() + fPz * p.Pz()) / v;
    if (v > 1)
      v = 1;
    else if (v < -1)
      v = -1; // just a precaution
    return TMath::ACos(v);
  }
 
  // ****** redefine several most oftenly used methods of LORENTZ_VECTOR
 
  Double_t Vx() const { return fVx; }
  Double_t Vy() const { return fVy; }
  Double_t Vz() const { return fVz; }
  Double_t T() const { return fVt; }
  Double_t R() const
  {
    return TMath::Sqrt(fVx * fVx + fVy * fVy);
  } // Radius of production vertex in cylindrical system
  Double_t Rho() const
  {
    return TMath::Sqrt(fVx * fVx + fVy * fVy + fVz * fVz);
  } // Radius of production vertex in spherical system
  Double_t Px() const { return fPx; }
  Double_t Py() const { return fPy; }
  Double_t Pz() const { return fPz; }
  Double_t P() const { return TMath::Sqrt(fPx * fPx + fPy * fPy + fPz * fPz); }
  Double_t Pt() const { return TMath::Sqrt(fPx * fPx + fPy * fPy); }
  Double_t Energy() const { return fE; }
  Double_t Eta() const
  {
    Double_t pmom = P();
    if (pmom != TMath::Abs(fPz))
      return 0.5 * TMath::Log((pmom + fPz) / (pmom - fPz));
    else
      return 1.e30;
  }
  Double_t Y() const
  {
    if (fE != TMath::Abs(fPz))
      return 0.5 * TMath::Log((fE + fPz) / (fE - fPz));
    else
      return 1.e30;
  }
 
  Double_t Phi() const
  {
    return TMath::Pi() + TMath::ATan2(-fPy, -fPx);
  } // note that Phi() returns an angle between 0 and 2pi
  Double_t Theta() const
  {
    return (fPz == 0) ? TMath::PiOver2() : TMath::ACos(fPz / P());
  }
 
  // setters
 
  void SetFirstMother(int code) { fMother[0] = code; }
  void SetMother(int i, int code) { fMother[i] = code; }
  void SetLastMother(int code) { fMother[1] = code; }
  void SetFirstDaughter(int code) { fDaughter[0] = code; }
  void SetDaughter(int i, int code) { fDaughter[i] = code; }
  void SetLastDaughter(int code) { fDaughter[1] = code; }
  void SetCalcMass(Double_t mass) { fCalcMass = mass; }
  void SetPdgCode(Int_t pdg);
  void SetPolarisation(Double_t polx, Double_t poly, Double_t polz);
  void SetPolarisation(const TVector3& v)
  {
    SetPolarisation(v.X(), v.Y(), v.Z());
  }
  void SetStatusCode(int status) { fStatusCode = status; }
  void SetWeight(Float_t weight = 1) { fWeight = weight; }
  void SetMomentum(Double_t px, Double_t py, Double_t pz, Double_t e)
  {
    fPx = px;
    fPy = py;
    fPz = pz;
    fE = e;
  }
  void SetMomentum(const TLorentzVector& p)
  {
    SetMomentum(p.Px(), p.Py(), p.Pz(), p.Energy());
  }
  void SetProductionVertex(Double_t vx, Double_t vy, Double_t vz, Double_t t)
  {
    fVx = vx;
    fVy = vy;
    fVz = vz;
    fVt = t;
  }
  void SetProductionVertex(const TLorentzVector& v)
  {
    SetProductionVertex(v.X(), v.Y(), v.Z(), v.T());
  }
 
 protected:
  // data members
 
  Int_t fPdgCode;     // PDG code of the particle
  Int_t fStatusCode;  // generation status code
  Int_t fMother[2];   // Indices of the mother particles
  Int_t fDaughter[2]; // Indices of the daughter particles
  Float_t fWeight;    // particle weight
 
  Double_t fCalcMass; // Calculated mass
 
  Double_t fPx; // x component of momentum
  Double_t fPy; // y component of momentum
  Double_t fPz; // z component of momentum
  Double_t fE;  // Energy
 
  Double_t fVx; // x of production vertex
  Double_t fVy; // y of production vertex
  Double_t fVz; // z of production vertex
  Double_t fVt; // t of production vertex
 
  Double_t fPolarTheta; // Polar angle of polarisation
  Double_t fPolarPhi;   // azymutal angle of polarisation
 
  Int_t fParticlePDG; // PDG number
  //----------------------------------------------------------------------------
  //  functions
  //----------------------------------------------------------------------------
 
  //   ClassDef(Ex03dParticle,1)  // Ex03dParticle vertex particle information
};
 
#endif