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 "RootTypesDef.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