KFParticleBase.h

Go to the documentation of this file.

//---------------------------------------------------------------------------------
// The KFParticleBase class
// .
// @author  S.Gorbunov, I.Kisel, I.Kulakov, M.Zyzak
// @version 1.0
// @since   13.05.07
//
// Class to reconstruct and store the decayed particle parameters.
// The method is described in CBM-SOFT note 2007-003,
// ``Reconstruction of decayed particles based on the Kalman filter'',
// http://www.gsi.de/documents/DOC-2007-May-14-1.pdf
//
// This class describes general mathematics which is used by KFParticle class
//
//  -= Copyright &copy ALICE HLT and CBM L1 Groups =-
//_________________________________________________________________________________
 
#ifndef ALIKFPARTICLEBASE_H
#define ALIKFPARTICLEBASE_H
 
#include "TObject.h"
 
#include <vector>
 
class KFParticleBase : public TObject
{
 
  public:
    //*
    //* ABSTRACT METHODS HAVE TO BE DEFINED IN USER CLASS
    //*
 
    //* Virtual method to access the magnetic field
 
    virtual void GetFieldValue(const Double_t xyz[], Double_t B[]) const = 0;
 
    //* Virtual methods needed for particle transportation
    //* One can use particular implementations for collider (only Bz component)
    //* geometry and for fixed-target (CBM-like) geometry which are provided below
    //* in TRANSPORT section
 
    //* Get dS to xyz[] space point
 
    virtual Double_t GetDStoPoint(const Double_t xyz[]) const = 0;
 
    //* Get dS to other particle p (dSp for particle p also returned)
 
    virtual void GetDStoParticle(const KFParticleBase& p, Double_t& DS, Double_t& DSp) const = 0;
 
    //* Transport on dS value along trajectory, output to P,C
 
    virtual void Transport(Double_t dS, Double_t P[], Double_t C[]) const = 0;
 
    //*
    //*  INITIALIZATION
    //*
 
    //* Constructor
 
    KFParticleBase();
 
    //* Destructor
 
    virtual ~KFParticleBase() { ; }
 
    //* Initialisation from "cartesian" coordinates ( X Y Z Px Py Pz )
    //* Parameters, covariance matrix, charge, and mass hypothesis should be provided
 
    void Initialize(const Double_t Param[], const Double_t Cov[], Int_t Charge, Double_t Mass);
 
    //* Initialise covariance matrix and set current parameters to 0.0
 
    void Initialize();
 
    //* Set decay vertex parameters for linearisation
 
    void SetVtxGuess(Double_t x, Double_t y, Double_t z);
 
    //* Set consruction method
 
    void SetConstructMethod(Int_t m) { fConstructMethod = m; }
 
    //* Set and get mass hypothesis of the particle
    void SetMassHypo(Double_t m) { fMassHypo = m; }
    const Double_t& GetMassHypo() const { return fMassHypo; }
 
    //* Returns the sum of masses of the daughters
    const Double_t& GetSumDaughterMass() const { return SumDaughterMass; }
 
    //*
    //*  ACCESSORS
    //*
 
    //* Simple accessors
 
    Double_t GetX() const { return fP[0]; }
    Double_t GetY() const { return fP[1]; }
    Double_t GetZ() const { return fP[2]; }
    Double_t GetPx() const { return fP[3]; }
    Double_t GetPy() const { return fP[4]; }
    Double_t GetPz() const { return fP[5]; }
    Double_t GetE() const { return fP[6]; }
    Double_t GetS() const { return fP[7]; }
    Int_t GetQ() const { return fQ; }
    Double_t GetChi2() const { return fChi2; }
    Int_t GetNDF() const { return fNDF; }
 
    const Double_t& X() const { return fP[0]; }
    const Double_t& Y() const { return fP[1]; }
    const Double_t& Z() const { return fP[2]; }
    const Double_t& Px() const { return fP[3]; }
    const Double_t& Py() const { return fP[4]; }
    const Double_t& Pz() const { return fP[5]; }
    const Double_t& E() const { return fP[6]; }
    const Double_t& S() const { return fP[7]; }
    const Int_t& Q() const { return fQ; }
    const Double_t& Chi2() const { return fChi2; }
    const Int_t& NDF() const { return fNDF; }
 
    Double_t GetParameter(Int_t i) const { return fP[i]; }
    Double_t GetCovariance(Int_t i) const { return fC[i]; }
    Double_t GetCovariance(Int_t i, Int_t j) const { return fC[IJ(i, j)]; }
 
    //* Accessors with calculations( &value, &estimated sigma )
    //* error flag returned (0 means no error during calculations)
 
    Int_t GetMomentum(Double_t& P, Double_t& SigmaP) const;
    Int_t GetPt(Double_t& Pt, Double_t& SigmaPt) const;
    Int_t GetEta(Double_t& Eta, Double_t& SigmaEta) const;
    Int_t GetPhi(Double_t& Phi, Double_t& SigmaPhi) const;
    Int_t GetMass(Double_t& M, Double_t& SigmaM) const;
    Int_t GetDecayLength(Double_t& L, Double_t& SigmaL) const;
    Int_t GetDecayLengthXY(Double_t& L, Double_t& SigmaL) const;
    Int_t GetLifeTime(Double_t& T, Double_t& SigmaT) const;
    Int_t GetR(Double_t& R, Double_t& SigmaR) const;
 
    //*
    //*  MODIFIERS
    //*
 
    Double_t& X() { return fP[0]; }
    Double_t& Y() { return fP[1]; }
    Double_t& Z() { return fP[2]; }
    Double_t& Px() { return fP[3]; }
    Double_t& Py() { return fP[4]; }
    Double_t& Pz() { return fP[5]; }
    Double_t& E() { return fP[6]; }
    Double_t& S() { return fP[7]; }
    Int_t& Q() { return fQ; }
    Double_t& Chi2() { return fChi2; }
    Int_t& NDF() { return fNDF; }
 
    Double_t& Parameter(Int_t i) { return fP[i]; }
    Double_t& Covariance(Int_t i) { return fC[i]; }
    Double_t& Covariance(Int_t i, Int_t j) { return fC[IJ(i, j)]; }
 
    //*
    //* CONSTRUCTION OF THE PARTICLE BY ITS DAUGHTERS AND MOTHER
    //* USING THE KALMAN FILTER METHOD
    //*
 
    //* Simple way to add daughter ex. D0+= Pion;
 
    void operator+=(const KFParticleBase& Daughter);
 
    //* Add daughter track to the particle
 
    void AddDaughter(const KFParticleBase& Daughter);
 
    void AddDaughterWithEnergyFit(const KFParticleBase& Daughter);
    void AddDaughterWithEnergyCalc(const KFParticleBase& Daughter);
    void AddDaughterWithEnergyFitMC(const KFParticleBase& Daughter);   // with mass constrained
 
    //* Set production vertex
 
    void SetProductionVertex(const KFParticleBase& Vtx);
 
    //* Set mass constraint
 
    void SetNonlinearMassConstraint(Double_t Mass);
    void SetMassConstraint(Double_t Mass, Double_t SigmaMass = 0);
 
    //* Set no decay length for resonances
 
    void SetNoDecayLength();
 
    //* Everything in one go
 
    void Construct(const KFParticleBase* vDaughters[],
                   Int_t nDaughters,
                   const KFParticleBase* ProdVtx = 0,
                   Double_t Mass = -1,
                   Bool_t IsConstrained = 0);
 
    //*
    //*                   TRANSPORT
    //*
    //*  ( main transportation parameter is S = SignedPath/Momentum )
    //*  ( parameters of decay & production vertices are stored locally )
    //*
 
    //* Transport the particle to its decay vertex
 
    void TransportToDecayVertex();
 
    //* Transport the particle to its production vertex
 
    void TransportToProductionVertex();
 
    //* Transport the particle on dS parameter (SignedPath/Momentum)
 
    void TransportToDS(Double_t dS);
 
    //* Particular extrapolators one can use
 
    Double_t GetDStoPointBz(Double_t Bz, const Double_t xyz[]) const;
    Double_t GetDStoPointBy(Double_t By, const Double_t xyz[]) const;
 
    void GetDStoParticleBz(Double_t Bz, const KFParticleBase& p, Double_t& dS, Double_t& dS1) const;
    void GetDStoParticleBy(Double_t B, const KFParticleBase& p, Double_t& dS, Double_t& dS1) const;
 
    Double_t GetDStoPointCBM(const Double_t xyz[]) const;
    void GetDStoParticleCBM(const KFParticleBase& p, Double_t& dS, Double_t& dS1) const;
 
    void TransportBz(Double_t Bz, Double_t dS, Double_t P[], Double_t C[]) const;
    void TransportCBM(Double_t dS, Double_t P[], Double_t C[]) const;
 
    //*
    //* OTHER UTILITIES
    //*
 
    //* Calculate distance from another object [cm]
 
    Double_t GetDistanceFromVertex(const Double_t vtx[]) const;
    Double_t GetDistanceFromVertex(const KFParticleBase& Vtx) const;
    Double_t GetDistanceFromParticle(const KFParticleBase& p) const;
 
    //* Calculate sqrt(Chi2/ndf) deviation from vertex
    //* v = [xyz], Cv=[Cxx,Cxy,Cyy,Cxz,Cyz,Czz]-covariance matrix
 
    Double_t GetDeviationFromVertex(const Double_t v[], const Double_t Cv[] = 0) const;
    Double_t GetDeviationFromVertex(const KFParticleBase& Vtx) const;
    Double_t GetDeviationFromParticle(const KFParticleBase& p) const;
 
    //* Subtract the particle from the vertex
 
    void SubtractFromVertex(KFParticleBase& Vtx) const;
 
    //* Special method for creating gammas
 
    void ConstructGammaBz(const KFParticleBase& daughter1, const KFParticleBase& daughter2, double Bz);
 
    //* return parameters for the Armenteros-Podolanski plot
    static void GetArmenterosPodolanski(KFParticleBase& positive, KFParticleBase& negative, Double_t QtAlfa[2]);
 
    //* Rotates the KFParticle object around OZ axis, OZ axis is set by the vertex position
    void RotateXY(Double_t angle, Double_t Vtx[3]);
 
    int Id() const { return fId; };
    int NDaughters() const { return fDaughtersIds.size(); };
    const std::vector<int>& DaughterIds() const { return fDaughtersIds; };
 
    void SetId(int id) { fId = id; };   // should be always used (manualy)
    void AddDaughterId(int id) { fDaughtersIds.push_back(id); };
 
    void SetPDG(int pdg) { fPDG = pdg; }
    int GetPDG() const { return fPDG; }
 
  protected:
    static Int_t IJ(Int_t i, Int_t j) { return (j <= i) ? i * (i + 1) / 2 + j : j * (j + 1) / 2 + i; }
 
    Double_t& Cij(Int_t i, Int_t j) { return fC[IJ(i, j)]; }
 
    void Convert(bool ToProduction);
    void TransportLine(Double_t S, Double_t P[], Double_t C[]) const;
    Double_t GetDStoPointLine(const Double_t xyz[]) const;
    void GetDStoParticleLine(const KFParticleBase& p, Double_t& dS, Double_t& dS1) const;
 
    void GetDSIter(const KFParticleBase& p, Double_t const& dS, Double_t x[3], Double_t dx[3], Double_t ddx[3]) const;
 
    static Bool_t InvertSym3(const Double_t A[], Double_t Ainv[]);
 
    static void MultQSQt(const Double_t Q[], const Double_t S[], Double_t SOut[]);
 
    static Double_t GetSCorrection(const Double_t Part[], const Double_t XYZ[]);
 
    void GetMeasurement(const Double_t XYZ[], Double_t m[], Double_t V[]) const;
 
    //* Mass constraint function. is needed for the nonlinear mass constraint and a fit with mass constraint
    void SetMassConstraint(Double_t* mP, Double_t* mC, Double_t mJ[7][7], Double_t mass);
 
    Double_t fP[8];    //* Main particle parameters {X,Y,Z,Px,Py,Pz,E,S[=DecayLength/P]}
    Double_t fC[36];   //* Low-triangle covariance matrix of fP
    Int_t fQ;          //* Particle charge
    Int_t fNDF;        //* Number of degrees of freedom
    Double_t fChi2;    //* Chi^2
 
    Double_t fSFromDecay;   //* Distance from decay vertex to current position
 
    Bool_t fAtProductionVertex;   //* Flag shows that the particle error along
                                  //* its trajectory is taken from production vertex
 
    Double_t fVtxGuess[3];   //* Guess for the position of the decay vertex
                             //* ( used for linearisation of equations )
 
    Bool_t fIsLinearized;   //* Flag shows that the guess is present
 
    Int_t fConstructMethod;   //* Determines the method for the particle construction.
    //* 0 - Energy considered as an independent veriable, fitted independently from momentum, without any constraints on
    // mass
    //* 1 - Energy considered as a dependent variable, calculated from the momentum and mass hypothesis
    //* 2 - Energy considered as an independent variable, fitted independently from momentum, with constraints on mass
    // of daughter particle
 
    Double_t SumDaughterMass;   //* sum of the daughter particles masses
    Double_t fMassHypo;         //* sum of the daughter particles masse
 
    int fId;                          // id of particle
    std::vector<int> fDaughtersIds;   // id of particles it created from. if size == 1 then this is id of track. TODO
                                      // use in functions. why unsigned short int doesn't work???
 
    int fPDG;   // pdg hypothesis
 
    ClassDef(KFParticleBase, 1);
};
 
#endif