Fixes for #89479: change in AliPIDResponse and AliPIDCombined (STEERBase). F.Noferini

diff --git a/STEER/STEERBase/AliPIDCombined.cxx b/STEER/STEERBase/AliPIDCombined.cxx
index 31db1f52be0759622c51a784fb9a426eaec57830..8a66c0428bc0e8dc0666013097adf93f48396efe 100644 (file)
--- a/STEER/STEERBase/AliPIDCombined.cxx
+++ b/STEER/STEERBase/AliPIDCombined.cxx
@@ -32,12 +32,23 @@
 
 #include "AliPIDCombined.h"
 
+#include "TMath.h"
+#include "TFile.h"
+
+#include "AliOADBContainer.h"
+
+// initialize static members
+TH2F* AliPIDCombined::fDefaultPriorsTPC[5];
+
+ClassImp(AliPIDCombined);
+
 AliPIDCombined::AliPIDCombined():
        TNamed("CombinedPID","CombinedPID"),
        fDetectorMask(0),
        fRejectMismatchMask(0x7F),
        fEnablePriors(kTRUE),
-       fSelectedSpecies(AliPID::kSPECIES)
+       fSelectedSpecies(AliPID::kSPECIES),
+       fUseDefaultTPCPriors(kFALSE)
 {      
   //
   // default constructor
@@ -52,7 +63,8 @@ AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
        fDetectorMask(0),
        fRejectMismatchMask(0x7F),
        fEnablePriors(kTRUE),
-       fSelectedSpecies(AliPID::kSPECIES)
+       fSelectedSpecies(AliPID::kSPECIES),
+       fUseDefaultTPCPriors(kFALSE)
 {
   //
   // default constructor with name and title
@@ -162,7 +174,26 @@ UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPID
          p[i] = pITS[i]*pTPC[i]*pTRD[i]*pTOF[i]*pHMPID[i]*pEMCAL[i]*pPHOS[i];
        }
        Double_t priors[fSelectedSpecies];
-       if (fEnablePriors) GetPriors(track,priors);
+       if (fEnablePriors){
+         GetPriors(track,priors,response->GetCurrentCentrality());
+         
+         // apply tof matching efficiency to priors if TOF joined PID for this track
+         if(fUseDefaultTPCPriors && (usedMask & AliPIDResponse::kDetTOF)){
+           Double_t pt=TMath::Abs(track->Pt());
+           Float_t kaonTOFfactor = 0.1;
+           if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
+           Float_t protonTOFfactor = 0.1;
+           if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
+           
+           if(track->Charge() < 0){ // for negative tracks
+             kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
+             protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
+           }
+           
+           p[3] *= kaonTOFfactor;
+           p[4] *= protonTOFfactor;
+         }
+       }
        else { for (Int_t i=0;i<fSelectedSpecies;i++) priors[i]=1.;}
        ComputeBayesProbabilities(bayesProbabilities,p,priors);
        return usedMask;
@@ -170,13 +201,26 @@ UInt_t AliPIDCombined::ComputeProbabilities(const AliVTrack *track, const AliPID
 
 
 //-------------------------------------------------------------------------------------------------
-void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p) const {
+void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p,Float_t centrality) const {
        
        //
        // get priors from distributions
        //
        
        Double_t pt=TMath::Abs(track->Pt());
+
+        if(fUseDefaultTPCPriors){ // use default priors if requested
+         Float_t usedCentr = centrality+5*(centrality>0);
+         if(usedCentr < -0.99) usedCentr = -0.99;
+         else if(usedCentr > 99) usedCentr = 99;
+         if(pt > 9.99) pt = 9.99;
+         else if(pt < 0.01)  pt = 0.01;
+         
+         for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
+
+         return;
+       }
+       
        Double_t sumPriors = 0;
        for (Int_t i=0;i<fSelectedSpecies;++i){
                if (fPriorsDistributions[i]){
@@ -195,7 +239,61 @@ void AliPIDCombined::GetPriors(const AliVTrack *track, Double_t* p) const {
        }
        return;
 }
+//-------------------------------------------------------------------------------------------------
+void AliPIDCombined::GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const{
+       
+       //
+       // get priors from distributions
+       //
+        Double_t pt=TMath::Abs(track->Pt());
+       
+        if(fUseDefaultTPCPriors){ // use default priors if requested
+         Float_t usedCentr = response->GetCurrentCentrality()+5*(response->GetCurrentCentrality()>0);
+         if(usedCentr < -0.99) usedCentr = -0.99;
+         else if(usedCentr > 99) usedCentr = 99;
+         if(pt > 9.99) pt = 9.99;
+         else if(pt < 0.01)  pt = 0.01;
+         
+         for(Int_t i=0;i<5;i++) p[i] = fDefaultPriorsTPC[i]->Interpolate(usedCentr,pt);
+
+         // Extra factor if TOF matching was required
+         if(detUsed & AliPIDResponse::kDetTOF){
+           Float_t kaonTOFfactor = 0.1;
+           if(pt > 0.35) kaonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,4.19618E-07)/5.68017E-01)*TMath::Power(pt,-1.50705);
+           Float_t protonTOFfactor = 0.1;
+           if(pt > 0.4) protonTOFfactor = 1 - TMath::Exp(-TMath::Power(pt,3.30978)/8.57396E-02)*TMath::Power(pt,-4.42661E-01);
+           
+           if(track->Charge() < 0){ // for negative tracks
+             kaonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,4.87912E-07)/3.26431E-01)*TMath::Power(pt,-1.22893);
+             protonTOFfactor *= 1 - TMath::Exp(-TMath::Power(pt,2.00575E-07)/4.95605E-01)*TMath::Power(pt,-6.71305E-01);
+           }
+           
+           p[3] *= kaonTOFfactor;
+           p[4] *= protonTOFfactor;
+         }
+
+         return;
+       }
+
+
+       Double_t sumPriors = 0;
+       for (Int_t i=0;i<fSelectedSpecies;++i){
+               if (fPriorsDistributions[i]){
+                       p[i]=fPriorsDistributions[i]->Interpolate(pt);
+               }
+               else {
+                       p[i]=0.;
+               }
+               sumPriors+=p[i];                
+       }
 
+       // normalizing priors
+       if (sumPriors == 0) return;
+       for (Int_t i=0;i<fSelectedSpecies;++i){
+          p[i]=p[i]/sumPriors;
+       }
+       return;
+}
 //-------------------------------------------------------------------------------------------------    
 void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Double_t* probDensity, const Double_t* prior) const {
 
@@ -208,6 +306,7 @@ void AliPIDCombined::ComputeBayesProbabilities(Double_t* probabilities, const Do
     sum += probDensity[i] * prior[i];
   }
   if (sum <= 0) {
+
     AliError("Invalid probability densities or priors");
     for (Int_t i = 0; i < fSelectedSpecies; i++) probabilities[i] = -1;
     return;
@@ -235,5 +334,25 @@ void AliPIDCombined::SetCombinedStatus(AliPIDResponse::EDetPidStatus status, UIn
 
 }
 
-ClassImp(AliPIDCombined);
-
+//----------------------------------------------------------------------------------------
+void AliPIDCombined::SetDefaultTPCPriors(){
+  fEnablePriors=kTRUE;
+  fUseDefaultTPCPriors = kTRUE;
+
+  TString oadbfilename("$ALICE_ROOT/OADB/COMMON/PID/data/");
+  oadbfilename += "/PIDdefaultPriors.root";
+  TFile * foadb = TFile::Open(oadbfilename.Data());
+  if(!foadb->IsOpen()) AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
+  
+  AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
+  if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
+  
+  const char *namespecies[5] = {"El","Mu","Pi","Ka","Pr"};
+  char name[100];
+
+  for(Int_t i=0;i < 5;i++){
+    sprintf(name,"hDefault%sPriors",namespecies[i]);
+    fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
+    if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
+  }
+}

diff --git a/STEER/STEERBase/AliPIDCombined.h b/STEER/STEERBase/AliPIDCombined.h
index 446449fd03ee39afe98be5e72b32d36136c2f4c5..ce1a52b2ec7844d115348a7620b2d2cd9f916bb8 100644 (file)
--- a/STEER/STEERBase/AliPIDCombined.h
+++ b/STEER/STEERBase/AliPIDCombined.h
@@ -18,6 +18,7 @@
 #include <AliPID.h>\r
 #include <AliPIDResponse.h>\r
 #include <TH1F.h>\r
+#include <TH2F.h>\r
 \r
 //class TH1;\r
 class AliPIDResponse;\r
@@ -37,13 +38,17 @@ public:
   void SetPriorDistribution(AliPID::EParticleType type,TH1F *prior);\r
   //  const TH1* GetPriorDistribution(AliPID::EParticleType type) const {return (TH1*)fPriorsDistributions[type];}\r
   TH1* GetPriorDistribution(AliPID::EParticleType type)  const {return (TH1*)fPriorsDistributions[type];}\r
+  \r
+  void GetPriors(const AliVTrack *track,Double_t* p,const AliPIDResponse *response,UInt_t detUsed) const;\r
+  \r
+  void SetDefaultTPCPriors();\r
        \r
   UInt_t ComputeProbabilities(const AliVTrack *track, const AliPIDResponse *response, Double_t* bayesProbabilities) const;\r
   void SetSelectedSpecies(Int_t selectedSpecies) {fSelectedSpecies = selectedSpecies;}\r
   Int_t GetSelectedSpecies() const {return fSelectedSpecies;}\r
 \r
 protected:\r
-  void GetPriors(const AliVTrack *track,Double_t* priors) const;\r
+  void GetPriors(const AliVTrack *track,Double_t* priors,Float_t centrality=-1) const;\r
   void ComputeBayesProbabilities(Double_t* bayesProbabilities,const Double_t* probDensity, const Double_t* priors) const;\r
   void SetCombinedStatus(const AliPIDResponse::EDetPidStatus status,UInt_t *mask, const AliPIDResponse::EDetCode bit, Double_t* p) const;\r
 \r
@@ -56,8 +61,10 @@ private:
   Bool_t fEnablePriors;      // Enable bayesian PID (if kFALSE priors set flat)\r
   Int_t fSelectedSpecies;    // Number of selected species to study\r
   TH1F *fPriorsDistributions[AliPID::kSPECIES+AliPID::kSPECIESLN]; // priors\r
+  Bool_t fUseDefaultTPCPriors; // switch to use Defaul TPC Priors\r
+  static TH2F *fDefaultPriorsTPC[5]; // Default priors for TPC tracks\r
 \r
-  ClassDef(AliPIDCombined,1);\r
+  ClassDef(AliPIDCombined,2);\r
 };\r
 \r
 #endif\r

diff --git a/STEER/STEERBase/AliPIDResponse.cxx b/STEER/STEERBase/AliPIDResponse.cxx
index ce9e0e9de479514b1d0b89e134e4fce7df4e36c4..5f7a77bb95bf782e32cd132c7a7c639df56fb862 100644 (file)
--- a/STEER/STEERBase/AliPIDResponse.cxx
+++ b/STEER/STEERBase/AliPIDResponse.cxx
@@ -40,6 +40,8 @@
 
 #include "AliPIDResponse.h"
 
+#include "AliCentrality.h"
+
 ClassImp(AliPIDResponse);
 
 AliPIDResponse::AliPIDResponse(Bool_t isMC/*=kFALSE*/) :
@@ -69,7 +71,8 @@ fTRDPIDReference(0x0),
 fTOFTimeZeroType(kBest_T0),
 fTOFres(100.),
 fEMCALPIDParams(0x0),
-fCurrentEvent(0x0)
+fCurrentEvent(0x0),
+fCurrCentrality(0.0)
 {
   //
   // default ctor
@@ -120,7 +123,8 @@ fTRDPIDReference(0x0),
 fTOFTimeZeroType(AliPIDResponse::kBest_T0),
 fTOFres(100.),
 fEMCALPIDParams(0x0),
-fCurrentEvent(0x0)
+fCurrentEvent(0x0),
+fCurrCentrality(0.0)
 {
   //
   // copy ctor
@@ -523,7 +527,15 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
   
   //TOF resolution
   SetTOFResponse(event, (AliPIDResponse::EStartTimeType_t)fTOFTimeZeroType);
-  
+
+  // Get and set centrality
+  AliCentrality *centrality = event->GetCentrality();
+  if(centrality){
+    fCurrCentrality = centrality->GetCentralityPercentile("V0M");
+  }
+  else{
+    fCurrCentrality = -1;
+  }
 }
 
 //______________________________________________________________________________

diff --git a/STEER/STEERBase/AliPIDResponse.h b/STEER/STEERBase/AliPIDResponse.h
index d6afdb7e0f47d9045ee8d8748e5baed96c71549c..e1a8306fe3e989cfb801a6ed410a1cb2ceb8793e 100644 (file)
--- a/STEER/STEERBase/AliPIDResponse.h
+++ b/STEER/STEERBase/AliPIDResponse.h
@@ -80,6 +80,7 @@ public:
   void SetTRDslicesForPID(UInt_t slice1, UInt_t slice2) {fTRDslicesForPID[0]=slice1;fTRDslicesForPID[1]=slice2;}
   
   void SetOADBPath(const char* path) {fOADBPath=path;}
+  const char *GetOADBPath() const {return fOADBPath.Data();}
   void InitialiseEvent(AliVEvent *event, Int_t pass);
   void SetCurrentFile(const char* file) { fCurrentFile=file; }
 
@@ -88,9 +89,13 @@ public:
   // User settings for the MC period and reco pass
   void SetMCperiod(const char *mcPeriod) {fMCperiodUser=mcPeriod;}
   void SetRecoPass(Int_t recoPass)       {fRecoPassUser=recoPass;}
-  
+
+  // event info
+  Float_t GetCurrentCentrality() const {return fCurrCentrality;};
+
   AliPIDResponse(const AliPIDResponse &other);
   AliPIDResponse& operator=(const AliPIDResponse &other);
+
   
 protected:
   AliITSPIDResponse fITSResponse;    //PID response function of the ITS
@@ -130,6 +135,8 @@ private:
   TObjArray *fEMCALPIDParams;             //! EMCAL PID Params
 
   AliVEvent *fCurrentEvent;            //! event currently being processed
+
+  Float_t fCurrCentrality;             //! current centrality
   
   void ExecNewRun();
   
@@ -158,7 +165,7 @@ private:
   //
   void SetRecoInfo();
   
-  ClassDef(AliPIDResponse,3);  //PID response handling
+  ClassDef(AliPIDResponse,4);  //PID response handling
 };
 
 inline Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const {