PWG/CaloTrackCorrBase/AliCaloPID.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 // Class for PID selection with calorimeters
  18 // The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster,
  19 // being kPhoton, kElectron, kPi0 ... as defined in the header file
  20 //   - GetIdentifiedParticleType(const AliVCluster * cluster)
  21 //      Assignes a PID tag to the cluster, right now there is the possibility to : use bayesian weights from reco,
  22 //      recalculate them (EMCAL) or use other procedures not used in reco.
  23 //      In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
  24 //      and do SwitchOnBayesianRecalculation().
  25 //      To change the PID parameters from Low to High like the ones by default, use the constructor
  26 //      AliCaloPID(flux)
  27 //      where flux is AliCaloPID::kLow or AliCaloPID::kHigh
  28 //      If it is necessary to change the parameters use the constructor
  29 //      AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
  30
  31 //   - GetGetIdentifiedParticleTypeFromBayesian(const Double_t * pid, const Float_t energy)
  32 //      Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle,
  33 //      executed when bayesian is ON by GetIdentifiedParticleType(const AliVCluster * cluster)
  34 //   - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the
  35 //     result of the PID bayesian a different PID bit is set.
  36 //
  37 //   - IsTrackMatched(): Independent method that needs to be combined with GetIdentifiedParticleType to know if the cluster was matched
  38 //
  39 //*-- Author: Gustavo Conesa (INFN-LNF)
  40 //////////////////////////////////////////////////////////////////////////////
  41
  42
  43 // --- ROOT system ---
  44 #include <TMath.h>
  45 #include <TString.h>
  46 #include <TList.h>
  47
  48 // ---- ANALYSIS system ----
  49 #include "AliCaloPID.h"
  50 #include "AliAODCaloCluster.h"
  51 #include "AliVCaloCells.h"
  52 #include "AliVTrack.h"
  53 #include "AliAODPWG4Particle.h"
  54 #include "AliCalorimeterUtils.h"
  55 #include "AliVEvent.h"
  56
  57 // ---- Detector ----
  58 #include "AliEMCALPIDUtils.h"
  59
  60 ClassImp(AliCaloPID)
  61
  62
  63 //________________________
  64 AliCaloPID::AliCaloPID() :
  65 TObject(),                fDebug(-1),                  fParticleFlux(kLow),
  66 //Bayesian
  67 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
  68 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),
  69 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
  70 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),
  71 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
  72 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
  73 fPHOSPhotonWeightFormulaExpression(""),
  74 fPHOSPi0WeightFormulaExpression(""),
  75 //PID calculation
  76 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),
  77 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
  78 fTOFCut(0.),
  79 fPHOSDispersionCut(1000), fPHOSRCut(1000),
  80 //Split
  81 fUseSimpleMassCut(kFALSE),
  82 fUseSimpleM02Cut(kFALSE),
  83 fUseSplitAsyCut(kFALSE),
  84 fUseSplitSSCut(kTRUE),
  85 fSplitM02MaxCut(0),       fSplitM02MinCut(0),          fSplitMinNCells(0),
  86 fMassEtaMin(0),           fMassEtaMax(0),
  87 fMassPi0Min(0),           fMassPi0Max(0),
  88 fMassPhoMin(0),           fMassPhoMax(0),
  89 fM02MaxParamShiftNLMN(0),
  90 fSplitWidthSigma(0),      fMassShiftHighECell(0)
  91 {
  92   //Ctor
  93
  94   //Initialize parameters
  95   InitParameters();
  96 }
  97
  98 //________________________________________
  99 AliCaloPID::AliCaloPID(Int_t flux) :
 100 TObject(),                fDebug(-1),                  fParticleFlux(flux),
 101 //Bayesian
 102 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
 103 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),
 104 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
 105 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),
 106 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
 107 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
 108 fPHOSPhotonWeightFormulaExpression(""),
 109 fPHOSPi0WeightFormulaExpression(""),
 110 //PID calculation
 111 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),
 112 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
 113 fTOFCut(0.),
 114 fPHOSDispersionCut(1000), fPHOSRCut(1000),
 115 //Split
 116 fUseSimpleMassCut(kFALSE),
 117 fUseSimpleM02Cut(kFALSE),
 118 fUseSplitAsyCut(kFALSE),
 119 fUseSplitSSCut(kTRUE),
 120 fSplitM02MaxCut(0),       fSplitM02MinCut(0),          fSplitMinNCells(0),
 121 fMassEtaMin(0),           fMassEtaMax(0),
 122 fMassPi0Min(0),           fMassPi0Max(0),
 123 fMassPhoMin(0),           fMassPhoMax(0),
 124 fM02MaxParamShiftNLMN(0),
 125 fSplitWidthSigma(0),      fMassShiftHighECell(0)
 126 {
 127   //Ctor
 128
 129   //Initialize parameters
 130   InitParameters();
 131
 132 }
 133
 134 //_______________________________________________
 135 AliCaloPID::AliCaloPID(const TNamed * emcalpid) :
 136 TObject(),                   fDebug(-1),                  fParticleFlux(kLow),
 137 //Bayesian
 138 fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),
 139 fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
 140 fEMCALPhotonWeight(0.),      fEMCALPi0Weight(0.),
 141 fEMCALElectronWeight(0.),    fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
 142 fPHOSPhotonWeight(0.),       fPHOSPi0Weight(0.),
 143 fPHOSElectronWeight(0.),     fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
 144 fPHOSWeightFormula(0),       fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
 145 fPHOSPhotonWeightFormulaExpression(""),
 146 fPHOSPi0WeightFormulaExpression(""),
 147 //PID calculation
 148 fEMCALL0CutMax(100.),        fEMCALL0CutMin(0),
 149 fEMCALDEtaCut(2000.),        fEMCALDPhiCut(2000.),
 150 fTOFCut(0.),
 151 fPHOSDispersionCut(1000),    fPHOSRCut(1000),
 152 //Split
 153 fUseSimpleMassCut(kFALSE),
 154 fUseSimpleM02Cut(kFALSE),
 155 fUseSplitAsyCut(kFALSE),
 156 fUseSplitSSCut(kTRUE),
 157 fSplitM02MaxCut(0),       fSplitM02MinCut(0),          fSplitMinNCells(0),
 158 fMassEtaMin(0),           fMassEtaMax(0),
 159 fMassPi0Min(0),           fMassPi0Max(0),
 160 fMassPhoMin(0),           fMassPhoMax(0),
 161 fM02MaxParamShiftNLMN(0),
 162 fSplitWidthSigma(0),      fMassShiftHighECell(0)
 163
 164 {
 165   //Ctor
 166
 167   //Initialize parameters
 168   InitParameters();
 169 }
 170
 171 //_______________________
 172 AliCaloPID::~AliCaloPID()
 173 {
 174   //Dtor
 175
 176   delete fPHOSPhotonWeightFormula ;
 177   delete fPHOSPi0WeightFormula ;
 178   delete fEMCALPIDUtils ;
 179
 180 }
 181
 182 //_______________________________
 183 void AliCaloPID::InitParameters()
 184 {
 185   //Initialize the parameters of the PID.
 186
 187   // Bayesian
 188   fEMCALPhotonWeight   = 0.6 ;
 189   fEMCALPi0Weight      = 0.6 ;
 190   fEMCALElectronWeight = 0.6 ;
 191   fEMCALChargeWeight   = 0.6 ;
 192   fEMCALNeutralWeight  = 0.6 ;
 193
 194   fPHOSPhotonWeight    = 0.6 ;
 195   fPHOSPi0Weight       = 0.6 ;
 196   fPHOSElectronWeight  = 0.6 ;
 197   fPHOSChargeWeight    = 0.6 ;
 198   fPHOSNeutralWeight   = 0.6 ;
 199
 200   //Formula to set the PID weight threshold for photon or pi0
 201   fPHOSWeightFormula       = kFALSE;
 202   fPHOSPhotonWeightFormulaExpression = "0.98*(x<40)+ 0.68*(x>=100)+(x>=40 && x<100)*(0.98+x*(6e-3)-x*x*(2e-04)+x*x*x*(1.1e-06))";
 203   fPHOSPi0WeightFormulaExpression    = "0.98*(x<65)+ 0.915*(x>=100)+(x>=65 && x-x*(1.95e-3)-x*x*(4.31e-05)+x*x*x*(3.61e-07))"   ;
 204
 205   if(fRecalculateBayesian){
 206     if(fParticleFlux == kLow){
 207       printf("AliCaloPID::Init() - SetLOWFluxParam\n");
 208       fEMCALPIDUtils->SetLowFluxParam() ;
 209     }
 210     else if (fParticleFlux == kHigh){
 211       printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
 212       fEMCALPIDUtils->SetHighFluxParam() ;
 213     }
 214   }
 215
 216   //PID recalculation, not bayesian
 217
 218   //EMCAL
 219   fEMCALL0CutMax = 0.3 ;
 220   fEMCALL0CutMin = 0.01;
 221
 222   fEMCALDPhiCut  = 0.05; // Same cut as in AliEMCALRecoUtils
 223   fEMCALDEtaCut  = 0.025;// Same cut as in AliEMCALRecoUtils
 224
 225   // PHOS / EMCAL, not used
 226   fTOFCut        = 1.e-6;
 227
 228   //PHOS
 229   fPHOSRCut          = 2. ;
 230   fPHOSDispersionCut = 2.5;
 231
 232   // Cluster splitting
 233
 234   fSplitM02MinCut = 0.3 ;
 235   fSplitM02MaxCut = 5   ;
 236   fSplitMinNCells = 4   ;
 237
 238   fMassEtaMin  = 0.4;
 239   fMassEtaMax  = 0.6;
 240
 241   fMassPi0Min  = 0.11;
 242   fMassPi0Max  = 0.18;
 243
 244   fMassPhoMin  = 0.0;
 245   fMassPhoMax  = 0.08;
 246
 247   fMassPi0Param[0][0] = 0     ; // Constant term on mass dependence
 248   fMassPi0Param[0][1] = 0     ; // slope term on mass dependence
 249   fMassPi0Param[0][2] = 0     ; // E function change
 250   fMassPi0Param[0][3] = 0.044 ; // constant term on mass dependence
 251   fMassPi0Param[0][4] = 0.0049; // slope term on mass dependence
 252   fMassPi0Param[0][5] = 0.070 ; // Absolute low mass cut
 253
 254   fMassPi0Param[1][0] = 0.115 ; // Constant term below 21 GeV
 255   fMassPi0Param[1][1] = 0.00096; // slope term below 21 GeV
 256   fMassPi0Param[1][2] = 21    ; // E function change
 257   fMassPi0Param[1][3] = 0.10  ; // constant term on mass dependence
 258   fMassPi0Param[1][4] = 0.0017; // slope term on mass dependence
 259   fMassPi0Param[1][5] = 0.070 ; // Absolute low mass cut
 260
 261   fWidthPi0Param[0][0] = 0.012 ; // Constant term on width dependence
 262   fWidthPi0Param[0][1] = 0.0   ; // Slope term on width dependence
 263   fWidthPi0Param[0][2] = 19    ; // E function change
 264   fWidthPi0Param[0][3] = 0.0012; // Constant term on width dependence
 265   fWidthPi0Param[0][4] = 0.0006; // Slope term on width dependence
 266   fWidthPi0Param[0][5] = 0.0   ; // xx term
 267
 268   fWidthPi0Param[1][0] = 0.009 ; // Constant term on width dependence
 269   fWidthPi0Param[1][1] = 0.000 ; // Slope term on width dependence
 270   fWidthPi0Param[1][2] = 10    ; // E function change
 271   fWidthPi0Param[1][3] = 0.0023 ; // Constant term on width dependence
 272   fWidthPi0Param[1][4] = 0.00067; // Slope term on width dependence
 273   fWidthPi0Param[1][5] = 0.000 ;// xx term
 274
 275   fMassShiftHighECell = 0; // Shift of cuts in case of higher energy threshold in cells, 5 MeV when Ecell>150 MeV
 276
 277   //TF1 *lM02MinNLM1 = new TF1("M02MinNLM1","exp(2.135-0.245*x)",6,13.6);
 278   fM02MinParam[0][0] = 2.135  ;
 279   fM02MinParam[0][1] =-0.245  ;
 280   fM02MinParam[0][2] = 0.0    ;
 281   fM02MinParam[0][3] = 0.0    ;
 282   fM02MinParam[0][4] = 0.0    ;
 283
 284   // Same as NLM=1 for NLM=2
 285   fM02MinParam[1][0] = 2.135  ;
 286   fM02MinParam[1][1] =-0.245  ;
 287   fM02MinParam[1][2] = 0.0    ;
 288   fM02MinParam[1][3] = 0.0    ;
 289   fM02MinParam[1][4] = 0.0    ;
 290
 291   //TF1 *lM02MaxNLM1 = new TF1("M02MaxNLM1","exp(0.0662-0.0201*x)-0.0955+0.00186*x[0]+9.91/x[0]",6,100);
 292   fM02MaxParam[0][0] = 0.0662 ;
 293   fM02MaxParam[0][1] =-0.0201 ;
 294   fM02MaxParam[0][2] =-0.0955 ;
 295   fM02MaxParam[0][3] = 0.00186;
 296   fM02MaxParam[0][4] = 9.91   ;
 297
 298   //TF1 *lM02MaxNLM2 = new TF1("M02MaxNLM2","exp(0.353-0.0264*x)-0.524+0.00559*x[0]+21.9/x[0]",6,100);
 299   fM02MaxParam[1][0] = 0.353  ;
 300   fM02MaxParam[1][1] =-0.0264 ;
 301   fM02MaxParam[1][2] =-0.524  ;
 302   fM02MaxParam[1][3] = 0.00559;
 303   fM02MaxParam[1][4] = 21.9   ;
 304
 305   fM02MaxParamShiftNLMN = 0.75;
 306
 307   //TF1 *lAsyNLM1 = new TF1("lAsyNLM1","0.96-879/(x*x*x)",5,100);
 308   fAsyMinParam[0][0] = 0.96 ;
 309   fAsyMinParam[0][1] = 0    ;
 310   fAsyMinParam[0][2] =-879  ;
 311   fAsyMinParam[0][3] = 0.96 ; // Absolute max
 312
 313   //TF1 *lAsyNLM2 = new TF1("lAsyNLM2","0.95+0.0015*x-233/(x*x*x)",5,100);
 314   fAsyMinParam[1][0] = 0.95  ;
 315   fAsyMinParam[1][1] = 0.0015;
 316   fAsyMinParam[1][2] =-233   ;
 317   fAsyMinParam[1][3] = 1.0   ; // Absolute max
 318
 319   fSplitEFracMin[0]   = 0.0 ; // 0.96
 320   fSplitEFracMin[1]   = 0.0 ; // 0.96
 321   fSplitEFracMin[2]   = 0.0 ; // 0.7
 322
 323   fSubClusterEMin[0]  = 0.0; // 3 GeV
 324   fSubClusterEMin[1]  = 0.0; // 1 GeV
 325   fSubClusterEMin[2]  = 0.0; // 1 GeV
 326
 327
 328   fSplitWidthSigma = 3. ;
 329
 330 }
 331
 332
 333 //_________________________________________________________________________________________
 334 Bool_t AliCaloPID::IsInPi0SplitAsymmetryRange(Float_t energy, Float_t asy, Int_t nlm) const
 335 {
 336   // Select the appropriate mass range for pi0 selection in splitting method
 337   // No used yet in splitting ID decision
 338
 339   if(!fUseSplitAsyCut) return kTRUE ;
 340
 341   Float_t abasy = TMath::Abs(asy);
 342
 343   Int_t inlm = nlm-1;
 344   if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
 345
 346   // Get the parametrized min cut of asymmetry for NLM=2 up to 11 GeV
 347
 348   Float_t cut = fAsyMinParam[inlm][0] + fAsyMinParam[inlm][1]*energy + fAsyMinParam[inlm][2]/energy/energy/energy ;
 349
 350   // In any case and beyond validity energy range of the function,
 351   // the parameter cannot be smaller than 1
 352   if( cut > fAsyMinParam[inlm][3] ) cut = fAsyMinParam[inlm][3];
 353
 354   //printf("energy %2.2f - nlm: %d (%d)- p0 %f, p1 %f, p2 %f, p3 %f ; cut: %2.2f\n",energy,nlm,inlm,
 355   //       fAsyMinParam[inlm][0],fAsyMinParam[inlm][1],fAsyMinParam[inlm][2],fAsyMinParam[inlm][3],cut);
 356
 357   if(abasy < cut) return kTRUE;
 358   else            return kFALSE;
 359
 360 }
 361
 362 //______________________________________________________________________________________
 363 Bool_t AliCaloPID::IsInPi0SplitMassRange(Float_t energy, Float_t mass, Int_t nlm) const
 364 {
 365   // Select the appropriate mass range for pi0 selection in splitting method
 366
 367   if(fUseSimpleMassCut)
 368   {
 369     if(mass < fMassPi0Max && mass > fMassPi0Min) return kTRUE;
 370     else                                         return kFALSE;
 371   }
 372
 373   // Get the selected mean value as reference for the mass
 374   Int_t inlm = nlm-1;
 375   if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
 376
 377   Float_t meanMass = energy * fMassPi0Param[inlm][1] + fMassPi0Param[inlm][0];
 378   if(energy > fMassPi0Param[inlm][2]) meanMass = energy * fMassPi0Param[inlm][4] + fMassPi0Param[inlm][3];
 379
 380   // In case of higher energy cell cut than 50 MeV, smaller mean mass 0-5 MeV, not really necessary
 381   meanMass -= fMassShiftHighECell;
 382
 383   // Get the parametrized width of the mass
 384   Float_t width   = 0.009;
 385   if     (energy > 8 && energy < fWidthPi0Param[inlm][2])
 386     width = energy * fWidthPi0Param[inlm][1] + fWidthPi0Param[inlm][0];
 387   else if(              energy > fWidthPi0Param[inlm][2])
 388     width = energy * energy * fWidthPi0Param[inlm][5] + energy * fWidthPi0Param[inlm][4] + fWidthPi0Param[inlm][3];
 389
 390   // Calculate the 2 sigma cut
 391   Float_t minMass = meanMass-fSplitWidthSigma*width;
 392   Float_t maxMass = meanMass+fSplitWidthSigma*width;
 393
 394   // In case of low energy, hard cut to avoid conversions
 395   if(energy < 10  && minMass < fMassPi0Param[inlm][5] ) minMass = fMassPi0Param[inlm][5];
 396
 397   //printf("E %2.2f, mass %1.1f, nlm %d: sigma %1.1f width %3.1f, mean Mass %3.0f, minMass %3.0f, maxMass %3.0f\n ",
 398   //       energy,mass *1000, inlm, fSplitWidthSigma, width*1000, meanMass*1000,minMass*1000,maxMass*1000);
 399
 400   if(mass < maxMass && mass > minMass) return kTRUE;
 401   else                                 return kFALSE;
 402
 403 }
 404
 405 //________________________________________________
 406 Bool_t AliCaloPID::IsInM02Range(Float_t m02) const
 407 {
 408   // Select the appropriate m02 range, fix cut, not E dependent
 409
 410   Float_t minCut = fSplitM02MinCut;
 411   Float_t maxCut = fSplitM02MaxCut;
 412
 413   if(m02 < maxCut && m02 > minCut) return kTRUE;
 414   else                             return kFALSE;
 415
 416 }
 417
 418 //_______________________________________________________________________________
 419 Bool_t AliCaloPID::IsInPi0M02Range(Float_t energy, Float_t m02,  Int_t nlm) const
 420 {
 421   // Select the appropriate m02 range in splitting method for pi0
 422
 423   if(!fUseSplitSSCut) return kTRUE ;
 424
 425   //First check the absolute minimum and maximum
 426   if(!IsInM02Range(m02)) return kFALSE ;
 427
 428   //If requested, check the E dependent cuts
 429   else if(!fUseSimpleM02Cut)
 430   {
 431     Int_t inlm = nlm-1;
 432     if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
 433
 434     Float_t minCut = fSplitM02MinCut;
 435     Float_t maxCut = fSplitM02MaxCut;
 436
 437     //e^{a+bx} + c + dx + e/x
 438     if(energy > 1) minCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*energy ) +
 439                             fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*energy + fM02MinParam[inlm][4]/energy;
 440
 441     if(energy > 1) maxCut = TMath::Exp( fM02MaxParam[inlm][0] + fM02MaxParam[inlm][1]*energy ) +
 442                             fM02MaxParam[inlm][2] + fM02MaxParam[inlm][3]*energy + fM02MaxParam[inlm][4]/energy;
 443
 444     // In any case and beyond validity energy range of the function,
 445     // the parameter cannot be smaller than 0.3 or larger than 4-5
 446     if( minCut < fSplitM02MinCut) minCut = fSplitM02MinCut;
 447     if( maxCut > fSplitM02MaxCut) maxCut = fSplitM02MaxCut;
 448     if( nlm > 2 ) maxCut+=fM02MaxParamShiftNLMN;
 449
 450     //if(energy > 7) printf("\t \t E %2.2f, nlm %d, m02 %2.2f, minM02 %2.2f, maxM02 %2.2f\n",energy, nlm, m02,minCut,maxCut);
 451
 452     if(m02 < maxCut && m02 > minCut) return kTRUE;
 453     else                             return kFALSE;
 454
 455   }
 456
 457   else return kTRUE;
 458
 459 }
 460
 461
 462 //______________________________________________________________________________
 463 Bool_t AliCaloPID::IsInEtaM02Range(Float_t energy, Float_t m02, Int_t nlm) const
 464 {
 465   // Select the appropriate m02 range in splitting method to select eta's
 466   // Use same parametrization as pi0, just shift the distributions (to be tuned)
 467
 468   if(!fUseSplitSSCut) return kTRUE ;
 469
 470   //First check the absolute minimum and maximum
 471   if(!IsInM02Range(m02)) return kFALSE ;
 472
 473   //DO NOT USE, study parametrization
 474
 475   //If requested, check the E dependent cuts
 476   else if(!fUseSimpleM02Cut)
 477   {
 478     Int_t inlm = nlm-1;
 479     if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
 480
 481     Float_t minCut = fSplitM02MinCut;
 482     Float_t maxCut = fSplitM02MaxCut;
 483
 484     Float_t shiftE = energy-20; // to be tuned
 485     if(nlm==1) shiftE=energy-28;
 486
 487     //e^{a+bx} + c + dx + e/x
 488     if(shiftE > 1) minCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*shiftE ) +
 489                   fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*shiftE + fM02MinParam[inlm][4]/shiftE;
 490
 491     // In any case the parameter cannot be smaller than 0.3
 492     if( minCut < fSplitM02MinCut) minCut = fSplitM02MinCut;
 493
 494     shiftE = energy+20; // to be tuned
 495
 496     if(shiftE > 1)  maxCut = 1 + TMath::Exp( fM02MaxParam[inlm][0] + fM02MaxParam[inlm][1]*shiftE ) +
 497                              fM02MaxParam[inlm][2] + fM02MaxParam[inlm][3]*shiftE + fM02MaxParam[inlm][4]/shiftE;
 498
 499     // In any case the parameter cannot be smaller than 4-5
 500     if( maxCut > fSplitM02MaxCut) maxCut = fSplitM02MaxCut;
 501     if( nlm > 2 ) maxCut+=fM02MaxParamShiftNLMN;
 502
 503     //if(energy>6)printf("\t \t E %2.2f, nlm %d, m02 %2.2f, minM02 %2.2f, maxM02 %2.2f\n",energy, nlm, m02,minCut,maxCut);
 504
 505     if(m02 < maxCut && m02 > minCut) return kTRUE;
 506     else                             return kFALSE;
 507
 508   }
 509
 510   else return kTRUE;
 511
 512 }
 513
 514 //______________________________________________________________________________
 515 Bool_t AliCaloPID::IsInConM02Range(Float_t energy, Float_t m02, Int_t nlm) const
 516 {
 517   // Select the appropriate m02 range in splitting method for converted photons
 518   // Just min limit for pi0s is max for conversion.
 519
 520   if(!fUseSplitSSCut) return kTRUE ;
 521
 522   Float_t minCut = 0.1;
 523   Float_t maxCut = fSplitM02MinCut;
 524
 525   if(!fUseSimpleM02Cut)
 526   {
 527     Int_t inlm = nlm-1;
 528     if(nlm > 2) inlm=1; // only 2 cases defined nlm=1 and nlm>=2
 529
 530     //e^{a+bx} + c + dx + e/x
 531     if(energy > 1) maxCut = TMath::Exp( fM02MinParam[inlm][0] + fM02MinParam[inlm][1]*energy ) +
 532                             fM02MinParam[inlm][2] + fM02MinParam[inlm][3]*energy + fM02MinParam[inlm][4]/energy;
 533
 534     if( maxCut < fSplitM02MinCut) maxCut = fSplitM02MinCut;
 535   }
 536
 537   if(m02 < maxCut && m02 > minCut) return kTRUE;
 538   else                             return kFALSE;
 539
 540 }
 541
 542 //______________________________________________
 543 AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils()
 544 {
 545   // return pointer to AliEMCALPIDUtils, create it if needed
 546
 547   if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ;
 548   return fEMCALPIDUtils ;
 549
 550 }
 551
 552
 553 //________________________________________________________________
 554 Int_t AliCaloPID::GetIdentifiedParticleType(AliVCluster * cluster)
 555 {
 556   // Returns a PDG number corresponding to the likely ID of the cluster
 557
 558   Float_t energy  = cluster->E();
 559   Float_t lambda0 = cluster->GetM02();
 560   Float_t lambda1 = cluster->GetM20();
 561
 562   // ---------------------
 563   // Use bayesian approach
 564   // ---------------------
 565
 566   if(fUseBayesianWeights)
 567   {
 568     Double_t weights[AliPID::kSPECIESCN];
 569
 570     if(cluster->IsEMCAL() && fRecalculateBayesian)
 571     {
 572       fEMCALPIDUtils->ComputePID(energy, lambda0);
 573       for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
 574     }
 575     else
 576     {
 577       for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = cluster->GetPID()[i];
 578     }
 579
 580     if(fDebug > 0)  PrintClusterPIDWeights(weights);
 581
 582     return GetIdentifiedParticleTypeFromBayesWeights(cluster->IsEMCAL(), weights, energy);
 583   }
 584
 585   // -------------------------------------------------------
 586   // Calculate PID SS from data, do not use bayesian weights
 587   // -------------------------------------------------------
 588
 589   if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: EMCAL %d?, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
 590                         cluster->IsEMCAL(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
 591                         cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
 592
 593   if(cluster->IsEMCAL())
 594   {
 595     if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax);
 596
 597     if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
 598     else                                                     return kNeutralUnknown ;
 599   }    // EMCAL
 600   else // PHOS
 601   {
 602     if(TestPHOSDispersion(energy,lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
 603     else                                                                return kNeutralUnknown;
 604   }
 605
 606 }
 607
 608 //_________________________________________________________________________________________________________
 609 Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(Bool_t isEMCAL, Double_t * pid, Float_t energy)
 610 {
 611   //Return most probable identity of the particle after bayesian weights calculated in reconstruction
 612
 613   if(!pid)
 614   {
 615     AliFatal("pid pointer not initialized!!!");
 616     return kNeutralUnknown; // not needed, added to content coverity
 617   }
 618
 619   Float_t wPh  =  fPHOSPhotonWeight ;
 620   Float_t wPi0 =  fPHOSPi0Weight ;
 621   Float_t wE   =  fPHOSElectronWeight ;
 622   Float_t wCh  =  fPHOSChargeWeight ;
 623   Float_t wNe  =  fPHOSNeutralWeight ;
 624
 625   if(!isEMCAL && fPHOSWeightFormula){
 626     wPh  = GetPHOSPhotonWeightFormula()->Eval(energy) ;
 627     wPi0 = GetPHOSPi0WeightFormula()   ->Eval(energy);
 628   }
 629   else
 630   {
 631     wPh  =  fEMCALPhotonWeight ;
 632     wPi0 =  fEMCALPi0Weight ;
 633     wE   =  fEMCALElectronWeight ;
 634     wCh  =  fEMCALChargeWeight ;
 635     wNe  =  fEMCALNeutralWeight ;
 636   }
 637
 638   if(fDebug > 0) PrintClusterPIDWeights(pid);
 639
 640   Int_t pdg = kNeutralUnknown ;
 641   Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
 642   pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
 643   Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
 644   Float_t allChargedWeight    = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
 645   Float_t allNeutralWeight    = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
 646
 647   //Select most probable ID
 648   if(!isEMCAL) // PHOS
 649   {
 650     if(pid[AliVCluster::kPhoton] > wPh)        pdg = kPhoton ;
 651     else if(pid[AliVCluster::kPi0] > wPi0)     pdg = kPi0 ;
 652     else if(pid[AliVCluster::kElectron] > wE)  pdg = kElectron ;
 653     else if(pid[AliVCluster::kEleCon] >  wE)   pdg = kEleCon ;
 654     else if(chargedHadronWeight > wCh)         pdg = kChargedHadron ;
 655     else if(neutralHadronWeight > wNe)         pdg = kNeutralHadron ;
 656     else if(allChargedWeight >  allNeutralWeight)
 657       pdg = kChargedUnknown ;
 658     else
 659       pdg = kNeutralUnknown ;
 660   }
 661   else //EMCAL
 662   {
 663     if(pid[AliVCluster::kPhoton]  > wPh)                     pdg = kPhoton ;
 664     else if(pid[AliVCluster::kElectron]  > wE)               pdg = kElectron ;
 665     else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
 666     else if(pid[AliVCluster::kPi0] > wPi0)                   pdg = kPi0 ;
 667     else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
 668     else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
 669     else                                                     pdg = kNeutralUnknown ;
 670   }
 671
 672   if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
 673
 674   return pdg ;
 675
 676 }
 677
 678 //____________________________________________________________________________________________________________
 679 Int_t AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting(AliVCluster* cluster,
 680                                                                 AliVCaloCells* cells,
 681                                                                 AliCalorimeterUtils * caloutils,
 682                                                                 Double_t   vertex[3],
 683                                                                 Int_t    & nMax,
 684                                                                 Double_t & mass,   Double_t & angle,
 685                                                                 TLorentzVector & l1, TLorentzVector & l2,
 686                                                                 Int_t   & absId1,   Int_t   & absId2,
 687                                                                 Float_t & distbad1, Float_t & distbad2,
 688                                                                 Bool_t  & fidcut1,  Bool_t  & fidcut2  ) const
 689 {
 690   // Split the cluster in 2, do invariant mass, get the mass and decide
 691   // if this is a photon, pi0, eta, ...
 692
 693   Float_t eClus  = cluster->E();
 694   Float_t m02    = cluster->GetM02();
 695   const Int_t nc = cluster->GetNCells();
 696   Int_t   absIdList[nc];
 697   Float_t maxEList [nc];
 698
 699   mass  = -1.;
 700   angle = -1.;
 701
 702   //If too low number of cells, skip it
 703   if ( nc < fSplitMinNCells)  return kNeutralUnknown ;
 704
 705   if(fDebug > 0) printf("\t pass nCells cut\n");
 706
 707   // Get Number of local maxima
 708   nMax  = caloutils->GetNumberOfLocalMaxima(cluster, cells, absIdList, maxEList) ;
 709
 710   if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting() - Cluster : E %1.1f, M02 %1.2f, NLM %d, N Cells %d\n",
 711                         eClus,m02,nMax,nc);
 712
 713   //---------------------------------------------------------------------
 714   // Get the 2 max indeces and do inv mass
 715   //---------------------------------------------------------------------
 716
 717   TString  calorimeter = "EMCAL";
 718   if(cluster->IsPHOS()) calorimeter = "PHOS";
 719
 720   if     ( nMax == 2 )
 721   {
 722     absId1 = absIdList[0];
 723     absId2 = absIdList[1];
 724
 725     //Order in energy
 726     Float_t en1 = cells->GetCellAmplitude(absId1);
 727     caloutils->RecalibrateCellAmplitude(en1,calorimeter,absId1);
 728     Float_t en2 = cells->GetCellAmplitude(absId2);
 729     caloutils->RecalibrateCellAmplitude(en2,calorimeter,absId2);
 730     if(en1 < en2)
 731     {
 732       absId2 = absIdList[0];
 733       absId1 = absIdList[1];
 734     }
 735   }
 736   else if( nMax == 1 )
 737   {
 738
 739     absId1 = absIdList[0];
 740
 741     //Find second highest energy cell
 742
 743     Float_t enmax = 0 ;
 744     for(Int_t iDigit = 0 ; iDigit < cluster->GetNCells() ; iDigit++)
 745     {
 746       Int_t absId = cluster->GetCellsAbsId()[iDigit];
 747       if( absId == absId1 ) continue ;
 748       Float_t endig = cells->GetCellAmplitude(absId);
 749       caloutils->RecalibrateCellAmplitude(endig,calorimeter,absId);
 750       if(endig > enmax)
 751       {
 752         enmax  = endig ;
 753         absId2 = absId ;
 754       }
 755     }// cell loop
 756   }// 1 maxima
 757   else
 758   {  // n max > 2
 759     // loop on maxima, find 2 highest
 760
 761     // First max
 762     Float_t enmax = 0 ;
 763     for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
 764     {
 765       Float_t endig = maxEList[iDigit];
 766       if(endig > enmax)
 767       {
 768         enmax  = endig ;
 769         absId1 = absIdList[iDigit];
 770       }
 771     }// first maxima loop
 772
 773     // Second max
 774     Float_t enmax2 = 0;
 775     for(Int_t iDigit = 0 ; iDigit < nMax ; iDigit++)
 776     {
 777       if(absIdList[iDigit]==absId1) continue;
 778       Float_t endig = maxEList[iDigit];
 779       if(endig > enmax2)
 780       {
 781         enmax2  = endig ;
 782         absId2 = absIdList[iDigit];
 783       }
 784     }// second maxima loop
 785
 786   } // n local maxima > 2
 787
 788   if(absId2<0 || absId1<0)
 789   {
 790     if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleTypeFromClusterSplitting() - Bad index for local maxima : N max %d, i1 %d, i2 %d, cluster E %2.2f, ncells %d, m02 %2.2f\n",
 791                           nMax,absId1,absId2,eClus,nc,m02);
 792     return kNeutralUnknown ;
 793   }
 794
 795   //---------------------------------------------------------------------
 796   // Split the cluster energy in 2, around the highest 2 local maxima
 797   //---------------------------------------------------------------------
 798
 799   AliAODCaloCluster cluster1(0, 0,NULL,0.,NULL,NULL,1,0);
 800   AliAODCaloCluster cluster2(1, 0,NULL,0.,NULL,NULL,1,0);
 801
 802   caloutils->SplitEnergy(absId1,absId2,cluster, cells, &cluster1, &cluster2,nMax); /*absIdList, maxEList,*/
 803
 804   fidcut1 = caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), &cluster1,cells);
 805   fidcut2 = caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), &cluster2,cells);
 806
 807   caloutils->GetEMCALRecoUtils()->RecalculateClusterDistanceToBadChannel(caloutils->GetEMCALGeometry(),cells,&cluster1);
 808   caloutils->GetEMCALRecoUtils()->RecalculateClusterDistanceToBadChannel(caloutils->GetEMCALGeometry(),cells,&cluster2);
 809
 810   distbad1 = cluster1.GetDistanceToBadChannel();
 811   distbad2 = cluster2.GetDistanceToBadChannel();
 812 //  if(!fidcut2 || !fidcut1 || distbad1 < 2 || distbad2 < 2)
 813 //    printf("*** Dist to bad channel cl %f, cl1 %f, cl2 %f; fid cut cl %d, cl1 %d, cl2 %d \n",
 814 //           cluster->GetDistanceToBadChannel(),distbad1,distbad2,
 815 //           caloutils->GetEMCALRecoUtils()->CheckCellFiducialRegion(caloutils->GetEMCALGeometry(), cluster,cells),fidcut1,fidcut2);
 816
 817   cluster1.GetMomentum(l1,vertex);
 818   cluster2.GetMomentum(l2,vertex);
 819
 820   mass  = (l1+l2).M();
 821   angle = l2.Angle(l1.Vect());
 822   Float_t e1 = cluster1.E();
 823   Float_t e2 = cluster2.E();
 824
 825   // Consider clusters with splitted energy not too different to original cluster energy
 826   Float_t splitFracCut = 0;
 827   if(nMax < 3)  splitFracCut = fSplitEFracMin[nMax-1];
 828   else          splitFracCut = fSplitEFracMin[2];
 829   if((e1+e2)/eClus < splitFracCut) return kNeutralUnknown ;
 830
 831   if(fDebug > 0) printf("\t pass Split E frac cut\n");
 832
 833   // Consider sub-clusters with minimum energy
 834   Float_t minECut = fSubClusterEMin[2];
 835   if     (nMax == 2)  minECut = fSubClusterEMin[1];
 836   else if(nMax == 1)  minECut = fSubClusterEMin[0];
 837   if(e1 < minECut || e2 < minECut)
 838   {
 839     //printf("Reject: e1 %2.1f, e2 %2.1f, cut %2.1f\n",e1,e2,minECut);
 840     return kNeutralUnknown ;
 841   }
 842
 843   if(fDebug > 0) printf("\t pass min sub-cluster E cut\n");
 844
 845   // Asymmetry of cluster
 846   Float_t asy =-10;
 847   if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
 848
 849   if( !IsInPi0SplitAsymmetryRange(eClus,asy,nMax) ) return kNeutralUnknown ;
 850
 851
 852   if (fDebug>0) printf("\t pass asymmetry cut\n");
 853
 854   Bool_t pi0OK = kFALSE;
 855   Bool_t etaOK = kFALSE;
 856   Bool_t conOK = kFALSE;
 857
 858   //If too small or big M02, skip it
 859   if     (IsInPi0M02Range(eClus,m02,nMax))  pi0OK = kTRUE;
 860   else if(IsInEtaM02Range(eClus,m02,nMax))  etaOK = kTRUE;
 861   else if(IsInConM02Range(eClus,m02,nMax))  conOK = kTRUE;
 862
 863   Float_t energy = eClus;
 864   if(nMax > 2) energy = e1+e2; // In case of NLM>2 use mass cut for NLM=2 but for the split sum not the cluster energy that is not the pi0 E.
 865
 866   // Check the mass, and set an ID to the splitted cluster
 867   if     ( conOK && mass < fMassPhoMax && mass > fMassPhoMin ) { if(fDebug > 0) printf("\t Split Conv \n"); return kPhoton ; }
 868   else if( etaOK && mass < fMassEtaMax && mass > fMassEtaMin ) { if(fDebug > 0) printf("\t Split Eta \n");  return kEta    ; }
 869   else if( pi0OK && IsInPi0SplitMassRange(energy,mass,nMax)  ) { if(fDebug > 0) printf("\t Split Pi0 \n");  return kPi0    ; }
 870   else                                                                                                      return kNeutralUnknown ;
 871
 872 }
 873
 874 //_________________________________________
 875 TString  AliCaloPID::GetPIDParametersList()
 876 {
 877   //Put data member values in string to keep in output container
 878
 879   TString parList ; //this will be list of parameters used for this analysis.
 880   const Int_t buffersize = 255;
 881   char onePar[buffersize] ;
 882   snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ;
 883   parList+=onePar ;
 884   if(fUseBayesianWeights)
 885   {
 886     snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
 887     parList+=onePar ;
 888     snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
 889     parList+=onePar ;
 890     snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
 891     parList+=onePar ;
 892     snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
 893     parList+=onePar ;
 894     snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
 895     parList+=onePar ;
 896     snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
 897     parList+=onePar ;
 898     snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
 899     parList+=onePar ;
 900     snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
 901     parList+=onePar ;
 902     snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
 903     parList+=onePar ;
 904     snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
 905     parList+=onePar ;
 906
 907     if(fPHOSWeightFormula)
 908     {
 909       snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ;
 910       parList+=onePar;
 911       snprintf(onePar,buffersize,"PHOS Pi0    Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data()    ) ;
 912       parList+=onePar;
 913     }
 914   }
 915   else
 916   {
 917     snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f  (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ;
 918     parList+=onePar ;
 919     snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f  (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ;
 920     parList+=onePar ;
 921     snprintf(onePar,buffersize,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
 922     parList+=onePar ;
 923     snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f  (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ;
 924     parList+=onePar ;
 925
 926   }
 927
 928   if(fUseSimpleM02Cut)
 929   {
 930     snprintf(onePar,buffersize,"%2.2f< M02 < %2.2f \n",    fSplitM02MinCut, fSplitM02MaxCut) ;
 931     parList+=onePar ;
 932   }
 933   snprintf(onePar,buffersize,"fMinNCells =%d \n",        fSplitMinNCells) ;
 934   parList+=onePar ;
 935   if(fUseSimpleMassCut)
 936   {
 937     snprintf(onePar,buffersize,"pi0 : %2.1f < m <%2.1f\n", fMassPi0Min,fMassPi0Max);
 938     parList+=onePar ;
 939   }
 940   snprintf(onePar,buffersize,"eta : %2.1f < m <%2.1f\n", fMassEtaMin,fMassEtaMax);
 941   parList+=onePar ;
 942   snprintf(onePar,buffersize,"conv: %2.1f < m <%2.1f\n", fMassPhoMin,fMassPhoMax);
 943   parList+=onePar ;
 944
 945
 946   return parList;
 947
 948 }
 949
 950 //________________________________________________
 951 void AliCaloPID::Print(const Option_t * opt) const
 952 {
 953
 954   //Print some relevant parameters set for the analysis
 955   if(! opt)
 956     return;
 957
 958   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 959
 960   if(fUseBayesianWeights)
 961   {
 962     printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
 963            fPHOSPhotonWeight,   fPHOSPi0Weight,
 964            fPHOSElectronWeight, fPHOSChargeWeight, fPHOSNeutralWeight) ;
 965     printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
 966            fEMCALPhotonWeight,   fEMCALPi0Weight,
 967            fEMCALElectronWeight, fEMCALChargeWeight, fEMCALNeutralWeight) ;
 968
 969     printf("PHOS Parametrized weight on?  =     %d\n",  fPHOSWeightFormula) ;
 970     if(fPHOSWeightFormula)
 971     {
 972       printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
 973       printf("Pi0    weight formula = %s\n", fPHOSPi0WeightFormulaExpression   .Data());
 974     }
 975     if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux?    = %d\n",fParticleFlux);
 976   }
 977   else
 978   {
 979     printf("TOF cut        = %e\n",                                   fTOFCut);
 980     printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n",            fEMCALL0CutMin,fEMCALL0CutMax);
 981     printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n",        fEMCALDEtaCut, fEMCALDPhiCut);
 982     printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut,     fPHOSDispersionCut) ;
 983
 984   }
 985
 986   printf("Min. N Cells =%d \n",         fSplitMinNCells) ;
 987   if(fUseSimpleM02Cut) printf("%2.2f < lambda_0^2 <%2.2f \n",fSplitM02MinCut,fSplitM02MaxCut);
 988   if(fUseSimpleMassCut)printf("pi0 : %2.2f<m<%2.2f \n",      fMassPi0Min,fMassPi0Max);
 989   printf("eta : %2.2f<m<%2.2f \n",      fMassEtaMin,fMassEtaMax);
 990   printf("phot: %2.2f<m<%2.2f \n",      fMassPhoMin,fMassPhoMax);
 991
 992   printf(" \n");
 993
 994 }
 995
 996 //_________________________________________________________________
 997 void AliCaloPID::PrintClusterPIDWeights(const Double_t * pid) const
 998 {
 999   // print PID of cluster, (AliVCluster*)cluster->GetPID()
1000
1001   printf("AliCaloPID::PrintClusterPIDWeights() \n \t ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, \n \t \
1002          pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
1003          pid[AliVCluster::kPhoton],    pid[AliVCluster::kPi0],
1004          pid[AliVCluster::kElectron],  pid[AliVCluster::kEleCon],
1005          pid[AliVCluster::kPion],      pid[AliVCluster::kKaon],
1006          pid[AliVCluster::kProton],
1007          pid[AliVCluster::kNeutron],   pid[AliVCluster::kKaon0]);
1008
1009 }
1010
1011 //___________________________________________________________________________
1012 void AliCaloPID::SetPIDBits(AliVCluster * cluster,
1013                             AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
1014                             AliVEvent* event)
1015 {
1016   //Set Bits for PID selection
1017
1018   //Dispersion/lambdas
1019   //Double_t disp= cluster->GetDispersion()  ;
1020   Double_t l1  = cluster->GetM20() ;
1021   Double_t l0  = cluster->GetM02() ;
1022   Bool_t isDispOK = kTRUE ;
1023   if(cluster->IsPHOS()){
1024     if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
1025     else                                                        isDispOK = kFALSE;
1026   }
1027   else{//EMCAL
1028
1029     if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
1030
1031   }
1032
1033   ph->SetDispBit(isDispOK) ;
1034
1035   //TOF
1036   Double_t tof=cluster->GetTOF()  ;
1037   ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
1038
1039   //Charged
1040   Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
1041
1042   ph->SetChargedBit(isNeutral);
1043
1044   //Set PID pdg
1045   ph->SetIdentifiedParticleType(GetIdentifiedParticleType(cluster));
1046
1047   if(fDebug > 0)
1048   {
1049     printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);
1050     printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
1051            ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
1052   }
1053 }
1054
1055 //_________________________________________________________
1056 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
1057                                   AliCalorimeterUtils * cu,
1058                                   AliVEvent* event) const
1059 {
1060   //Check if there is any track attached to this cluster
1061
1062   Int_t nMatches = cluster->GetNTracksMatched();
1063   AliVTrack * track = 0;
1064   Double_t p[3];
1065
1066   if(nMatches > 0)
1067   {
1068     //In case of ESDs, by default without match one entry with negative index, no match, reject.
1069     if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
1070     {
1071       Int_t iESDtrack = cluster->GetTrackMatchedIndex();
1072       if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
1073       else return kFALSE;
1074
1075       if (!track)
1076       {
1077         if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n");
1078         return kFALSE;
1079       }
1080     }
1081     else { // AOD
1082       track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
1083       if (!track)
1084       {
1085         if(fDebug > 0) printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n");
1086         return kFALSE;
1087       }
1088     }
1089
1090     Float_t dZ  = cluster->GetTrackDz();
1091     Float_t dR  = cluster->GetTrackDx();
1092
1093     // if track matching was recalculated
1094     if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn())
1095     {
1096       dR = 2000., dZ = 2000.;
1097       cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
1098     }
1099
1100     if(cluster->IsPHOS())
1101     {
1102       track->GetPxPyPz(p) ;
1103       TLorentzVector trackmom(p[0],p[1],p[2],0);
1104       Int_t charge = track->Charge();
1105       Double_t mf  = event->GetMagneticField();
1106       if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
1107       else                                                                    return kFALSE;
1108
1109     }//PHOS
1110     else //EMCAL
1111     {
1112     if(fDebug > 1)
1113         printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut);
1114
1115       if(TMath::Abs(dR) < fEMCALDPhiCut &&
1116          TMath::Abs(dZ) < fEMCALDEtaCut)   return kTRUE;
1117       else                                 return kFALSE;
1118
1119     }//EMCAL cluster
1120
1121
1122   } // more than 1 match, at least one track in array
1123   else return kFALSE;
1124
1125 }
1126
1127 //___________________________________________________________________________________________________
1128 Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const
1129 {
1130   //Check if cluster photon-like. Uses photon cluster parameterization in real pp data
1131   //Returns distance in sigmas. Recommended cut 2.5
1132
1133   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1134   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1135   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1136   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1137   Double_t c       =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1138   Double_t r2      = 0.5*  (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1139                      0.5*  (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1140                      0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1141
1142   if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut);
1143
1144   return TMath::Sqrt(r2) ;
1145
1146 }
1147
1148 //_______________________________________________________________________________________________
1149 Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx,  const Double_t dz, const Double_t pt,
1150                                         const Int_t charge, const Double_t mf) const
1151 {
1152   //Checks distance to the closest track. Takes into account
1153   //non-perpendicular incidence of tracks.
1154   //returns distance in sigmas. Recommended cut: 2.
1155   //Requires (sign) of magnetic filed. onc can find it for example as following
1156   //  Double_t mf=0. ;
1157   //  AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
1158   //  if(event)
1159   //    mf = event->GetMagneticField(); //Positive for ++ and negative for --
1160
1161
1162   Double_t meanX = 0.;
1163   Double_t meanZ = 0.;
1164   Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1165                            6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
1166                            1.59219);
1167   Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
1168                            1.60) ;
1169
1170   if(mf<0.){ //field --
1171     meanZ = -0.468318 ;
1172     if(charge>0)
1173       meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
1174                          0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1175     else
1176       meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
1177                          1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1178   }
1179   else{ //Field ++
1180     meanZ = -0.468318;
1181     if(charge>0)
1182       meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
1183                          0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1184     else
1185       meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
1186                          0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1187   }
1188
1189   Double_t rz = (dz-meanZ)/sz ;
1190   Double_t rx = (dx-meanX)/sx ;
1191
1192   if(fDebug > 0)
1193     printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut);
1194
1195   return TMath::Sqrt(rx*rx+rz*rz) ;
1196
1197 }