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