PWG/CaloTrackCorrBase/AliCaloPID.h

   1 #ifndef ALICALOPID_H
   2 #define ALICALOPID_H
   3 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   4  * See cxx source for full Copyright notice     */
   5
   6 //_________________________________________________________________________
   7 // Class for PID selection with calorimeters
   8 // The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster,
   9 // being kPhoton, kElectron, kPi0 ... as defined in the header file
  10 //   - GetIdentifiedParticleType(const AliVCluster * cluster)
  11 //      Assignes a PID tag to the cluster, right now there is the possibility to : use bayesian weights from reco,
  12 //      recalculate them (EMCAL) or use other procedures not used in reco.
  13 //      In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
  14 //      and do SwitchOnBayesianRecalculation().
  15 //      To change the PID parameters from Low to High like the ones by default, use the constructor
  16 //      AliCaloPID(flux)
  17 //      where flux is AliCaloPID::kLow or AliCaloPID::kHigh
  18 //      If it is necessary to change the parameters use the constructor
  19 //      AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
  20
  21 //   - GetGetIdentifiedParticleTypeFromBayesian(const Double_t * pid, const Float_t energy)
  22 //      Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle,
  23 //      executed when bayesian is ON by GetIdentifiedParticleType(const AliVCluster * cluster)
  24 //   - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the
  25 //     result of the PID bayesian a different PID bit is set.
  26 //
  27 //
  28 //*-- Author: Gustavo Conesa (INFN-LNF)
  29
  30 // --- ROOT system ---
  31 #include <TObject.h>
  32 class TString ;
  33 class TLorentzVector ;
  34 #include <TFormula.h>
  35 class TList;
  36 class TH2F ;
  37
  38 //--- AliRoot system ---
  39 class AliVCluster;
  40 class AliVCaloCells;
  41 class AliAODPWG4Particle;
  42 class AliEMCALPIDUtils;
  43 class AliCalorimeterUtils;
  44 class AliVEvent;
  45
  46 class AliCaloPID : public TObject {
  47
  48  public:
  49
  50   AliCaloPID() ; // ctor
  51   AliCaloPID(Int_t particleFlux) ; // ctor, to be used when recalculating bayesian PID
  52   AliCaloPID(const TNamed * emcalpid) ; // ctor, to be used when recalculating bayesian PID and need different parameters
  53   virtual ~AliCaloPID() ;//virtual dtor
  54
  55   enum PidType
  56   {
  57     kPhoton         = 22,
  58     kPi0            = 111,
  59     kEta            = 221,
  60     kElectron       = 11,
  61     kEleCon         =-11,
  62     kNeutralHadron  = 2112,
  63     kChargedHadron  = 211,
  64     kNeutralUnknown = 130,
  65     kChargedUnknown = 321
  66   };
  67
  68   enum TagType {kPi0Decay, kEtaDecay, kOtherDecay, kConversion, kNoTag = -1};
  69
  70   // Main methods
  71
  72   //  TList *   GetCreateOutputObjects(); // Not implemented
  73
  74   void      InitParameters();
  75
  76   Bool_t    IsInPi0SplitAsymmetryRange(Float_t energy, Float_t asy,  Int_t nlm) const;
  77
  78   Bool_t    IsInPi0SplitMassRange     (Float_t energy, Float_t mass, Int_t nlm) const;
  79
  80   Bool_t    IsInM02Range              (Float_t m02) const;
  81   Bool_t    IsInPi0M02Range           (Float_t energy, Float_t m02,  Int_t nlm) const;
  82   Bool_t    IsInEtaM02Range           (Float_t energy, Float_t m02,  Int_t nlm) const;
  83   Bool_t    IsInConM02Range           (Float_t energy, Float_t m02,  Int_t nlm) const;
  84
  85
  86   Int_t     GetIdentifiedParticleTypeFromBayesWeights(Bool_t isEMCAL, Double_t * pid, Float_t energy) ;
  87
  88   Int_t     GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster * cluster, AliVCaloCells* cells,
  89                                                           AliCalorimeterUtils * caloutils,
  90                                                           Double_t vertex[3],
  91                                                           Int_t & nLocMax, Double_t & mass, Double_t & angle,
  92                                                           TLorentzVector & l1  , TLorentzVector & l2,
  93                                                           Int_t   & absId1,   Int_t   & absId2,
  94                                                           Float_t & distbad1, Float_t & distbad2,
  95                                                           Bool_t  & fidcut1,  Bool_t  & fidcut2  ) const;
  96
  97   Int_t     GetIdentifiedParticleType(AliVCluster * cluster) ;
  98
  99   TString   GetPIDParametersList();
 100
 101   Bool_t    IsTrackMatched(AliVCluster * cluster, AliCalorimeterUtils* cu, AliVEvent* event) const ;
 102
 103   void      SetPIDBits(AliVCluster * cluster, AliAODPWG4Particle *aodph,
 104                        AliCalorimeterUtils* cu, AliVEvent* event);
 105
 106   void      Print(const Option_t * opt)const;
 107
 108   void      PrintClusterPIDWeights(const Double_t * pid) const;
 109
 110   //Check if cluster photon-like. Uses photon cluster parameterization in real pp data
 111   //Returns distance in sigmas. Recommended cut 2.5
 112   Float_t TestPHOSDispersion(Double_t pt, Double_t m20, Double_t m02) const ;
 113   //Checks distance to the closest track. Takes into account
 114   //non-perpendicular incidence of tracks.
 115   Float_t TestPHOSChargedVeto(Double_t dx, Double_t dz, Double_t ptTrack,
 116                               Int_t chargeTrack, Double_t mf) const ;
 117
 118   // Setters, getters
 119
 120   void    SetDebug(Int_t deb)                  { fDebug = deb                 ; }
 121   Int_t   GetDebug()                     const { return fDebug                ; }
 122
 123   enum    eventType{kLow,kHigh};
 124   void    SetLowParticleFlux()                 { fParticleFlux        = kLow  ; }
 125   void    SetHighParticleFlux()                { fParticleFlux        = kHigh ; }
 126   // not really used, only for bayesian recalculation in EMCAL, but could be useful in future
 127
 128   // Bayesian
 129
 130   void    SwitchOnBayesian()                   { fUseBayesianWeights  = kTRUE ; }
 131   void    SwitchOffBayesian()                  { fUseBayesianWeights  = kFALSE; }
 132   void    SwitchOnBayesianRecalculation()      { fRecalculateBayesian = kTRUE ; fUseBayesianWeights  = kTRUE ;} // EMCAL
 133   void    SwitchOffBayesianRecalculation()     { fRecalculateBayesian = kFALSE; } // EMCAL
 134
 135   AliEMCALPIDUtils * GetEMCALPIDUtils() ;
 136
 137   //Weight getters
 138   Float_t GetEMCALPhotonWeight()         const { return fEMCALPhotonWeight    ; }
 139   Float_t GetEMCALPi0Weight()            const { return fEMCALPi0Weight       ; }
 140   Float_t GetEMCALElectronWeight()       const { return fEMCALElectronWeight  ; }
 141   Float_t GetEMCALChargeWeight()         const { return fEMCALChargeWeight    ; }
 142   Float_t GetEMCALNeutralWeight()        const { return fEMCALNeutralWeight   ; }
 143   Float_t GetPHOSPhotonWeight()          const { return fPHOSPhotonWeight     ; }
 144   Float_t GetPHOSPi0Weight()             const { return fPHOSPi0Weight        ; }
 145   Float_t GetPHOSElectronWeight()        const { return fPHOSElectronWeight   ; }
 146   Float_t GetPHOSChargeWeight()          const { return fPHOSChargeWeight     ; }
 147   Float_t GetPHOSNeutralWeight()         const { return fPHOSNeutralWeight    ; }
 148
 149   Bool_t  IsPHOSPIDWeightFormulaOn()     const { return fPHOSWeightFormula    ; }
 150
 151   TFormula * GetPHOSPhotonWeightFormula()      {
 152     if(!fPHOSPhotonWeightFormula)
 153       fPHOSPhotonWeightFormula = new TFormula("phos_photon_weight",
 154                                               fPHOSPhotonWeightFormulaExpression);
 155     return fPHOSPhotonWeightFormula                                           ; }
 156
 157   TFormula * GetPHOSPi0WeightFormula()         {
 158     if(!fPHOSPi0WeightFormula)
 159       fPHOSPi0WeightFormula = new TFormula("phos_pi0_weight",
 160                                            fPHOSPi0WeightFormulaExpression);
 161     return fPHOSPi0WeightFormula                                              ; }
 162
 163   TString  GetPHOSPhotonWeightFormulaExpression() const { return fPHOSPhotonWeightFormulaExpression ; }
 164   TString  GetPHOSPi0WeightFormulaExpression()    const { return fPHOSPi0WeightFormulaExpression    ; }
 165
 166   //Weight setters
 167   void    SetEMCALPhotonWeight  (Float_t  w)   { fEMCALPhotonWeight   = w     ; }
 168   void    SetEMCALPi0Weight     (Float_t  w)   { fEMCALPi0Weight      = w     ; }
 169   void    SetEMCALElectronWeight(Float_t  w)   { fEMCALElectronWeight = w     ; }
 170   void    SetEMCALChargeWeight  (Float_t  w)   { fEMCALChargeWeight   = w     ; }
 171   void    SetEMCALNeutralWeight (Float_t  w)   { fEMCALNeutralWeight  = w     ; }
 172   void    SetPHOSPhotonWeight   (Float_t  w)   { fPHOSPhotonWeight    = w     ; }
 173   void    SetPHOSPi0Weight      (Float_t  w)   { fPHOSPi0Weight       = w     ; }
 174   void    SetPHOSElectronWeight (Float_t  w)   { fPHOSElectronWeight  = w     ; }
 175   void    SetPHOSChargeWeight   (Float_t  w)   { fPHOSChargeWeight    = w     ; }
 176   void    SetPHOSNeutralWeight  (Float_t  w)   { fPHOSNeutralWeight   = w     ; }
 177
 178   void    UsePHOSPIDWeightFormula   (Bool_t ok )           { fPHOSWeightFormula                 = ok ; }
 179   void    SetPHOSPhotonWeightFormulaExpression(TString ph) { fPHOSPhotonWeightFormulaExpression = ph ; }
 180   void    SetPHOSPi0WeightFormulaExpression   (TString pi) { fPHOSPi0WeightFormulaExpression    = pi ; }
 181
 182   //PID cuts
 183
 184   void    SetEMCALLambda0CutMax(Float_t lcut ) { fEMCALL0CutMax     = lcut    ; }
 185   Float_t GetEMCALLambda0CutMax()        const { return fEMCALL0CutMax        ; }
 186
 187   void    SetEMCALLambda0CutMin(Float_t lcut ) { fEMCALL0CutMin     = lcut    ; }
 188   Float_t GetEMCALLambda0CutMin()        const { return fEMCALL0CutMin        ; }
 189
 190   void    SetEMCALDEtaCut(Float_t dcut )       { fEMCALDEtaCut      = dcut    ; }
 191   Float_t GetEMCALDEtaCut()              const { return fEMCALDEtaCut         ; }
 192
 193   void    SetEMCALDPhiCut(Float_t dcut )       { fEMCALDPhiCut      = dcut    ; }
 194   Float_t GetEMCALDPhiCut()              const { return fEMCALDPhiCut         ; }
 195
 196   void    SetTOFCut(Float_t tcut )             { fTOFCut            = tcut    ; }
 197   Float_t GetTOFCut()                    const { return fTOFCut               ; }
 198
 199   void    SetPHOSRCut(Float_t rcut )           { fPHOSRCut          = rcut    ; }
 200   Float_t GetPHOSRCut()                  const { return fPHOSRCut             ; }
 201
 202   void    SetPHOSDispersionCut(Float_t dcut )  { fPHOSDispersionCut = dcut    ; }
 203   Float_t GetPHOSDispersionCut()         const { return fPHOSDispersionCut    ; }
 204
 205   // Cluster splitting analysis
 206
 207   void    SwitchOnSimpleSplitMassCut()         { fUseSimpleMassCut   = kTRUE  ; }
 208   void    SwitchOffSimpleSplitMassCut()        { fUseSimpleMassCut   = kFALSE ; }
 209
 210   void    SwitchOnSimpleSplitM02Cut()          { fUseSimpleM02Cut    = kTRUE  ; }
 211   void    SwitchOffSimpleSplitM02Cut()         { fUseSimpleM02Cut    = kFALSE ; }
 212
 213   void    SwitchOnSplitAsymmetryCut()          { fUseSplitAsyCut     = kTRUE  ; }
 214   void    SwitchOffSplitAsymmetryCut()         { fUseSplitAsyCut     = kFALSE ; }
 215   Bool_t  IsSplitAsymmetryCutOn()              { return fUseSplitAsyCut       ; }
 216
 217   void    SwitchOnSplitShowerShapeCut()        { fUseSplitSSCut      = kTRUE  ; }
 218   void    SwitchOffSplitShowerShapeCut()       { fUseSplitSSCut      = kFALSE ; }
 219   Bool_t  IsSplitShowerShapeCutOn()            { return fUseSplitSSCut        ; }
 220
 221   void    SetClusterSplittingM02Cut(Float_t min=0, Float_t max=100)
 222   { fSplitM02MinCut   = min ; fSplitM02MaxCut  = max ; }
 223
 224   void    SetClusterSplittingMinNCells(Int_t c) { fSplitMinNCells = c         ; }
 225   Int_t   GetClusterSplittingMinNCells() const  { return fSplitMinNCells      ; }
 226
 227   void    SetSplitEnergyFractionMinimum(Int_t i, Float_t min){ if (i < 3 && i >=0 ) fSplitEFracMin[i]  = min ; }
 228   Float_t GetSplitEnergyFractionMinimum(Int_t i) const       { if( i < 3 && i >=0 ) return fSplitEFracMin[i] ;  else return 0 ; }
 229
 230   void    SetSubClusterEnergyMinimum   (Int_t i, Float_t min){ if (i < 3 && i >=0 ) fSubClusterEMin[i] = min ; }
 231   Float_t GetSubClusterEnergyMinimum   (Int_t i) const       { if( i < 3 && i >=0 ) return fSubClusterEMin[i];  else return 0 ; }
 232
 233   Float_t GetPi0MinMass()                const { return fMassPi0Min           ; } // Simple cut case
 234   Float_t GetEtaMinMass()                const { return fMassEtaMin           ; } // Simple cut case
 235   Float_t GetPhotonMinMass()             const { return fMassPhoMin           ; }
 236   Float_t GetPi0MaxMass()                const { return fMassPi0Max           ; }
 237   Float_t GetEtaMaxMass()                const { return fMassEtaMax           ; }
 238   Float_t GetPhotonMaxMass()             const { return fMassPhoMax           ; }
 239
 240   void    SetSplitWidthSigma(Float_t s)        { fSplitWidthSigma        = s  ; }
 241
 242   void    SetPi0MassShiftHighECell(Float_t s)  { fMassShiftHighECell     = s  ; }
 243
 244   void    SetPi0MassSelectionParameters    (Int_t inlm, Int_t iparam, Float_t param)
 245   { if(iparam < 6 ) fMassPi0Param[inlm][iparam] = param ; }
 246
 247   void    SetPi0WidthSelectionParameters    (Int_t inlm, Int_t iparam, Float_t param)
 248   { if(iparam < 6 ) fWidthPi0Param[inlm][iparam] = param ; }
 249
 250   void    SetM02MaximumSelectionParameters      (Int_t inlm, Int_t iparam, Float_t param)
 251   { if(iparam < 5 && inlm < 2) fM02MaxParam[inlm][iparam] = param ; }
 252
 253   void    SetM02MaximumShiftForNLMN(Int_t shift) { fM02MaxParamShiftNLMN = shift ; }
 254
 255   void    SetM02MinimumSelectionParameters      (Int_t inlm, Int_t iparam, Float_t param)
 256   { if(iparam < 5 && inlm < 2) fM02MinParam[inlm][iparam] = param ; }
 257
 258   void    SetAsymmetryMinimumSelectionParameters(Int_t inlm, Int_t iparam, Float_t param)
 259   { if(iparam < 4 && inlm < 2) fAsyMinParam[inlm][iparam] = param ; }
 260
 261   void    SetPi0MassRange(Float_t min, Float_t max)    { fMassPi0Min    = min ; fMassPi0Max = max ; } // Simple case
 262   void    SetEtaMassRange(Float_t min, Float_t max)    { fMassEtaMin    = min ; fMassEtaMax = max ; }
 263   void    SetPhotonMassRange(Float_t min, Float_t max) { fMassPhoMin    = min ; fMassPhoMax = max ; }
 264
 265 private:
 266
 267   Int_t     fDebug;                             // Debug level
 268   Int_t     fParticleFlux;                      // Particle flux for setting PID parameters
 269
 270   // Bayesian
 271   AliEMCALPIDUtils * fEMCALPIDUtils;            // Pointer to EMCALPID to redo the PID Bayesian calculation
 272   Bool_t    fUseBayesianWeights;                // Select clusters based on weights calculated in reconstruction
 273   Bool_t    fRecalculateBayesian;               // Recalculate PID bayesian or use simple PID?
 274
 275   Float_t   fEMCALPhotonWeight;                 // Bayesian PID weight for photons in EMCAL
 276   Float_t   fEMCALPi0Weight;                    // Bayesian PID weight for pi0 in EMCAL
 277   Float_t   fEMCALElectronWeight;               // Bayesian PID weight for electrons in EMCAL
 278   Float_t   fEMCALChargeWeight;                 // Bayesian PID weight for charged hadrons in EMCAL
 279   Float_t   fEMCALNeutralWeight;                // Bayesian PID weight for neutral hadrons in EMCAL
 280   Float_t   fPHOSPhotonWeight;                  // Bayesian PID weight for photons in PHOS
 281   Float_t   fPHOSPi0Weight;                     // Bayesian PID weight for pi0 in PHOS
 282   Float_t   fPHOSElectronWeight;                // Bayesian PID weight for electrons in PHOS
 283   Float_t   fPHOSChargeWeight;                  // Bayesian PID weight for charged hadrons in PHOS
 284   Float_t   fPHOSNeutralWeight;                 // Bayesian PID weight for neutral hadrons in PHOS
 285
 286   Bool_t    fPHOSWeightFormula ;                // Use parametrized weight threshold, function of energy
 287   TFormula *fPHOSPhotonWeightFormula ;          // Formula for photon weight
 288   TFormula *fPHOSPi0WeightFormula ;             // Formula for pi0 weight
 289   TString   fPHOSPhotonWeightFormulaExpression; // Photon weight formula in string
 290   TString   fPHOSPi0WeightFormulaExpression;    // Pi0 weight formula in string
 291
 292   // PID calculation
 293   Float_t   fEMCALL0CutMax;                     // Max Cut on shower shape lambda0, used in PID evaluation, only EMCAL
 294   Float_t   fEMCALL0CutMin;                     // Min Cut on shower shape lambda0, used in PID evaluation, only EMCAL
 295   Float_t   fEMCALDEtaCut;                      // Track matching cut on Dz
 296   Float_t   fEMCALDPhiCut;                      // Track matching cut on Dx
 297
 298   Float_t   fTOFCut;                            // Cut on TOF, used in PID evaluation
 299
 300   Float_t   fPHOSDispersionCut;                 // Shower shape elipse radious cut
 301   Float_t   fPHOSRCut;                          // Track-Cluster distance cut for track matching in PHOS
 302
 303   // Cluster splitting mass ranges
 304   Bool_t    fUseSimpleMassCut;                  // Use simple min-max pi0 mass cut
 305   Bool_t    fUseSimpleM02Cut;                   // Use simple min-max M02 cut
 306   Bool_t    fUseSplitAsyCut ;                   // Remove splitted clusters with too large asymmetry
 307   Bool_t    fUseSplitSSCut  ;                   // Remove splitted clusters out of shower shape band
 308   Float_t   fSplitM02MaxCut ;                   // Study clusters with l0 smaller than cut
 309   Float_t   fSplitM02MinCut ;                   // Study clusters with l0 larger than cut  // simple case
 310   Int_t     fSplitMinNCells ;                   // Study clusters with ncells larger than cut
 311   Float_t   fMassEtaMin  ;                      // Min Eta mass
 312   Float_t   fMassEtaMax  ;                      // Max Eta mass
 313   Float_t   fMassPi0Min  ;                      // Min Pi0 mass // simple cut case
 314   Float_t   fMassPi0Max  ;                      // Min Pi0 mass // simple cut case
 315   Float_t   fMassPhoMin  ;                      // Min Photon mass
 316   Float_t   fMassPhoMax  ;                      // Min Photon mass
 317   Float_t   fMassPi0Param [2][6] ;              // mean mass param, 2 regions in energy
 318   Float_t   fWidthPi0Param[2][6] ;              // width param, 2 regions in energy
 319   Float_t   fM02MinParam[2][5] ;                // 5 param for expo + pol fit on M02 minimum for pi0 selection (maximum for conversions)
 320   Float_t   fM02MaxParam[2][5] ;                // 5 param for expo + pol fit on M02 maximum for pi0 selection
 321   Float_t   fM02MaxParamShiftNLMN;              // shift of max M02 for NLM>2
 322   Float_t   fAsyMinParam[2][4] ;                // 3 param for fit on asymmetry minimum, for 2 cases, NLM=1 and NLM>=2
 323   Float_t   fSplitEFracMin[3]  ;                // Do not use clusters with too large energy in cluster compared
 324                                                 // to energy in splitted clusters, depeding on NLM
 325   Float_t   fSubClusterEMin[3]  ;               // Do not use sub-clusters with too low energy depeding on NLM
 326   Float_t   fSplitWidthSigma;                   // Cut on mass+-width*fSplitWidthSigma
 327   Float_t   fMassShiftHighECell;                // Shift cuts 5 MeV for Ecell > 150 MeV, default Ecell > 50 MeV
 328
 329   AliCaloPID & operator = (const AliCaloPID & cpid) ; // cpy assignment
 330   AliCaloPID(              const AliCaloPID & cpid) ; // cpy ctor
 331
 332   ClassDef(AliCaloPID,22)
 333
 334 } ;
 335
 336
 337 #endif //ALICALOPID_H
 338
 339
 340