Updated code to configure Bayeisan PID (Rossella)
authorprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Sep 2012 15:36:39 +0000 (15:36 +0000)
committerprino <prino@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 12 Sep 2012 15:36:39 +0000 (15:36 +0000)
PWGHF/vertexingHF/AliAODPidHF.cxx
PWGHF/vertexingHF/AliAODPidHF.h
PWGHF/vertexingHF/AliRDHFCuts.cxx
PWGHF/vertexingHF/macros/makeTFile4CutsLctopKpi.C

index 01b1c8a..12c6559 100644 (file)
 // Authors: D. Caffarri caffarri@pd.infn.it, A.Dainese andrea.dainese@pd.infn.it, S. Dash dash@to.infn.it, F. Prino prino@to.infn.it, R. Romita r.romita@gsi.de, Y. Wang yifei@pi0.physi.uni-heidelberg.de P. Antonioli pietro.antonioli@bo.infn.it
 //***********************************************************
 #include <TCanvas.h>
+#include <TString.h>
+#include <TH1F.h>
+#include <TF1.h>
+#include <TFile.h>
 
 #include "AliAODPidHF.h"
 #include "AliAODPid.h"
@@ -63,7 +67,10 @@ AliAODPidHF::AliAODPidHF():
   fPtThresholdTPC(999999.),
   fPidResponse(0),
   fPidCombined(new AliPIDCombined()),
-  fTPCResponse(new AliTPCPIDResponse())
+  fTPCResponse(new AliTPCPIDResponse()),
+  fPriorsH(),
+  fCombDetectors(kTPCTOF),
+  fUseCombined(kFALSE)
 {
  //
  // Default constructor
@@ -95,6 +102,9 @@ AliAODPidHF::~AliAODPidHF()
  //   if(fnSigma)  delete fnSigma;
  //   if(fPriors)  delete fPriors;
   delete fTPCResponse;
+  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
+    delete fPriorsH[ispecies];
+  }
 }
 //------------------------
 AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
@@ -126,7 +136,9 @@ AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
   fPtThresholdTPC(pid.fPtThresholdTPC),
   fPidResponse(pid.fPidResponse),
   fPidCombined(pid.fPidCombined),
-  fTPCResponse(0x0)
+  fTPCResponse(0x0),
+  fCombDetectors(pid.fCombDetectors),
+  fUseCombined(pid.fUseCombined)
 {
   
   fnSigma = new Double_t[fnNSigma];
@@ -141,6 +153,13 @@ AliAODPidHF::AliAODPidHF(const AliAODPidHF& pid) :
   for(Int_t i=0;i<fnPLimit;i++){
     fPLimit[i]=pid.fPLimit[i];
   }
+  fPriors = new Double_t[fnPriors];
+  for(Int_t i=0;i<fnPriors;i++){
+    fPriors[i]=pid.fPriors[i];
+  }
+  for(Int_t i=0;i<AliPID::kSPECIES;i++){
+    fPriorsH[i]=pid.fPriorsH[i];
+  }
 
   if(pid.fTPCResponse) fTPCResponse = new AliTPCPIDResponse(*(pid.fTPCResponse));
   //fPidResponse = new AliPIDResponse(*(pid.fPidResponse));
@@ -743,6 +762,7 @@ void AliAODPidHF::SetBetheBloch() {
 }
 //-----------------------
 Bool_t AliAODPidHF::IsTOFPiKexcluded(AliAODTrack *track,Double_t nsigmaK){
+  // TOF proton compatibility
 
   if(!CheckTOFPIDStatus(track)) return 0;
 
@@ -798,7 +818,8 @@ void AliAODPidHF::DrawPrior(AliPID::EParticleType type){
 
 //--------------------------------------------------------------------------
 Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsigma) const{
+  // get n sigma for TPC 
+
   if(!CheckTPCPIDStatus(track)) return -1;
   
   Double_t nsigmaTPC=-999;
@@ -825,6 +846,7 @@ Int_t AliAODPidHF::GetnSigmaTPC(AliAODTrack *track, Int_t species, Double_t &nsi
 //-----------------------------
 
 Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsigma) const{
+  // get n sigma for TOF
 
   if(!CheckTOFPIDStatus(track)) return -1;
 
@@ -839,8 +861,8 @@ Int_t AliAODPidHF::GetnSigmaTOF(AliAODTrack *track,Int_t species, Double_t &nsig
 }
 
 //-----------------------
-Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut,
-                              TString detectors) {
+Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t nsigmaCut, TString detectors) {
+  // Exclude a given hypothesis (labelTracks) in detector
 
   if (detectors.Contains("ITS")) {
 
@@ -875,3 +897,48 @@ Bool_t AliAODPidHF::IsExcluded(AliAODTrack *track, Int_t labelTrack, Double_t ns
 
 }
 //-----------------------------
+void AliAODPidHF::SetPriorsHistos(TString priorFileName){
+  // Set histograms with priors
+  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
+    if(fPriorsH[ispecies]) delete fPriorsH[ispecies];
+    TString nt ="name";
+    nt+="_prior_";
+    nt+=AliPID::ParticleName(ispecies);
+    fPriorsH[ispecies]=new TH1F(nt,nt,100,0,10);
+  }
+  TFile *priorFile=TFile::Open(priorFileName);
+  if (priorFile) {
+    fPriorsH[AliPID::kProton]->Add(static_cast<TH1*>(priorFile->Get("priors3step9")));
+    fPriorsH[AliPID::kKaon  ]->Add(static_cast<TH1*>(priorFile->Get("priors2step9")));
+    fPriorsH[AliPID::kPion  ]->Add(static_cast<TH1*>(priorFile->Get("priors1step9")));
+    delete priorFile;
+    TF1 *salt=new TF1("salt","1.e-10",0,10);
+    fPriorsH[AliPID::kProton]->Add(salt);
+    fPriorsH[AliPID::kKaon  ]->Add(salt);
+    fPriorsH[AliPID::kPion  ]->Add(salt);
+    delete salt;
+  }
+}
+//----------------------------------
+void AliAODPidHF::SetUpCombinedPID(){
+  // Configuration of combined Bayesian PID
+
+ fPidCombined->SetSelectedSpecies(AliPID::kSPECIES);
+  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
+    fPidCombined->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),fPriorsH[ispecies]);
+  }
+ switch (fCombDetectors){
+  case kTPCTOF:
+   fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetTOF);
+  break;
+  case kTPCITS:
+   fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC|AliPIDResponse::kDetITS);
+  break;
+  case kTPC:
+   fPidCombined->SetDetectorMask(AliPIDResponse::kDetTPC);
+  break;
+  case kTOF:
+   fPidCombined->SetDetectorMask(AliPIDResponse::kDetTOF);
+  break;
+ }
+}
index fe43073..2e93f35 100644 (file)
 //// Authors: D. Caffarri caffarri@pd.infn.it, A.Dainese andrea.dainese@pd.infn.it, S. Dash dash@to.infn.it, F. Prino prino@to.infn.it, R. Romita r.romita@gsi.de, Y. Wang yifei@pi0.physi.uni-heidelberg.de
 ////***********************************************************
 
+#include <TString.h>
+#include <TH1F.h>
 #include "AliAODPid.h"
 #include "AliAODTrack.h"
 #include "AliPIDResponse.h"
 #include "AliPIDCombined.h"
+#include "AliPID.h"
 
 class AliAODPidHF : public AliAODPid{
 
  public:
 
+ enum ECombDetectors {
+  kTPC,
+  kTOF,
+  kTPCTOF,
+  kTPCITS
+ };
+
  AliAODPidHF();
  AliAODPidHF(const AliAODPidHF& pid);
  AliAODPidHF& operator=(const AliAODPidHF& pid);
@@ -57,6 +67,9 @@ class AliAODPidHF : public AliAODPid{
  void SetOldPid(Bool_t oldPid){fOldPid=oldPid;return;}
  void SetPtThresholdTPC(Double_t ptThresholdTPC){fPtThresholdTPC=ptThresholdTPC;return;}
  void SetPidResponse(AliPIDResponse *pidResp) {fPidResponse=pidResp;return;}
+ void SetCombDetectors(ECombDetectors pidComb) {
+   fCombDetectors=pidComb;
+  }
  
  //Getters
  
@@ -88,6 +101,10 @@ class AliAODPidHF : public AliAODPid{
  Double_t GetPtThresholdTPC(){return fPtThresholdTPC;}
  AliPIDResponse *GetPidResponse() const {return fPidResponse;}
  AliPIDCombined *GetPidCombined() const {return fPidCombined;}
+ ECombDetectors GetCombDetectors() const {
+   return fCombDetectors;
+  }
+ Bool_t GetUseCombined() {return fUseCombined;}
 
  Int_t RawSignalPID (AliAODTrack *track, TString detector) const;
  Bool_t IsKaonRaw (AliAODTrack *track, TString detector) const;
@@ -117,6 +134,9 @@ class AliAODPidHF : public AliAODPid{
  void SetSelectedSpecies(Int_t ispecies = AliPID::kSPECIES){GetPidCombined()->SetSelectedSpecies(ispecies);};
  void SetPriorDistribution(AliPID::EParticleType type,TH1F *prior);
  void DrawPrior(AliPID::EParticleType type);
+ void SetPriorsHistos(TString priorFileName);
+ void SetUpCombinedPID();
+ void SetUseCombined(Bool_t useCombined=kTRUE) {fUseCombined=useCombined;}
 
  Int_t ApplyPidTPCRaw(AliAODTrack *track,Int_t specie) const;
  Int_t ApplyPidTOFRaw(AliAODTrack *track,Int_t specie) const;
@@ -156,7 +176,11 @@ class AliAODPidHF : public AliAODPid{
 
  AliTPCPIDResponse* fTPCResponse; //! TPC response 
 
- ClassDef(AliAODPidHF,17) // AliAODPid for heavy flavor PID
+ TH1F* fPriorsH[AliPID::kSPECIES]; // priors histos
+ ECombDetectors fCombDetectors; // detectors to be involved for combined PID
+ Bool_t fUseCombined; // detectors to be involved for combined PID
+
+ ClassDef(AliAODPidHF,18) // AliAODPid for heavy flavor PID
 
 };
 
index eaa741b..96496d3 100644 (file)
@@ -305,6 +305,7 @@ Bool_t AliRDHFCuts::IsEventSelected(AliVEvent *event) {
       AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
       fPidHF->SetPidResponse(pidResp);
     }
+      if(fPidHF->GetUseCombined()) fPidHF->SetUpCombinedPID();
     if(fPidHF->GetOldPid()) {
       // pp, from LHC10d onwards
       if((event->GetRunNumber()>121693 && event->GetRunNumber()<136851) ||
index 4fef0cc..e125380 100644 (file)
@@ -5,6 +5,7 @@
 #include <TClonesArray.h>
 #include <TParameter.h>
 #include <TF1.h>
+#include <TCanvas.h>
 
 
 //Use:
 //macro to make a .root file which contains an AliRDHFCutsLctopKpi for AliAnalysisTaskSELambdac task
 
 void SetupCombinedPID(AliRDHFCutsLctopKpi *cutsObj,Double_t threshold,TString priorFileName="noferini-priors.root") {
-  AliPIDCombined *pid=cutsObj->GetPidHF()->GetPidCombined();
-  pid->SetSelectedSpecies(AliPID::kSPECIES);
-  pid->SetDetectorMask(AliPIDResponse::kDetITS
-                      |AliPIDResponse::kDetTPC
-                      |AliPIDResponse::kDetTOF);
-  TH1F *priors[5]; // 5 to make PROOF happy
-  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
-    TString nt ="name";
-    nt+="_prior_";
-    nt+=AliPID::ParticleName(ispecies);
-    priors[ispecies]=new TH1F(nt,nt,100,0,10);
-  }
-  TFile *priorFile=TFile::Open(priorFileName);
-  if (priorFile) {
-    priors[AliPID::kProton]->Add(static_cast<TH1*>(priorFile->Get("priors3step9")));
-    priors[AliPID::kKaon  ]->Add(static_cast<TH1*>(priorFile->Get("priors2step9")));
-    priors[AliPID::kPion  ]->Add(static_cast<TH1*>(priorFile->Get("priors1step9")));
-    delete priorFile;
-    TF1 *salt=new TF1("salt","1.e-10",0,10);
-    priors[AliPID::kProton]->Add(salt);
-    priors[AliPID::kKaon  ]->Add(salt);
-    priors[AliPID::kPion  ]->Add(salt);
-    delete salt;
-  }
-  else {
-    TF1 *flat=new TF1("flat","1",0,10);
-    priors[AliPID::kProton]->Add(flat,1.0); // ... who likes 
-    priors[AliPID::kKaon  ]->Add(flat,1.0); //  these priors
-    priors[AliPID::kPion  ]->Add(flat,1.0); //   anyways?? - Rossella
-    //    priors[AliPID::kProton]->Add(flat,0.162); // from 900 GeV identified particle 
-    //    priors[AliPID::kKaon  ]->Add(flat,0.366); // paper (pp)
-    //    priors[AliPID::kPion  ]->Add(flat,2.977); // dN/dy
-    delete flat;
-  }
-  for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies) {
-    pid->SetPriorDistribution(static_cast<AliPID::EParticleType>(ispecies),priors[ispecies]);
-  }
 
+  cutsObj->GetPidHF()->SetUseCombined(kTRUE);
+  cutsObj->GetPidHF()->SetPriorsHistos(priorFileName);
+  cutsObj->GetPidHF()->SetCombDetectors(AliAODPidHF::kTPCTOF);
   for (Int_t ispecies=0;ispecies<AliPID::kSPECIES;++ispecies)
     cutsObj->SetPIDThreshold(static_cast<AliPID::EParticleType>(ispecies),threshold);
 }
@@ -238,11 +205,17 @@ void makeInputAliAnalysisTaskSELctopKpi(){
   RDHFLctopKpiAn->SetPidprot(pidObjp);
 
   // uncomment these lines for Baysian PID:
-  // Double_t threshold=0.3;
-  // SetupCombinedPID(RDHFLctopKpiAn  ,threshold);
-  // SetupCombinedPID(RDHFLctopKpiProd,threshold);
-  // RDHFLctopKpiAn  ->SetPIDStrategy(AliRDHFCutsLctopKpi::kCombined);
-  // RDHFLctopKpiProd->SetPIDStrategy(AliRDHFCutsLctopKpi::kCombined);
+   Double_t threshold=0.3;
+   SetupCombinedPID(RDHFLctopKpiAn  ,threshold);
+   SetupCombinedPID(RDHFLctopKpiProd,threshold);
+   RDHFLctopKpiAn  ->SetPIDStrategy(AliRDHFCutsLctopKpi::kCombined);
+   RDHFLctopKpiProd->SetPIDStrategy(AliRDHFCutsLctopKpi::kCombined);
+
+
+   AliPIDCombined *pid=RDHFLctopKpiAn->GetPidHF()->GetPidCombined();
+   TH1* hplot=pid->GetPriorDistribution(static_cast<AliPID::EParticleType>(3));
+   TCanvas *c1=new TCanvas();
+   hplot->Draw();
   //
 
 
@@ -265,7 +238,7 @@ void makeInputAliAnalysisTaskSELctopKpi(){
   cout<<"This is the object I'm going to save:"<<endl;
   RDHFLctopKpiProd->PrintAll();
   RDHFLctopKpiAn->PrintAll();
-  TFile* fout=new TFile("cuts4LctopKpi.root","RECREATE"); 
+  TFile* fout=new TFile("prova.root","RECREATE"); 
   fout->cd();
   RDHFLctopKpiProd->Write();
   RDHFLctopKpiAn->Write();