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 TString calo, const TLorentzVector mom, 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 TString calo, 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 TString calo, const TLorentzVector mom, 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 "AliVCluster.h"
  51 #include "AliVTrack.h"
  52 #include "AliAODPWG4Particle.h"
  53 #include "AliCalorimeterUtils.h"
  54 #include "AliVEvent.h"
  55
  56 // ---- Detector ----
  57 #include "AliEMCALPIDUtils.h"
  58
  59 ClassImp(AliCaloPID)
  60
  61
  62 //________________________
  63 AliCaloPID::AliCaloPID() :
  64 TObject(),                fDebug(-1),                  fParticleFlux(kLow),
  65 //Bayesian
  66 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
  67 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),
  68 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
  69 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),
  70 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
  71 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
  72 fPHOSPhotonWeightFormulaExpression(""),
  73 fPHOSPi0WeightFormulaExpression(""),
  74 //PID calculation
  75 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),
  76 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
  77 fTOFCut(0.),
  78 fPHOSDispersionCut(1000), fPHOSRCut(1000)
  79 {
  80   //Ctor
  81
  82   //Initialize parameters
  83   InitParameters();
  84 }
  85
  86 //________________________________________
  87 AliCaloPID::AliCaloPID(const Int_t flux) :
  88 TObject(),                fDebug(-1),                  fParticleFlux(flux),
  89 //Bayesian
  90 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
  91 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),
  92 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
  93 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),
  94 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
  95 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
  96 fPHOSPhotonWeightFormulaExpression(""),
  97 fPHOSPi0WeightFormulaExpression(""),
  98 //PID calculation
  99 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),
 100 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
 101 fTOFCut(0.),
 102 fPHOSDispersionCut(1000), fPHOSRCut(1000)
 103 {
 104   //Ctor
 105
 106   //Initialize parameters
 107   InitParameters();
 108
 109 }
 110
 111 //_______________________________________________
 112 AliCaloPID::AliCaloPID(const TNamed * emcalpid) :
 113 TObject(),                   fDebug(-1),                  fParticleFlux(kLow),
 114 //Bayesian
 115 fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),
 116 fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
 117 fEMCALPhotonWeight(0.),      fEMCALPi0Weight(0.),
 118 fEMCALElectronWeight(0.),    fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
 119 fPHOSPhotonWeight(0.),       fPHOSPi0Weight(0.),
 120 fPHOSElectronWeight(0.),     fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
 121 fPHOSWeightFormula(0),       fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
 122 fPHOSPhotonWeightFormulaExpression(""),
 123 fPHOSPi0WeightFormulaExpression(""),
 124 //PID calculation
 125 fEMCALL0CutMax(100.),        fEMCALL0CutMin(0),
 126 fEMCALDEtaCut(2000.),        fEMCALDPhiCut(2000.),
 127 fTOFCut(0.),
 128 fPHOSDispersionCut(1000),    fPHOSRCut(1000)
 129 {
 130   //Ctor
 131
 132   //Initialize parameters
 133   InitParameters();
 134 }
 135
 136 //_______________________
 137 AliCaloPID::~AliCaloPID()
 138 {
 139   //Dtor
 140
 141   delete fPHOSPhotonWeightFormula ;
 142   delete fPHOSPi0WeightFormula ;
 143   delete fEMCALPIDUtils ;
 144
 145 }
 146
 147 //_______________________________
 148 void AliCaloPID::InitParameters()
 149 {
 150   //Initialize the parameters of the PID.
 151
 152   // Bayesian
 153   fEMCALPhotonWeight   = 0.6 ;
 154   fEMCALPi0Weight      = 0.6 ;
 155   fEMCALElectronWeight = 0.6 ;
 156   fEMCALChargeWeight   = 0.6 ;
 157   fEMCALNeutralWeight  = 0.6 ;
 158
 159   fPHOSPhotonWeight    = 0.6 ;
 160   fPHOSPi0Weight       = 0.6 ;
 161   fPHOSElectronWeight  = 0.6 ;
 162   fPHOSChargeWeight    = 0.6 ;
 163   fPHOSNeutralWeight   = 0.6 ;
 164
 165   //Formula to set the PID weight threshold for photon or pi0
 166   fPHOSWeightFormula       = kFALSE;
 167   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))";
 168   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))"   ;
 169
 170   if(fRecalculateBayesian){
 171     if(fParticleFlux == kLow){
 172       printf("AliCaloPID::Init() - SetLOWFluxParam\n");
 173       fEMCALPIDUtils->SetLowFluxParam() ;
 174     }
 175     else if (fParticleFlux == kHigh){
 176       printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
 177       fEMCALPIDUtils->SetHighFluxParam() ;
 178     }
 179   }
 180
 181   //PID recalculation, not bayesian
 182
 183   //EMCAL
 184   fEMCALL0CutMax = 0.3 ;
 185   fEMCALL0CutMin = 0.01;
 186
 187   fEMCALDPhiCut  = 0.05; // Same cut as in AliEMCALRecoUtils
 188   fEMCALDEtaCut  = 0.025;// Same cut as in AliEMCALRecoUtils
 189
 190   // PHOS / EMCAL, not used
 191   fTOFCut        = 1.e-6;
 192
 193   //PHOS
 194   fPHOSRCut          = 2. ;
 195   fPHOSDispersionCut = 2.5;
 196
 197 }
 198
 199 //______________________________________________
 200 AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils()
 201 {
 202   // return pointer to AliEMCALPIDUtils, create it if needed
 203
 204   if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ;
 205   return fEMCALPIDUtils ;
 206
 207 }
 208
 209
 210 //______________________________________________________________________
 211 Int_t AliCaloPID::GetIdentifiedParticleType(const TString calo,
 212                                             const TLorentzVector mom,
 213                                             const AliVCluster * cluster)
 214 {
 215   // Returns a PDG number corresponding to the likely ID of the cluster
 216
 217   Float_t energy  = mom.E();
 218   Float_t lambda0 = cluster->GetM02();
 219   Float_t lambda1 = cluster->GetM20();
 220
 221   // ---------------------
 222   // Use bayesian approach
 223   // ---------------------
 224
 225   if(fUseBayesianWeights){
 226
 227     Double_t weights[AliPID::kSPECIESN];
 228
 229     if(calo == "EMCAL"&& fRecalculateBayesian){
 230       fEMCALPIDUtils->ComputePID(energy, lambda0);
 231       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
 232     }
 233     else {
 234       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i];
 235     }
 236
 237     if(fDebug > 0)  {
 238       printf("AliCaloPID::GetIdentifiedParticleType: BEFORE calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
 239              calo.Data(),
 240              weights[AliVCluster::kPhoton],    weights[AliVCluster::kPi0],
 241              weights[AliVCluster::kElectron],  weights[AliVCluster::kEleCon],
 242              weights[AliVCluster::kPion],      weights[AliVCluster::kKaon],
 243              weights[AliVCluster::kProton],
 244              weights[AliVCluster::kNeutron],   weights[AliVCluster::kKaon0]);
 245     }
 246
 247     return GetIdentifiedParticleTypeFromBayesWeights(calo, weights, energy);
 248   }
 249
 250   // -------------------------------------------------------
 251   // Calculate PID SS from data, do not use bayesian weights
 252   // -------------------------------------------------------
 253
 254   if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
 255                         calo.Data(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
 256                         cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
 257
 258   if(cluster->IsEMCAL()){
 259
 260     if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax);
 261
 262     if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
 263     else                                                     return kNeutralUnknown ;
 264   }//EMCAL
 265   else {//PHOS
 266     if(TestPHOSDispersion(mom.Pt(),lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
 267     else                                                                  return kNeutralUnknown;
 268   }
 269
 270 }
 271
 272 //_______________________________________________________________________________
 273 Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const TString calo,
 274                                                             const Double_t * pid,
 275                                                             const Float_t energy)
 276 {
 277   //Return most probable identity of the particle after bayesian weights calculated in reconstruction
 278
 279   if(!pid){
 280     printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n");
 281     abort();
 282   }
 283
 284   Float_t wPh  =  fPHOSPhotonWeight ;
 285   Float_t wPi0 =  fPHOSPi0Weight ;
 286   Float_t wE   =  fPHOSElectronWeight ;
 287   Float_t wCh  =  fPHOSChargeWeight ;
 288   Float_t wNe  =  fPHOSNeutralWeight ;
 289
 290   if(calo == "PHOS" && fPHOSWeightFormula){
 291     wPh  = GetPHOSPhotonWeightFormula()->Eval(energy) ;
 292     wPi0 = GetPHOSPi0WeightFormula()   ->Eval(energy);
 293   }
 294
 295   if(calo == "EMCAL"){
 296
 297     wPh  =  fEMCALPhotonWeight ;
 298     wPi0 =  fEMCALPi0Weight ;
 299     wE   =  fEMCALElectronWeight ;
 300     wCh  =  fEMCALChargeWeight ;
 301     wNe  =  fEMCALNeutralWeight ;
 302
 303   }
 304
 305   if(fDebug > 0)  printf("AliCaloPID::GetIdentifiedParticleType: calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
 306                          calo.Data(),pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
 307                          pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
 308                          pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton],
 309                          pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
 310
 311   Int_t pdg = kNeutralUnknown ;
 312   Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
 313   pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
 314   Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
 315   Float_t allChargedWeight    = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
 316   Float_t allNeutralWeight    = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
 317
 318   //Select most probable ID
 319   if(calo=="PHOS"){
 320     if(pid[AliVCluster::kPhoton] > wPh)        pdg = kPhoton ;
 321     else if(pid[AliVCluster::kPi0] > wPi0)     pdg = kPi0 ;
 322     else if(pid[AliVCluster::kElectron] > wE)  pdg = kElectron ;
 323     else if(pid[AliVCluster::kEleCon] >  wE)   pdg = kEleCon ;
 324     else if(chargedHadronWeight > wCh)         pdg = kChargedHadron ;
 325     else if(neutralHadronWeight > wNe)         pdg = kNeutralHadron ;
 326     else if(allChargedWeight >  allNeutralWeight)
 327       pdg = kChargedUnknown ;
 328     else
 329       pdg = kNeutralUnknown ;
 330   }
 331   else{//EMCAL
 332
 333     if(pid[AliVCluster::kPhoton]  > wPh)                     pdg = kPhoton ;
 334     else if(pid[AliVCluster::kElectron]  > wE)               pdg = kElectron ;
 335     else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
 336     else if(pid[AliVCluster::kPi0] > wPi0)                   pdg = kPi0 ;
 337     else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
 338     else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
 339     else                                                     pdg = kNeutralUnknown ;
 340   }
 341
 342   if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
 343
 344   return pdg ;
 345
 346 }
 347
 348 //_________________________________________
 349 TString  AliCaloPID::GetPIDParametersList()
 350 {
 351   //Put data member values in string to keep in output container
 352
 353   TString parList ; //this will be list of parameters used for this analysis.
 354   const Int_t buffersize = 255;
 355   char onePar[buffersize] ;
 356   snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ;
 357   parList+=onePar ;
 358   if(fUseBayesianWeights){
 359     snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
 360     parList+=onePar ;
 361     snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
 362     parList+=onePar ;
 363     snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
 364     parList+=onePar ;
 365     snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
 366     parList+=onePar ;
 367     snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
 368     parList+=onePar ;
 369     snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
 370     parList+=onePar ;
 371     snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
 372     parList+=onePar ;
 373     snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
 374     parList+=onePar ;
 375     snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
 376     parList+=onePar ;
 377     snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
 378     parList+=onePar ;
 379
 380     if(fPHOSWeightFormula){
 381       snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ;
 382       parList+=onePar;
 383       snprintf(onePar,buffersize,"PHOS Pi0    Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data()    ) ;
 384       parList+=onePar;
 385     }
 386   }
 387   else {
 388     snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f  (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ;
 389     parList+=onePar ;
 390     snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f  (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ;
 391     parList+=onePar ;
 392     snprintf(onePar,buffersize,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
 393     parList+=onePar ;
 394     snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f  (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ;
 395     parList+=onePar ;
 396
 397   }
 398
 399   return parList;
 400
 401 }
 402
 403 //________________________________________________
 404 void AliCaloPID::Print(const Option_t * opt) const
 405 {
 406
 407   //Print some relevant parameters set for the analysis
 408   if(! opt)
 409     return;
 410
 411   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 412
 413   if(fUseBayesianWeights){
 414     printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
 415            fPHOSPhotonWeight,  fPHOSPi0Weight,
 416            fPHOSElectronWeight,  fPHOSChargeWeight,   fPHOSNeutralWeight) ;
 417     printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
 418            fEMCALPhotonWeight,  fEMCALPi0Weight,
 419            fEMCALElectronWeight,  fEMCALChargeWeight,  fEMCALNeutralWeight) ;
 420
 421     printf("PHOS Parametrized weight on?  =     %d\n",  fPHOSWeightFormula) ;
 422     if(fPHOSWeightFormula){
 423       printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
 424       printf("Pi0    weight formula = %s\n", fPHOSPi0WeightFormulaExpression   .Data());
 425     }
 426     if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux?    = %d\n",fParticleFlux);
 427   }
 428   else {
 429     printf("TOF cut        = %e\n",fTOFCut);
 430     printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n",fEMCALL0CutMin, fEMCALL0CutMax);
 431     printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n",fEMCALDEtaCut, fEMCALDPhiCut);
 432     printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut,fPHOSDispersionCut) ;
 433
 434   }
 435
 436   printf(" \n");
 437
 438 }
 439
 440 //___________________________________________________________________________
 441 void AliCaloPID::SetPIDBits(const TString calo, AliVCluster * cluster,
 442                             AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
 443                             AliVEvent* event)
 444 {
 445   //Set Bits for PID selection
 446
 447   //Dispersion/lambdas
 448   //Double_t disp= cluster->GetDispersion()  ;
 449   Double_t l1  = cluster->GetM20() ;
 450   Double_t l0  = cluster->GetM02() ;
 451   Bool_t isDispOK = kTRUE ;
 452   if(cluster->IsPHOS()){
 453     if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
 454     else                                                        isDispOK = kFALSE;
 455   }
 456   else{//EMCAL
 457
 458     if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
 459
 460   }
 461
 462   ph->SetDispBit(isDispOK) ;
 463
 464   //TOF
 465   Double_t tof=cluster->GetTOF()  ;
 466   ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
 467
 468   //Charged
 469   Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
 470
 471   ph->SetChargedBit(isNeutral);
 472
 473   //Set PID pdg
 474   ph->SetIdentifiedParticleType(GetIdentifiedParticleType(calo,*ph->Momentum(),cluster));
 475
 476   if(fDebug > 0){
 477     printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);
 478     printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
 479            ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
 480   }
 481 }
 482
 483 //_________________________________________________________
 484 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
 485                                   AliCalorimeterUtils * cu,
 486                                   AliVEvent* event) const
 487 {
 488   //Check if there is any track attached to this cluster
 489
 490   Int_t nMatches = cluster->GetNTracksMatched();
 491   AliVTrack * track = 0;
 492   Double_t p[3];
 493
 494   if(nMatches > 0){
 495
 496     //In case of ESDs, by default without match one entry with negative index, no match, reject.
 497     if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
 498     {
 499       Int_t iESDtrack = cluster->GetTrackMatchedIndex();
 500       if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
 501       else return kFALSE;
 502
 503       if (!track){
 504         printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n");
 505         return kFALSE;
 506       }
 507     }
 508     else { // AOD
 509       track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
 510       if (!track){
 511         printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n");
 512         return kFALSE;
 513       }
 514     }
 515
 516     Float_t dZ  = cluster->GetTrackDz();
 517     Float_t dR  = cluster->GetTrackDx();
 518
 519     // if track matching was recalculated
 520     if(cluster->IsEMCAL() && cu && cu->IsRecalculationOfClusterTrackMatchingOn()){
 521       dR = 2000., dZ = 2000.;
 522       cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 523     }
 524
 525     if(cluster->IsPHOS()) {
 526
 527       track->GetPxPyPz(p) ;
 528       TLorentzVector trackmom(p[0],p[1],p[2],0);
 529       Int_t charge = track->Charge();
 530       Double_t mf  = event->GetMagneticField();
 531       if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
 532       else                                                                    return kFALSE;
 533
 534     }//PHOS
 535     else {//EMCAL
 536
 537     if(fDebug > 0)
 538         printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut);
 539
 540       if(TMath::Abs(dR) < fEMCALDPhiCut &&
 541          TMath::Abs(dZ) < fEMCALDEtaCut)   return kTRUE;
 542       else                                 return kFALSE;
 543
 544     }//EMCAL cluster
 545
 546
 547   } // more than 1 match, at least one track in array
 548   else return kFALSE;
 549
 550 }
 551
 552 //___________________________________________________________________________________________________
 553 Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const
 554 {
 555   //Check if cluster photon-like. Uses photon cluster parameterization in real pp data
 556   //Returns distance in sigmas. Recommended cut 2.5
 557
 558   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 559   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 560   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 561   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 562   Double_t c       =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 563   Double_t r2      = 0.5*  (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 564                      0.5*  (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 565                      0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 566
 567   if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(r2), fPHOSDispersionCut);
 568
 569   return TMath::Sqrt(r2) ;
 570
 571 }
 572
 573 //_______________________________________________________________________________________________
 574 Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx,  const Double_t dz, const Double_t pt,
 575                                         const Int_t charge, const Double_t mf) const
 576 {
 577   //Checks distance to the closest track. Takes into account
 578   //non-perpendicular incidence of tracks.
 579   //returns distance in sigmas. Recommended cut: 2.
 580   //Requires (sign) of magnetic filed. onc can find it for example as following
 581   //  Double_t mf=0. ;
 582   //  AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 583   //  if(event)
 584   //    mf = event->GetMagneticField(); //Positive for ++ and negative for --
 585
 586
 587   Double_t meanX = 0.;
 588   Double_t meanZ = 0.;
 589   Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 590                            6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
 591                            1.59219);
 592   Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
 593                            1.60) ;
 594
 595   if(mf<0.){ //field --
 596     meanZ = -0.468318 ;
 597     if(charge>0)
 598       meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
 599                          0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
 600     else
 601       meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
 602                          1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
 603   }
 604   else{ //Field ++
 605     meanZ = -0.468318;
 606     if(charge>0)
 607       meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
 608                          0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
 609     else
 610       meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
 611                          0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
 612   }
 613
 614   Double_t rz = (dz-meanZ)/sz ;
 615   Double_t rx = (dx-meanX)/sx ;
 616
 617   if(fDebug > 0)
 618     printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut);
 619
 620   return TMath::Sqrt(rx*rx+rz*rz) ;
 621
 622 }