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