]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWG3/hfe/AliHFEpid.cxx
compilation warning
[u/mrichter/AliRoot.git] / PWG3 / hfe / AliHFEpid.cxx
index 578fd2040ba5c5f7bbf24c440e0881f76835f462..4ba66efccb26a4c6d32020d934fcc30313f18e91 100644 (file)
 #include <TF1.h>
 #include <TIterator.h>
 #include <TList.h>
+#include <TMath.h>
 #include <TObjArray.h>
 #include <TObjString.h>
 #include <TString.h>
 
-#include "AliAODpidUtil.h"
-#include "AliESDpid.h"
 #include "AliLog.h"
 #include "AliPID.h"
 #include "AliVParticle.h"
 
 #include "AliHFEcontainer.h"
 #include "AliHFEpid.h"
-#include "AliHFEpidBase.h"
 #include "AliHFEpidQAmanager.h"
 #include "AliHFEpidITS.h"
 #include "AliHFEpidTPC.h"
 #include "AliHFEpidTRD.h"
 #include "AliHFEpidTOF.h"
+#include "AliHFEpidEMCAL.h"
 #include "AliHFEpidMC.h"
+#include "AliHFEpidBayes.h"
 #include "AliHFEvarManager.h"
 
 ClassImp(AliHFEpid)
 
 const Char_t* AliHFEpid::fgkDetectorName[AliHFEpid::kNdetectorPID + 1] = {
   "MCPID",
-  "ESDPID",
+  "BAYESPID",
   "ITSPID",
   "TPCPID",
   "TRDPID",
@@ -74,7 +74,7 @@ AliHFEpid::AliHFEpid():
   //
   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
   memset(fDetectorOrder, kUndefined, sizeof(UInt_t) * kNdetectorPID);
-  memset(fDetectorOrder, 0, sizeof(UInt_t) * kNdetectorPID);
+  memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
 }
 
 //____________________________________________________________
@@ -94,9 +94,11 @@ AliHFEpid::AliHFEpid(const Char_t *name):
   memset(fSortedOrder, 0, sizeof(UInt_t) * kNdetectorPID);
 
   fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
+  fDetectorPID[kBAYESpid] = new AliHFEpidBayes("BAYESPID");
   fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPCPID");
   fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
   fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
+  fDetectorPID[kEMCALpid] = new AliHFEpidEMCAL("EMCALPID");
 
 }
 
@@ -189,6 +191,7 @@ void AliHFEpid::AddDetector(TString detector, UInt_t position){
   UInt_t detectorID = kUndefined;
   detector.ToUpper();
   if(!detector.CompareTo("MC")) detectorID = kMCpid;
+  else if(!detector.CompareTo("BAYES")) detectorID = kBAYESpid;
   else if(!detector.CompareTo("TPC")) detectorID = kTPCpid;
   else if(!detector.CompareTo("TRD")) detectorID = kTRDpid;
   else if(!detector.CompareTo("TOF")) detectorID = kTOFpid;
@@ -203,27 +206,27 @@ void AliHFEpid::AddDetector(TString detector, UInt_t position){
 }
 
 //____________________________________________________________
-Bool_t AliHFEpid::InitializePID(){
+Bool_t AliHFEpid::InitializePID(Int_t run){
   //
   // Initializes the PID object
   //
 
-  TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
-  // Initialize PID Objects
+  if(!TestBit(kDetectorsSorted)) SortDetectors();
   Bool_t status = kTRUE;
   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
     if(!IsDetectorOn(idet)) continue;
     if(fDetectorPID[idet]){ 
-      status &= fDetectorPID[idet]->InitializePID();
+      status &= fDetectorPID[idet]->InitializePID(run);
       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
     }
   }
+  SetBit(kIsInit);
   PrintStatus();
   return status;
 }
 
 //____________________________________________________________
-Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
+Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
   //
   // Select Tracks
   //
@@ -262,23 +265,38 @@ Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont, cons
 }
 
 //____________________________________________________________
-void AliHFEpid::SetESDpid(AliESDpid *pid){
+void AliHFEpid::SortDetectors(){
+  //
+  // Make sorted list of detectors
+  //
+  if(TestBit(kDetectorsSorted)) return; // Don't sort detectors when they are already sorted
+  TMath::Sort(static_cast<UInt_t>(kNdetectorPID), fDetectorOrder, fSortedOrder, kFALSE);
+  SetBit(kDetectorsSorted);
+}
+
+//____________________________________________________________
+void AliHFEpid::SetPIDResponse(const AliPIDResponse * const pid){
   //
   // Set ESD PID to the Detector PID objects
   //
   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
-    if(fDetectorPID[idet]) fDetectorPID[idet]->SetESDpid(pid);
+    if(fDetectorPID[idet]) fDetectorPID[idet]->SetPIDResponse(pid);
   }
 }
 
 //____________________________________________________________
-void AliHFEpid::SetAODpid(AliAODpidUtil *pid){
+const AliPIDResponse *AliHFEpid::GetPIDResponse() const {
   //
-  // Set ESD PID to the Detector PID objects
+  // Return PID response function
   //
+  const AliPIDResponse *response = NULL;
   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
-    if(fDetectorPID[idet]) fDetectorPID[idet]->SetAODpid(pid);
-  }
+    if(fDetectorPID[idet]){
+      response = fDetectorPID[idet]->GetPIDResponse();
+      break;
+    }
+  } 
+  return response;
 }
 
 //____________________________________________________________
@@ -306,58 +324,110 @@ void AliHFEpid::ConfigureTPCrejectionSimple(){
 }
 
 //____________________________________________________________
-void AliHFEpid::ConfigureTPCrejection(const char *lowerCutParam, Double_t *params){
+void AliHFEpid::ConfigureTOF(Float_t TOFCut){
   //
-  // Combined TPC-TOF PID
+  // Set Number of sigmas for TOF PID
+  //
+  AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
+  if(tofpid) tofpid->SetTOFnSigma(TOFCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+  //
+  // Cofigure centrality dependent cut function for TPC PID 
+  //
+  ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+  //
+  // Cofigure default cut function for TPC PID 
+  //
+  ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+  //
+  // Cofigure cut function for TPC PID 
   // if no function parameterizaion is given, then the default one (exponential) is chosen
   //
+  
   if(HasMCData()) AliInfo("Configuring TPC for MC\n");
   AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
-  AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
-  if(tofpid) tofpid->SetTOFnSigma(3);
-
   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
-  TF1 *upperCut = new TF1("upperCut", "[0]", 0, 20); // Use constant upper cut
-  TF1 *lowerCut = new TF1("lowerCut", lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
-  upperCut->SetParameter(0, 3.);
-  //upperCut->SetParameter(0, 2.7);
-  //upperCut->SetParameter(1, -0.4357);
+  TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
+  TF1 *lowerCut = new TF1(Form("lowerCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
+
+  upperCut->SetParameter(0, upperTPCCut); // pp
+
   if(params){
-    for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++) lowerCut->SetParameter(ipar, params[ipar]);
+      for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
+      {
+         lowerCut->SetParameter(ipar, params[ipar]);
+        //  printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
+      }
   } else {
     // Set default parameterization
     if(HasMCData()) lowerCut->SetParameter(0, -2.5);
-    else lowerCut->SetParameter(0, -3.71769);
-    //else lowerCut->SetParameter(0, -3.7);
-
-    lowerCut->SetParameter(1, -0.40263);
-    //lowerCut->SetParameter(1, -0.8);
-
+    else lowerCut->SetParameter(0, -4.03);  //pp
+    lowerCut->SetParameter(1, -0.22); // pp
+    
     if(HasMCData()) lowerCut->SetParameter(2, -2.2);
-    else lowerCut->SetParameter(2, 0.267857);
-    //else lowerCut->SetParameter(2, -0.35);
+    else lowerCut->SetParameter(2, 0.92); //pp
   }
-  
+
+
   if(tpcpid){
     tpcpid->SetTPCnSigma(2);
-    tpcpid->SetUpperSigmaCut(upperCut);
-    tpcpid->SetLowerSigmaCut(lowerCut);
+    if(centralityBin < 0){
+      tpcpid->SetUpperSigmaCutDefault(upperCut);
+      tpcpid->SetLowerSigmaCutDefault(lowerCut);
+    } else {
+      tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
+      tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
+    }
   }
   AddCommonObject(upperCut);
   AddCommonObject(lowerCut);
 }
 
 //____________________________________________________________
-void AliHFEpid::ConfigureTPCstrategyParis(){
+void AliHFEpid::ConfigureBayesDetectorMask(Int_t detmask){
   //
-  // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
-  //   
-  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
-  if(pid){
-    pid->SetTPCnSigma(2);
-    pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
-    pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
+  // Configure detector mask for Bayes PID
+  // if no detector mask is set the default mask is chosen
+  //
+  
+  if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
+  AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
+
+  if(bayespid)
+  {
+      bayespid->SetBayesDetectorMask(detmask);
+      printf("detector mask in pid class %i \n",detmask);
+  }
+
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureBayesPIDThreshold(Float_t pidthres){
+  //
+  // Configure pid threshold for Bayes PID
+  // if no threshold is set the default threshold is chosen
+  //
+  
+  if(HasMCData()) AliInfo("Configuring Bayes for MC\n");
+  AliHFEpidBayes *bayespid = dynamic_cast<AliHFEpidBayes *>(fDetectorPID[kBAYESpid]);
+
+  if(bayespid)
+  {
+      bayespid->SetBayesPIDThreshold(pidthres);
+      printf("combined pidthreshold %f \n",pidthres);
   }
+
 }
 
 //____________________________________________________________
@@ -369,7 +439,7 @@ void AliHFEpid::PrintStatus() const {
   printf("===============================================\n");
   printf("PID Detectors: \n");
   Int_t npid = 0;
-  TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF"};
+  TString detectors[kNdetectorPID] = {"MC", "BAYES", "ITS", "TPC", "TRD", "TOF", "EMCAL"};
   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
     if(IsDetectorOn(idet)){
       printf("\t%s\n", detectors[idet].Data());