PWG4/PartCorrBase/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 /* $Id: AliCaloPID.cxx 21839 2007-10-29 13:49:42Z gustavo $ */
  16
  17 //_________________________________________________________________________
  18 // Class for PID selection with calorimeters
  19 // The Output of the main method GetIdentifiedParticleType is a PDG number identifying the cluster,
  20 // being kPhoton, kElectron, kPi0 ... as defined in the header file
  21 //   - GetIdentifiedParticleType(const TString calo, const TLorentzVector mom, const AliVCluster * cluster)
  22 //      Assignes a PID tag to the cluster, right now there is the possibility to : use bayesian weights from reco,
  23 //      recalculate them (EMCAL) or use other procedures not used in reco.
  24 //      In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
  25 //      and do SwitchOnBayesianRecalculation().
  26 //      To change the PID parameters from Low to High like the ones by default, use the constructor
  27 //      AliCaloPID(flux)
  28 //      where flux is AliCaloPID::kLow or AliCaloPID::kHigh
  29 //      If it is necessary to change the parameters use the constructor
  30 //      AliCaloPID(AliEMCALPIDUtils *utils) and set the parameters before.
  31
  32 //   - GetGetIdentifiedParticleTypeFromBayesian(const TString calo, const Double_t * pid, const Float_t energy)
  33 //      Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle,
  34 //      executed when bayesian is ON by GetIdentifiedParticleType(const TString calo, const TLorentzVector mom, const AliVCluster * cluster)
  35 //   - SetPIDBits: Simple PID, depending on the thresholds fLOCut fTOFCut and even the
  36 //     result of the PID bayesian a different PID bit is set.
  37 //
  38 //   - IsTrackMatched(): Independent method that needs to be combined with GetIdentifiedParticleType to know if the cluster was matched
  39 //
  40 //*-- Author: Gustavo Conesa (INFN-LNF)
  41 //////////////////////////////////////////////////////////////////////////////
  42
  43
  44 // --- ROOT system ---
  45 #include <TMath.h>
  46 #include <TString.h>
  47 #include <TList.h>
  48
  49 // ---- ANALYSIS system ----
  50 #include "AliCaloPID.h"
  51 #include "AliVCluster.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 // Histogram settings
  81 fHistoNEBins(100),        fHistoEMax(100.),            fHistoEMin(0.),
  82 fHistoNDEtaBins(100),     fHistoDEtaMax(0.15),         fHistoDEtaMin(-0.15),
  83 fHistoNDPhiBins(100),     fHistoDPhiMax(0.15),         fHistoDPhiMin(-0.15),
  84 fhTrackMatchedDEta(0x0),  fhTrackMatchedDPhi(0x0),     fhTrackMatchedDEtaDPhi(0x0)
  85 {
  86   //Ctor
  87
  88   //Initialize parameters
  89   InitParameters();
  90 }
  91
  92 //________________________________________
  93 AliCaloPID::AliCaloPID(const Int_t flux) :
  94 TObject(),                fDebug(-1),                  fParticleFlux(flux),
  95 //Bayesian
  96 fEMCALPIDUtils(),         fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
  97 fEMCALPhotonWeight(0.),   fEMCALPi0Weight(0.),
  98 fEMCALElectronWeight(0.), fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
  99 fPHOSPhotonWeight(0.),    fPHOSPi0Weight(0.),
 100 fPHOSElectronWeight(0.),  fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
 101 fPHOSWeightFormula(0),    fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
 102 fPHOSPhotonWeightFormulaExpression(""),
 103 fPHOSPi0WeightFormulaExpression(""),
 104 //PID calculation
 105 fEMCALL0CutMax(100.),     fEMCALL0CutMin(0),
 106 fEMCALDEtaCut(2000.),     fEMCALDPhiCut(2000.),
 107 fTOFCut(0.),
 108 fPHOSDispersionCut(1000), fPHOSRCut(1000),
 109 // Histogram settings
 110 fHistoNEBins(100),        fHistoEMax(100.),            fHistoEMin(0.),
 111 fHistoNDEtaBins(100),     fHistoDEtaMax(0.15),         fHistoDEtaMin(-0.15),
 112 fHistoNDPhiBins(100),     fHistoDPhiMax(0.15),         fHistoDPhiMin(-0.15),
 113 fhTrackMatchedDEta(0x0),  fhTrackMatchedDPhi(0x0),     fhTrackMatchedDEtaDPhi(0x0)
 114 {
 115   //Ctor
 116
 117   //Initialize parameters
 118   InitParameters();
 119
 120 }
 121
 122 //_______________________________________________
 123 AliCaloPID::AliCaloPID(const TNamed * emcalpid) :
 124 TObject(),                   fDebug(-1),                  fParticleFlux(kLow),
 125 //Bayesian
 126 fEMCALPIDUtils((AliEMCALPIDUtils*)emcalpid),
 127 fUseBayesianWeights(kFALSE), fRecalculateBayesian(kFALSE),
 128 fEMCALPhotonWeight(0.),      fEMCALPi0Weight(0.),
 129 fEMCALElectronWeight(0.),    fEMCALChargeWeight(0.),      fEMCALNeutralWeight(0.),
 130 fPHOSPhotonWeight(0.),       fPHOSPi0Weight(0.),
 131 fPHOSElectronWeight(0.),     fPHOSChargeWeight(0.) ,      fPHOSNeutralWeight(0.),
 132 fPHOSWeightFormula(0),       fPHOSPhotonWeightFormula(0), fPHOSPi0WeightFormula(0),
 133 fPHOSPhotonWeightFormulaExpression(""),
 134 fPHOSPi0WeightFormulaExpression(""),
 135 //PID calculation
 136 fEMCALL0CutMax(100.),        fEMCALL0CutMin(0),
 137 fEMCALDEtaCut(2000.),        fEMCALDPhiCut(2000.),
 138 fTOFCut(0.),
 139 fPHOSDispersionCut(1000),    fPHOSRCut(1000),
 140 // Histogram settings
 141 fHistoNEBins(100),           fHistoEMax(100.),            fHistoEMin(0.),
 142 fHistoNDEtaBins(100),        fHistoDEtaMax(0.15),         fHistoDEtaMin(-0.15),
 143 fHistoNDPhiBins(100),        fHistoDPhiMax(0.15),         fHistoDPhiMin(-0.15),
 144 fhTrackMatchedDEta(0x0),     fhTrackMatchedDPhi(0x0),     fhTrackMatchedDEtaDPhi(0x0)
 145 {
 146   //Ctor
 147
 148   //Initialize parameters
 149   InitParameters();
 150 }
 151
 152 //_______________________
 153 AliCaloPID::~AliCaloPID()
 154 {
 155   //Dtor
 156
 157   delete fPHOSPhotonWeightFormula ;
 158   delete fPHOSPi0WeightFormula ;
 159   delete fEMCALPIDUtils ;
 160
 161 }
 162
 163 //___________________________________________
 164 TList *  AliCaloPID::GetCreateOutputObjects()
 165 {
 166   // Create histograms to be saved in output file and
 167   // store them in outputContainer of the analysis class that calls this class.
 168
 169   TList * outputContainer = new TList() ;
 170   outputContainer->SetName("CaloPIDHistos") ;
 171
 172   outputContainer->SetOwner(kFALSE);
 173
 174   fhTrackMatchedDEta  = new TH2F
 175   ("TrackMatchedDEta",
 176    "d#eta of cluster-track vs cluster energy",
 177    fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNDEtaBins,fHistoDEtaMin,fHistoDEtaMax);
 178   fhTrackMatchedDEta->SetYTitle("d#eta");
 179   fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
 180
 181   fhTrackMatchedDPhi  = new TH2F
 182   ("TrackMatchedDPhi",
 183    "d#phi of cluster-track vs cluster energy"
 184    ,fHistoNEBins,fHistoEMin,fHistoEMax,fHistoNDPhiBins,fHistoDPhiMin,fHistoDPhiMax);
 185   fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
 186   fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
 187
 188   fhTrackMatchedDEtaDPhi  = new TH2F
 189   ("TrackMatchedDEtaDPhi",
 190    "d#eta vs d#phi of cluster-track vs cluster energy"
 191    ,fHistoNDEtaBins,fHistoDEtaMin,fHistoDEtaMax,fHistoNDPhiBins,fHistoDPhiMin,fHistoDPhiMax);
 192   fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
 193   fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
 194
 195   outputContainer->Add(fhTrackMatchedDEta) ;
 196   outputContainer->Add(fhTrackMatchedDPhi) ;
 197   outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
 198
 199   return outputContainer;
 200 }
 201
 202
 203
 204 //_______________________________
 205 void AliCaloPID::InitParameters()
 206 {
 207   //Initialize the parameters of the PID.
 208
 209   // Bayesian
 210   fEMCALPhotonWeight   = 0.6 ;
 211   fEMCALPi0Weight      = 0.6 ;
 212   fEMCALElectronWeight = 0.6 ;
 213   fEMCALChargeWeight   = 0.6 ;
 214   fEMCALNeutralWeight  = 0.6 ;
 215
 216   fPHOSPhotonWeight    = 0.6 ;
 217   fPHOSPi0Weight       = 0.6 ;
 218   fPHOSElectronWeight  = 0.6 ;
 219   fPHOSChargeWeight    = 0.6 ;
 220   fPHOSNeutralWeight   = 0.6 ;
 221
 222   //Formula to set the PID weight threshold for photon or pi0
 223   fPHOSWeightFormula       = kFALSE;
 224   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))";
 225   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))"   ;
 226
 227   if(fRecalculateBayesian){
 228     if(fParticleFlux == kLow){
 229       printf("AliCaloPID::Init() - SetLOWFluxParam\n");
 230       fEMCALPIDUtils->SetLowFluxParam() ;
 231     }
 232     else if (fParticleFlux == kHigh){
 233       printf("AliCaloPID::Init() - SetHIGHFluxParam\n");
 234       fEMCALPIDUtils->SetHighFluxParam() ;
 235     }
 236   }
 237
 238   //PID recalculation, not bayesian
 239
 240   //EMCAL
 241   fEMCALL0CutMax = 0.3 ;
 242   fEMCALL0CutMin = 0.01;
 243
 244   fEMCALDPhiCut  = 0.05; // Same cut as in AliEMCALRecoUtils
 245   fEMCALDEtaCut  = 0.025;// Same cut as in AliEMCALRecoUtils
 246
 247   // PHOS / EMCAL, not used
 248   fTOFCut        = 1.e-6;
 249
 250   //PHOS
 251   fPHOSRCut          = 2. ;
 252   fPHOSDispersionCut = 2.5;
 253
 254 }
 255
 256 //______________________________________________
 257 AliEMCALPIDUtils *AliCaloPID::GetEMCALPIDUtils()
 258 {
 259   // return pointer to AliEMCALPIDUtils, create it if needed
 260
 261   if(!fEMCALPIDUtils) fEMCALPIDUtils = new AliEMCALPIDUtils ;
 262   return fEMCALPIDUtils ;
 263
 264 }
 265
 266
 267 //______________________________________________________________________
 268 Int_t AliCaloPID::GetIdentifiedParticleType(const TString calo,
 269                                             const TLorentzVector mom,
 270                                             const AliVCluster * cluster)
 271 {
 272   // Returns a PDG number corresponding to the likely ID of the cluster
 273
 274   Float_t energy  = mom.E();
 275   Float_t lambda0 = cluster->GetM02();
 276   Float_t lambda1 = cluster->GetM20();
 277
 278   // ---------------------
 279   // Use bayesian approach
 280   // ---------------------
 281
 282   if(fUseBayesianWeights){
 283
 284     Double_t weights[AliPID::kSPECIESN];
 285
 286     if(calo == "EMCAL"&& fRecalculateBayesian){
 287       fEMCALPIDUtils->ComputePID(energy, lambda0);
 288       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
 289     }
 290     else {
 291       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i];
 292     }
 293
 294     if(fDebug > 0)  {
 295       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",
 296              calo.Data(),
 297              weights[AliVCluster::kPhoton],    weights[AliVCluster::kPi0],
 298              weights[AliVCluster::kElectron],  weights[AliVCluster::kEleCon],
 299              weights[AliVCluster::kPion],      weights[AliVCluster::kKaon],
 300              weights[AliVCluster::kProton],
 301              weights[AliVCluster::kNeutron],   weights[AliVCluster::kKaon0]);
 302     }
 303
 304     return GetIdentifiedParticleTypeFromBayesWeights(calo, weights, energy);
 305   }
 306
 307   // -------------------------------------------------------
 308   // Calculate PID SS from data, do not use bayesian weights
 309   // -------------------------------------------------------
 310
 311   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",
 312                         calo.Data(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(),
 313                         cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
 314
 315   if(cluster->IsEMCAL()){
 316
 317     if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType() - EMCAL SS %f <%f < %f?\n",fEMCALL0CutMin, lambda0, fEMCALL0CutMax);
 318
 319     if(lambda0 < fEMCALL0CutMax && lambda0 > fEMCALL0CutMin) return kPhoton ;
 320     else                                                     return kNeutralUnknown ;
 321   }//EMCAL
 322   else {//PHOS
 323     if(TestPHOSDispersion(mom.Pt(),lambda0,lambda1) < fPHOSDispersionCut) return kPhoton;
 324     else                                                                  return kNeutralUnknown;
 325   }
 326
 327 }
 328
 329 //_______________________________________________________________________________
 330 Int_t AliCaloPID::GetIdentifiedParticleTypeFromBayesWeights(const TString calo,
 331                                                             const Double_t * pid,
 332                                                             const Float_t energy)
 333 {
 334   //Return most probable identity of the particle after bayesian weights calculated in reconstruction
 335
 336   if(!pid){
 337     printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n");
 338     abort();
 339   }
 340
 341   Float_t wPh  =  fPHOSPhotonWeight ;
 342   Float_t wPi0 =  fPHOSPi0Weight ;
 343   Float_t wE   =  fPHOSElectronWeight ;
 344   Float_t wCh  =  fPHOSChargeWeight ;
 345   Float_t wNe  =  fPHOSNeutralWeight ;
 346
 347   if(calo == "PHOS" && fPHOSWeightFormula){
 348     wPh  = GetPHOSPhotonWeightFormula()->Eval(energy) ;
 349     wPi0 = GetPHOSPi0WeightFormula()   ->Eval(energy);
 350   }
 351
 352   if(calo == "EMCAL"){
 353
 354     wPh  =  fEMCALPhotonWeight ;
 355     wPi0 =  fEMCALPi0Weight ;
 356     wE   =  fEMCALElectronWeight ;
 357     wCh  =  fEMCALChargeWeight ;
 358     wNe  =  fEMCALNeutralWeight ;
 359
 360   }
 361
 362   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",
 363                          calo.Data(),pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
 364                          pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
 365                          pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton],
 366                          pid[AliVCluster::kNeutron], pid[AliVCluster::kKaon0]);
 367
 368   Int_t pdg = kNeutralUnknown ;
 369   Float_t chargedHadronWeight = pid[AliVCluster::kProton]+pid[AliVCluster::kKaon]+
 370   pid[AliVCluster::kPion]+pid[AliVCluster::kMuon];
 371   Float_t neutralHadronWeight = pid[AliVCluster::kNeutron]+pid[AliVCluster::kKaon0];
 372   Float_t allChargedWeight    = pid[AliVCluster::kElectron]+pid[AliVCluster::kEleCon]+ chargedHadronWeight;
 373   Float_t allNeutralWeight    = pid[AliVCluster::kPhoton]+pid[AliVCluster::kPi0]+ neutralHadronWeight;
 374
 375   //Select most probable ID
 376   if(calo=="PHOS"){
 377     if(pid[AliVCluster::kPhoton] > wPh)        pdg = kPhoton ;
 378     else if(pid[AliVCluster::kPi0] > wPi0)     pdg = kPi0 ;
 379     else if(pid[AliVCluster::kElectron] > wE)  pdg = kElectron ;
 380     else if(pid[AliVCluster::kEleCon] >  wE)   pdg = kEleCon ;
 381     else if(chargedHadronWeight > wCh)         pdg = kChargedHadron ;
 382     else if(neutralHadronWeight > wNe)         pdg = kNeutralHadron ;
 383     else if(allChargedWeight >  allNeutralWeight)
 384       pdg = kChargedUnknown ;
 385     else
 386       pdg = kNeutralUnknown ;
 387   }
 388   else{//EMCAL
 389
 390     if(pid[AliVCluster::kPhoton]  > wPh)                     pdg = kPhoton ;
 391     else if(pid[AliVCluster::kElectron]  > wE)               pdg = kElectron ;
 392     else if(pid[AliVCluster::kPhoton]+pid[AliVCluster::kElectron]  > wPh) pdg = kPhoton ; //temporal sollution until track matching for electrons is considered
 393     else if(pid[AliVCluster::kPi0] > wPi0)                   pdg = kPi0 ;
 394     else if(chargedHadronWeight + neutralHadronWeight > wCh) pdg = kChargedHadron ;
 395     else if(neutralHadronWeight + chargedHadronWeight > wNe) pdg = kNeutralHadron ;
 396     else                                                     pdg = kNeutralUnknown ;
 397   }
 398
 399   if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
 400
 401   return pdg ;
 402
 403 }
 404
 405 //_________________________________________
 406 TString  AliCaloPID::GetPIDParametersList()
 407 {
 408   //Put data member values in string to keep in output container
 409
 410   TString parList ; //this will be list of parameters used for this analysis.
 411   const Int_t buffersize = 255;
 412   char onePar[buffersize] ;
 413   snprintf(onePar,buffersize,"--- AliCaloPID ---\n") ;
 414   parList+=onePar ;
 415   if(fUseBayesianWeights){
 416     snprintf(onePar,buffersize,"fEMCALPhotonWeight =%2.2f (EMCAL bayesian weight for photons)\n",fEMCALPhotonWeight) ;
 417     parList+=onePar ;
 418     snprintf(onePar,buffersize,"fEMCALPi0Weight =%2.2f (EMCAL bayesian weight for pi0)\n",fEMCALPi0Weight) ;
 419     parList+=onePar ;
 420     snprintf(onePar,buffersize,"fEMCALElectronWeight =%2.2f(EMCAL bayesian weight for electrons)\n",fEMCALElectronWeight) ;
 421     parList+=onePar ;
 422     snprintf(onePar,buffersize,"fEMCALChargeWeight =%2.2f (EMCAL bayesian weight for charged hadrons)\n",fEMCALChargeWeight) ;
 423     parList+=onePar ;
 424     snprintf(onePar,buffersize,"fEMCALNeutralWeight =%2.2f (EMCAL bayesian weight for neutral hadrons)\n",fEMCALNeutralWeight) ;
 425     parList+=onePar ;
 426     snprintf(onePar,buffersize,"fPHOSPhotonWeight =%2.2f (PHOS bayesian weight for photons)\n",fPHOSPhotonWeight) ;
 427     parList+=onePar ;
 428     snprintf(onePar,buffersize,"fPHOSPi0Weight =%2.2f (PHOS bayesian weight for pi0)\n",fPHOSPi0Weight) ;
 429     parList+=onePar ;
 430     snprintf(onePar,buffersize,"fPHOSElectronWeight =%2.2f(PHOS bayesian weight for electrons)\n",fPHOSElectronWeight) ;
 431     parList+=onePar ;
 432     snprintf(onePar,buffersize,"fPHOSChargeWeight =%2.2f (PHOS bayesian weight for charged hadrons)\n",fPHOSChargeWeight) ;
 433     parList+=onePar ;
 434     snprintf(onePar,buffersize,"fPHOSNeutralWeight =%2.2f (PHOS bayesian weight for neutral hadrons)\n",fPHOSNeutralWeight) ;
 435     parList+=onePar ;
 436
 437     if(fPHOSWeightFormula){
 438       snprintf(onePar,buffersize,"PHOS Photon Weight Formula: %s\n",fPHOSPhotonWeightFormulaExpression.Data() ) ;
 439       parList+=onePar;
 440       snprintf(onePar,buffersize,"PHOS Pi0    Weight Formula: %s\n",fPHOSPi0WeightFormulaExpression.Data()    ) ;
 441       parList+=onePar;
 442     }
 443   }
 444   else {
 445     snprintf(onePar,buffersize,"EMCAL: fEMCALL0CutMin =%2.2f, fEMCALL0CutMax =%2.2f  (Cut on Shower Shape) \n",fEMCALL0CutMin, fEMCALL0CutMax) ;
 446     parList+=onePar ;
 447     snprintf(onePar,buffersize,"EMCAL: fEMCALDEtaCut =%2.2f, fEMCALDPhiCut =%2.2f  (Cut on track matching) \n",fEMCALDEtaCut, fEMCALDPhiCut) ;
 448     parList+=onePar ;
 449     snprintf(onePar,buffersize,"fTOFCut  =%e (Cut on TOF, used in PID evaluation) \n",fTOFCut) ;
 450     parList+=onePar ;
 451     snprintf(onePar,buffersize,"fPHOSRCut =%2.2f, fPHOSDispersionCut =%2.2f  (Cut on Shower Shape and CPV) \n",fPHOSRCut,fPHOSDispersionCut) ;
 452     parList+=onePar ;
 453
 454   }
 455
 456   return parList;
 457
 458 }
 459
 460 //________________________________________________
 461 void AliCaloPID::Print(const Option_t * opt) const
 462 {
 463
 464   //Print some relevant parameters set for the analysis
 465   if(! opt)
 466     return;
 467
 468   printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ;
 469
 470   if(fUseBayesianWeights){
 471     printf("PHOS PID weight , photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f \n",
 472            fPHOSPhotonWeight,  fPHOSPi0Weight,
 473            fPHOSElectronWeight,  fPHOSChargeWeight,   fPHOSNeutralWeight) ;
 474     printf("EMCAL PID weight, photon %0.2f, pi0 %0.2f, e %0.2f, charge %0.2f, neutral %0.2f\n",
 475            fEMCALPhotonWeight,  fEMCALPi0Weight,
 476            fEMCALElectronWeight,  fEMCALChargeWeight,  fEMCALNeutralWeight) ;
 477
 478     printf("PHOS Parametrized weight on?  =     %d\n",  fPHOSWeightFormula) ;
 479     if(fPHOSWeightFormula){
 480       printf("Photon weight formula = %s\n", fPHOSPhotonWeightFormulaExpression.Data());
 481       printf("Pi0    weight formula = %s\n", fPHOSPi0WeightFormulaExpression   .Data());
 482     }
 483     if(fRecalculateBayesian) printf(" Recalculate bayesian with Particle Flux?    = %d\n",fParticleFlux);
 484   }
 485   else {
 486     printf("TOF cut        = %e\n",fTOFCut);
 487     printf("EMCAL Lambda0 cut min = %2.2f; max = %2.2f\n",fEMCALL0CutMin, fEMCALL0CutMax);
 488     printf("EMCAL cluster-track dEta < %2.3f; dPhi < %2.3f\n",fEMCALDEtaCut, fEMCALDPhiCut);
 489     printf("PHOS Treac matching cut =%2.2f, Dispersion Cut =%2.2f \n",fPHOSRCut,fPHOSDispersionCut) ;
 490
 491   }
 492
 493   printf(" \n");
 494
 495 }
 496
 497 //___________________________________________________________________________
 498 void AliCaloPID::SetPIDBits(const TString calo, AliVCluster * cluster,
 499                             AliAODPWG4Particle * ph, AliCalorimeterUtils* cu,
 500                             AliVEvent* event)
 501 {
 502   //Set Bits for PID selection
 503
 504   //Dispersion/lambdas
 505   //Double_t disp= cluster->GetDispersion()  ;
 506   Double_t l1  = cluster->GetM20() ;
 507   Double_t l0  = cluster->GetM02() ;
 508   Bool_t isDispOK = kTRUE ;
 509   if(cluster->IsPHOS()){
 510     if(TestPHOSDispersion(ph->Pt(),l0,l1) < fPHOSDispersionCut) isDispOK = kTRUE;
 511     else                                                        isDispOK = kFALSE;
 512   }
 513   else{//EMCAL
 514
 515     if(l0 > fEMCALL0CutMin && l0 < fEMCALL0CutMax) isDispOK = kTRUE;
 516
 517   }
 518
 519   ph->SetDispBit(isDispOK) ;
 520
 521   //TOF
 522   Double_t tof=cluster->GetTOF()  ;
 523   ph->SetTOFBit(TMath::Abs(tof)<fTOFCut) ;
 524
 525   //Charged
 526   Bool_t isNeutral = IsTrackMatched(cluster,cu,event);
 527
 528   ph->SetChargedBit(isNeutral);
 529
 530   //Set PID pdg
 531   ph->SetIdentifiedParticleType(GetIdentifiedParticleType(calo,*ph->Momentum(),cluster));
 532
 533   if(fDebug > 0){
 534     printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);
 535     printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
 536            ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit());
 537   }
 538 }
 539
 540 //_____________________________________________________________________
 541 Bool_t AliCaloPID::IsTrackMatched(AliVCluster* cluster,
 542                                   AliCalorimeterUtils * cu,
 543                                   AliVEvent* event) const
 544 {
 545   //Check if there is any track attached to this cluster
 546
 547   Int_t nMatches = cluster->GetNTracksMatched();
 548   AliVTrack * track = 0;
 549   Double_t p[3];
 550
 551   if(nMatches > 0){
 552
 553     //In case of ESDs, by default without match one entry with negative index, no match, reject.
 554     if(!strcmp("AliESDCaloCluster",Form("%s",cluster->ClassName())))
 555     {
 556       Int_t iESDtrack = cluster->GetTrackMatchedIndex();
 557       if(iESDtrack >= 0) track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
 558       else return kFALSE;
 559
 560       if (!track){
 561         printf("AliCaloPID::IsTrackMatched() - Null matched track in ESD when index is OK!\n");
 562         return kFALSE;
 563       }
 564     }
 565     else { // AOD
 566       track = dynamic_cast<AliVTrack*> (cluster->GetTrackMatched(0));
 567       if (!track){
 568         printf("AliCaloPID::IsTrackMatched() - Null matched track in AOD!\n");
 569         return kFALSE;
 570       }
 571     }
 572
 573     Float_t clE = cluster->E();
 574     Float_t dZ  = cluster->GetTrackDz();
 575     Float_t dR  = cluster->GetTrackDx();
 576
 577     // if track matching was recalculated
 578     if(cluster->IsEMCAL() &&cu && cu->IsRecalculationOfClusterTrackMatchingOn()){
 579       dR = 2000., dZ = 2000.;
 580       cu->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dR,dZ);
 581     }
 582
 583     // Fill control histograms
 584     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999){
 585       fhTrackMatchedDEta->Fill(clE,dZ);
 586       fhTrackMatchedDPhi->Fill(clE,dR);
 587       if(clE > 0.5)fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
 588       //printf("AliCaloPID::IsTrackMatched - %d dR %f , dZ %f \n",cluster->IsEMCAL(),dR, dZ);
 589     }
 590
 591     if(cluster->IsPHOS()) {
 592
 593       track->GetPxPyPz(p) ;
 594       TLorentzVector trackmom(p[0],p[1],p[2],0);
 595       Int_t charge = track->Charge();
 596       Double_t mf  = event->GetMagneticField();
 597       if(TestPHOSChargedVeto(dR, dZ, trackmom.Pt(), charge, mf ) < fPHOSRCut) return kTRUE;
 598       else                                                                    return kFALSE;
 599
 600     }//PHOS
 601     else {//EMCAL
 602
 603       if(fDebug > 0)
 604         printf("AliCaloPID::IsTrackMatched - EMCAL dR %f < %f, dZ %f < %f \n",dR, fEMCALDPhiCut, dZ, fEMCALDEtaCut);
 605
 606       if(TMath::Abs(dR) < fEMCALDPhiCut &&
 607          TMath::Abs(dZ) < fEMCALDEtaCut)   return kTRUE;
 608       else                                 return kFALSE;
 609
 610     }//EMCAL cluster
 611
 612
 613   } // more than 1 match, at least one track in array
 614   else return kFALSE;
 615
 616 }
 617
 618 //___________________________________________________________________________________________________
 619 Float_t AliCaloPID::TestPHOSDispersion(const Double_t pt, const Double_t l1, const Double_t l2) const
 620 {
 621   //Check if cluster photon-like. Uses photon cluster parameterization in real pp data
 622   //Returns distance in sigmas. Recommended cut 2.5
 623
 624   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
 625   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
 626   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
 627   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
 628   Double_t c       =-0.35-0.550*TMath::Exp(-0.390730*pt) ;
 629   Double_t R2      = 0.5*  (l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
 630                      0.5*  (l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
 631                      0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
 632
 633   if(fDebug > 0) printf("AliCaloPID::TestPHOSDispersion() - PHOS SS R %f < %f?\n", TMath::Sqrt(R2), fPHOSDispersionCut);
 634
 635   return TMath::Sqrt(R2) ;
 636
 637
 638 }
 639
 640 //_______________________________________________________________________________________________
 641 Float_t AliCaloPID::TestPHOSChargedVeto(const Double_t dx,  const Double_t dz, const Double_t pt,
 642                                         const Int_t charge, const Double_t mf) const
 643 {
 644   //Checks distance to the closest track. Takes into account
 645   //non-perpendicular incidence of tracks.
 646   //returns distance in sigmas. Recommended cut: 2.
 647   //Requires (sign) of magnetic filed. onc can find it for example as following
 648   //  Double_t mf=0. ;
 649   //  AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 650   //  if(event)
 651   //    mf = event->GetMagneticField(); //Positive for ++ and negative for --
 652
 653
 654   Double_t meanX = 0.;
 655   Double_t meanZ = 0.;
 656   Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
 657                            6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+
 658                            1.59219);
 659   Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+
 660                            1.60) ;
 661
 662   if(mf<0.){ //field --
 663     meanZ = -0.468318 ;
 664     if(charge>0)
 665       meanX = TMath::Min(7.3, 3.89994*1.20679 *1.20679 /(pt*pt+1.20679*1.20679)+
 666                          0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
 667     else
 668       meanX =-TMath::Min(7.7, 3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+
 669                          1.23114 +4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
 670   }
 671   else{ //Field ++
 672     meanZ = -0.468318;
 673     if(charge>0)
 674       meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+
 675                          0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
 676     else
 677       meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-
 678                          0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
 679   }
 680
 681   Double_t rz = (dz-meanZ)/sz ;
 682   Double_t rx = (dx-meanX)/sx ;
 683
 684   if(fDebug > 0)
 685     printf("AliCaloPID::TestPHOSDispersion() - PHOS Matching R %f < %f\n",TMath::Sqrt(rx*rx+rz*rz), fPHOSRCut);
 686
 687   return TMath::Sqrt(rx*rx+rz*rz) ;
 688
 689 }