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