MpdPid.h 9.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189
  1. //------------------------------------------------------------------------------------------------------------------------
  2. #ifndef MPD_PID_H
  3. #define MPD_PID_H
  4. //------------------------------------------------------------------------------------------------------------------------
  5. /// \class MpdPid
  6. ///
  7. /// \brief
  8. /// \author Alexander Mudrokh (LHEP, JINR, Dubna)
  9. //------------------------------------------------------------------------------------------------------------------------
  10. #include "MpdPidUtils.h"
  11. #include "MpdTrack.h"
  12. #include "MpdTpcKalmanTrack.h" /// MiniDST
  13. #include "FairRunAna.h" /// MiniDST
  14. #include <TFile.h>
  15. #include <TF2.h>
  16. #include <TMatrixD.h> /// MiniDST
  17. //------------------------------------------------------------------------------------------------------------------------
  18. using namespace std;
  19. typedef vector<TF1*> vecTF1ptrs;
  20. class MpdPid : public TObject
  21. {
  22. public:
  23. //--------------------------------------------------------------------//
  24. //--------------------- Constructor / Destructor ---------------------//
  25. //--------------------------------------------------------------------//
  26. MpdPid(); /// default ctor
  27. MpdPid(Double_t sigmaTof, Double_t sigmaEloss, Double_t sqrts, /// generators: "PHSD", "URQMD", "LAQGSM" ("QGSM"), "DEFAULT", "EPOS" (for pp collisions), "NSIG" (for n-sigma method)
  28. Double_t EnLossCoef = 1., TString Generator = "DEFAULT", TString Tracking = "CFHM", /// tracking: "HP" (Hit Producer), "CF" (Cluster Finder)
  29. TString NSigPart = "pikapr"); /// possible expressions: pi, ka, pr, el, mu, de, tr, he3, he4
  30. virtual ~MpdPid(); /// destructor
  31. //--------------------------------------------------------------------//
  32. //---------------------- Array of probabilities ----------------------//
  33. //--------------------------------------------------------------------//
  34. /// Fill array of probabilities, otherwise return kFALSE
  35. Bool_t FillProbs(MpdTrack*);
  36. Bool_t FillProbs(Double_t, Double_t, Int_t); /// variables: full momentum, dE/dx, charge
  37. Bool_t FillProbs(Double_t, Double_t, Double_t, Int_t); /// variables: full momentum, dE/dx, mass squared, charge
  38. /// Return probabilities
  39. Double_t GetProbEl(void) {return fProb[MpdPidUtils::kElectron];}
  40. Double_t GetProbMu(void) {return fProb[MpdPidUtils::kMuon];}
  41. Double_t GetProbPi(void) {return fProb[MpdPidUtils::kPion];}
  42. Double_t GetProbKa(void) {return fProb[MpdPidUtils::kKaon];}
  43. Double_t GetProbPr(void) {return fProb[MpdPidUtils::kProton];}
  44. Double_t GetProbDe(void) {return fProb[MpdPidUtils::kDeuteron];}
  45. Double_t GetProbTr(void) {return fProb[MpdPidUtils::kTriton];}
  46. Double_t GetProbHe3(void) {return fProb[MpdPidUtils::kHe3];}
  47. Double_t GetProbHe4(void) {return fProb[MpdPidUtils::kHe4];}
  48. Double_t GetProb(MpdPidUtils::ePartType iType) {return fProb[iType];}
  49. Long_t GetMaxProb(); /// Returns the most probable pdg code
  50. //--------------------------------------------------------------------//
  51. //------------------------------ dE/dx -------------------------------//
  52. //--------------------------------------------------------------------//
  53. /// Returns expected <dE/dx> value (variable: full momentum)
  54. Double_t GetDedxElParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kElectron); }
  55. Double_t GetDedxMuParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kMuon); }
  56. Double_t GetDedxPiParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kPion); }
  57. Double_t GetDedxKaParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kKaon); }
  58. Double_t GetDedxPrParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kProton); }
  59. Double_t GetDedxDeParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kDeuteron); }
  60. Double_t GetDedxTrParam(Double_t p) { return GetDedxParam(p, MpdPidUtils::kTriton); }
  61. Double_t GetDedxHe3Param(Double_t p) { return GetDedxParam(p, MpdPidUtils::kHe3); }
  62. Double_t GetDedxHe4Param(Double_t p) { return GetDedxParam(p, MpdPidUtils::kHe4); }
  63. Double_t GetDedxParam(Double_t, MpdPidUtils::ePartType);
  64. /// Returns expected dE/dx width (variables: full momentum, species)
  65. Double_t GetDedxWidthValue(Double_t, MpdPidUtils::ePartType);
  66. /// Returns expected asymmetry value delta (variables: full momentum, species)
  67. Double_t GetTailValue(Double_t, MpdPidUtils::ePartType);
  68. /// Returns expexted <dE/dx>, dE/dx width and asymmetry parameterizations as vectors of TF1*
  69. vecTF1ptrs GetVecdEdxMean(MpdPidUtils::ePartType);
  70. vecTF1ptrs GetVecdEdxWidth(MpdPidUtils::ePartType);
  71. vecTF1ptrs GetVecdEdxAsym(MpdPidUtils::ePartType);
  72. //--------------------------------------------------------------------//
  73. //------------------------------- m^2 --------------------------------//
  74. //--------------------------------------------------------------------//
  75. /// Returns expected m^2 width (variables: full momentum, species)
  76. Double_t Getm2WidthParam(Double_t, MpdPidUtils::ePartType);
  77. /// Returns expected m^2 width parameterization as vector of TF1*
  78. vecTF1ptrs GetVecm2Width(MpdPidUtils::ePartType);
  79. //--------------------------------------------------------------------//
  80. //----------------------------- Yileds -------------------------------//
  81. //--------------------------------------------------------------------//
  82. Double_t GetBayesCoefficient(Double_t, MpdPidUtils::ePartType, MpdPidUtils::ePartCharge);
  83. Double_t GetPrRat() {return fPrRatio;} /// Returns pos./neg. ratio for protons
  84. void SetPrRat(Double_t); /// Set pos./neg. ratio for protons
  85. /// Returns expected yield parameterization as vector of TF1*
  86. /// use via vec[MpdPidUtils::kPos/kNeg]
  87. vecTF1ptrs GetVecYield(MpdPidUtils::ePartType);
  88. //--------------------------------------------------------------------//
  89. //------------------------------ Other -------------------------------//
  90. //--------------------------------------------------------------------//
  91. void Print(const char* comment, ostream& os);
  92. /// Returns distance to the most probable dE/dx and m2 value in terms of sigmas
  93. /// possible expressions: pi, ka, pr, el, mu, de, tr, he3, he4
  94. /// use it after MpdPid::FillProbs usage only!!!
  95. Double_t GetNsigmaToBetheBloch(TString);
  96. Double_t GetNsigmaToBetheBloch(MpdPidUtils::ePartType);
  97. Double_t GetNsigmaToAverageMass2(TString);
  98. Double_t GetNsigmaToAverageMass2(MpdPidUtils::ePartType);
  99. /// Includes particle species iType to n-sigma method
  100. void SetParticleEnabledNsig(MpdPidUtils::ePartType iType) { fNSigSpecies[iType] = kTRUE; }
  101. /// Additional functions
  102. Double_t parElBB(Double_t*, Double_t*);
  103. Double_t parMuBB(Double_t*, Double_t*);
  104. Double_t parPiBB(Double_t*, Double_t*);
  105. Double_t parKaBB(Double_t*, Double_t*);
  106. Double_t parPrBB(Double_t*, Double_t*);
  107. Double_t parDeBB(Double_t*, Double_t*);
  108. Double_t parTrBB(Double_t*, Double_t*);
  109. Double_t parHe3BB(Double_t*, Double_t*);
  110. Double_t parHe4BB(Double_t*, Double_t*);
  111. Double_t AsymGaus(Double_t*, Double_t*);
  112. Double_t AsymGaus2(Double_t*, Double_t*);
  113. Double_t MomPi(Double_t* , Double_t*);
  114. Double_t MomPr(Double_t*, Double_t*);
  115. protected:
  116. void Init(TString, TString, TString, Double_t);
  117. /// GetDedxProb_asym variables: cut_dedx, p, dedx, n, emean, sige, species
  118. Double_t ComputeDedxProb_asym(Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, MpdPidUtils::ePartType);
  119. /// GetCombProb_asym variables: cut_dedx, cut_m2, p, dedx, m2, n, emean, mmean, sige, sigm, species
  120. Double_t ComputeCombProb_asym(Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, Double_t, MpdPidUtils::ePartType);
  121. void ComputeBayesCoefficients(Double_t);
  122. Double_t ComputeEnLossSigma(Double_t, Double_t, Double_t, Double_t, MpdPidUtils::ePartType);
  123. Double_t ComputeMSquaredSigma(Double_t, Double_t, Double_t, MpdPidUtils::ePartType);
  124. /// Variables
  125. TF1 *fGaus;
  126. TF1 *fAsymGaus;
  127. TF2 *fGaus2;
  128. TF2 *fAsymGaus2;
  129. Double_t fProb[MpdPidUtils::kNSpecies+1]; ///< the probability to identify track as a species <i>
  130. Double_t fEnLossSigmasArray[MpdPidUtils::kNSpecies]; ///< the deviation of the measured energy loss from that expected for the species <i>, in terms of the detector resolution
  131. Double_t fMSquaredSigmasArray[MpdPidUtils::kNSpecies]; ///< the deviation of the measured mass squared from that expected for the species <i>, in terms of the detector resolution
  132. Double_t fPrRatio; ///< proton ratio is pos./neg.
  133. Double_t fEnergy; ///< collision energy
  134. Double_t fSigmaTof; ///< non-zero distance from the average mass-squared value (in terms of standard deviations)
  135. Double_t fSigmaEloss; ///< non-zero distance from the average dE/dx value (in terms of standard deviations)
  136. MpdPidUtils::eTrackingState fTrackingState;
  137. Bool_t fMethod; ///< PID method, kTRUE - Bayesian approach, kFALSE - n-sigma method
  138. MpdPidUtils::ePartCharge fCharge; ///< track charge
  139. Double_t fBayesCoefficients[MpdPidUtils::kNSpecies]; ///< Bayes Coefficients assigned to the track
  140. Bool_t fNSigSpecies[MpdPidUtils::kNSpecies]; ///< array of flags to the particle species involved in n-sigma method
  141. map <MpdPidUtils::ePartType,vecTF1ptrs> fdEdxBBMap; /// Bethe-Bloch functions for mean energy deposit description
  142. map <MpdPidUtils::ePartType,vecTF1ptrs> fdEdxSigmaMap; /// Polynomials for dE/dx width description
  143. map <MpdPidUtils::ePartType,vecTF1ptrs> fdEdxDeltaMap; /// Polynomials for dE/dx asymmetry description
  144. map <MpdPidUtils::ePartType,vecTF1ptrs> fParM2Map; /// Functions for m2 width description
  145. map <MpdPidUtils::ePartType,vecTF1ptrs> fPartYieldMap; /// Functions for multiplicity description
  146. ClassDef(MpdPid,4);
  147. };
  148. #endif