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