- Reshuffling of the particle codes in AliPID. Now the light nuclei are between the
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 23 Aug 2012 16:08:32 +0000 (16:08 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 23 Aug 2012 16:08:32 +0000 (16:08 +0000)
  five standard charged particles and the neutral
- Added basic support for light nuclei in ITS, TPC, TOF (to be refined)
- Added caching of PID information
- code cleanup
Jens.Wiechula@cern.ch
pietro.antonioli@bo.infn.it

48 files changed:
ANALYSIS/AliAnalysisTaskESDfilter.cxx
ANALYSIS/AliAnalysisTaskPIDResponse.cxx
ANALYSIS/AliAnalysisTaskPIDResponse.h
ANALYSIS/AliAnalysisTaskPIDqa.cxx
ANALYSIS/TenderSupplies/AddTaskTender.C
ANALYSIS/TenderSupplies/AliPIDTenderSupply.cxx
ANALYSIS/TenderSupplies/AliPIDTenderSupply.h
CORRFW/AliCFTrackCutPid.cxx
CORRFW/AliCFTrackCutPid.h
EMCAL/AliEMCALPID.cxx
EMCAL/AliEMCALPIDUtils.cxx
EMCAL/AliEMCALPIDUtils.h
EMCAL/AliEMCALRecoUtils.cxx
HLT/BASE/util/AliHLTCaloClusterDataStruct.h
HLT/CALO/AliHLTCaloClusterAnalyser.cxx
HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx
PHOS/AliPHOSFastRecParticle.cxx
PHOS/AliPHOSFastRecParticle.h
PHOS/AliPHOSPIDv1.cxx
PHOS/AliPHOSPIDv1.h
PHOS/AliPHOSRecParticle.cxx
PWG/CaloTrackCorrBase/AliCaloPID.cxx
PWGHF/hfe/AliHFEdca.cxx
STEER/AOD/AliAODPid.cxx
STEER/AOD/AliAODPid.h
STEER/AOD/AliAODTrack.cxx
STEER/AOD/AliAODTrack.h
STEER/AOD/AliAODpidUtil.cxx
STEER/AOD/AliAODpidUtil.h
STEER/CMakelibSTEERBase.pkg
STEER/ESD/AliESDCaloCluster.cxx
STEER/ESD/AliESDCaloCluster.h
STEER/ESD/AliESDpid.cxx
STEER/ESD/AliESDpid.h
STEER/ESD/AliESDtrack.cxx
STEER/ESD/AliESDtrack.h
STEER/STEERBase/AliEMCALPIDResponse.cxx
STEER/STEERBase/AliEMCALPIDResponse.h
STEER/STEERBase/AliITSPIDResponse.h
STEER/STEERBase/AliPID.cxx
STEER/STEERBase/AliPID.h
STEER/STEERBase/AliPIDCombined.cxx
STEER/STEERBase/AliPIDCombined.h
STEER/STEERBase/AliPIDResponse.cxx
STEER/STEERBase/AliPIDResponse.h
STEER/STEERBase/AliTPCPIDResponse.cxx
STEER/STEERBase/AliVTrack.h
STEER/STEERBaseLinkDef.h

index 4dac86a..90e4ae6 100644 (file)
@@ -2369,7 +2369,7 @@ void AliAnalysisTaskESDfilter::SetDetectorRawSignals(AliAODPid *aodpid, AliESDtr
  aodpid->SetTRDChi2(track->GetTRDchi2());
 
  //TOF PID  
- Double_t times[AliAODPid::kSPECIES]; track->GetIntegratedTimes(times);
+ Double_t times[AliPID::kSPECIES]; track->GetIntegratedTimes(times);
  aodpid->SetIntegratedTimes(times);
 
    //  Float_t tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
index 4ff0d91..f0a2f08 100644 (file)
@@ -27,7 +27,6 @@
 #include <AliPIDResponse.h>
 #include <AliESDpid.h>
 
-
 #include "AliAnalysisTaskPIDResponse.h"
 
 ClassImp(AliAnalysisTaskPIDResponse)
@@ -36,6 +35,7 @@ ClassImp(AliAnalysisTaskPIDResponse)
 AliAnalysisTaskPIDResponse::AliAnalysisTaskPIDResponse():
 AliAnalysisTaskSE(),
 fIsMC(kFALSE),
+fCachePID(kTRUE),
 fOADBPath(),
 fPIDResponse(0x0),
 fRun(0),
@@ -53,6 +53,7 @@ fRecoPassTuned(0)
 AliAnalysisTaskPIDResponse::AliAnalysisTaskPIDResponse(const char* name):
 AliAnalysisTaskSE(name),
 fIsMC(kFALSE),
+fCachePID(kTRUE),
 fOADBPath(),
 fPIDResponse(0x0),
 fRun(0),
@@ -81,8 +82,9 @@ void AliAnalysisTaskPIDResponse::UserCreateOutputObjects()
   //
   // Create the output QA objects
   //
-    
+  
   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDResponse",10);
+  
   //input hander
   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
@@ -118,6 +120,10 @@ void AliAnalysisTaskPIDResponse::UserExec(Option_t */*option*/)
   if(pidresp && AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
       pidresp->SetEventHandler(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
   }
+  //create and attach transient PID object
+  if (fCachePID) {
+    fPIDResponse->FillTrackDetectorPID();
+  }
 }
 
 //______________________________________________________________________________
index 2894ae7..71cdd66 100644 (file)
@@ -32,7 +32,9 @@ public:
   AliAnalysisTaskPIDResponse(const char *name);
   virtual ~AliAnalysisTaskPIDResponse();
 
-  void SetIsMC(Bool_t isMC=kTRUE) { fIsMC=isMC; }
+  void SetIsMC(Bool_t isMC=kTRUE)   { fIsMC=isMC; }
+  void SetCachePID(Bool_t cachePID) { fCachePID=cachePID; }
+  Bool_t GetCachePID() const { return fCachePID; }
   
   virtual void UserCreateOutputObjects();
   
@@ -43,7 +45,8 @@ public:
   void SetTuneOnData(Bool_t flag,Int_t recopass){fIsTunedOnData=flag;fRecoPassTuned=recopass;};
 
 private:
-  Bool_t fIsMC;                        //  If we run on MC data
+  Bool_t fIsMC;                        // If we run on MC data
+  Bool_t fCachePID;                    // Cache PID values in transient object
   TString fOADBPath;                   // OADB path to use
   
   AliPIDResponse *fPIDResponse;        //! PID response Handler
index e1eff32..54b8ff1 100644 (file)
@@ -381,13 +381,13 @@ void AliAnalysisTaskPIDqa::FillITSqa()
     }\r
     \r
     \r
-    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
       TH2 *h=(TH2*)theList->At(ispecie);\r
       if (!h) continue;\r
       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);\r
       h->Fill(mom,nSigma);\r
     }\r
-    TH2 *h=(TH2*)theList->At(AliPID::kSPECIES);\r
+    TH2 *h=(TH2*)theList->At(AliPID::kSPECIESC);\r
     if (h) {\r
       Double_t sig=track->GetITSsignal();\r
       h->Fill(mom,sig);\r
@@ -429,14 +429,14 @@ void AliAnalysisTaskPIDqa::FillTPCqa()
     \r
     Double_t mom=track->GetTPCmomentum();\r
     \r
-    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
       TH2 *h=(TH2*)fListQAtpc->At(ispecie);\r
       if (!h) continue;\r
       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);\r
       h->Fill(mom,nSigma);\r
     }\r
     \r
-    TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);\r
+    TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIESC);\r
     if (h) {\r
       Double_t sig=track->GetTPCsignal();\r
       h->Fill(mom,sig);\r
@@ -529,7 +529,7 @@ void AliAnalysisTaskPIDqa::FillTOFqa()
 \r
     Double_t mom=track->P();\r
 \r
-    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
       TH2 *h=(TH2*)fListQAtof->At(ispecie);\r
       if (!h) continue;\r
       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);\r
@@ -875,7 +875,7 @@ void AliAnalysisTaskPIDqa::FillTPCTOFqa()
     Double_t mom=track->P();\r
     Double_t momTPC=track->GetTPCmomentum();\r
 \r
-    for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+    for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
       //TOF nSigma\r
       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);\r
       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);\r
@@ -885,11 +885,11 @@ void AliAnalysisTaskPIDqa::FillTPCTOFqa()
       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);\r
 \r
       //TOF after TPC cut\r
-      h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);\r
+      h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);\r
       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);\r
 \r
       //EMCAL after TOF and TPC cut\r
-      h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);\r
+      h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);\r
       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){\r
 \r
        Int_t nMatchClus = track->GetEMCALcluster();\r
@@ -929,7 +929,7 @@ void AliAnalysisTaskPIDqa::SetupITSqa()
   TVectorD *vX=MakeLogBinning(200,.1,30);\r
   \r
   //ITS+TPC tracks\r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),\r
                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
@@ -943,7 +943,7 @@ void AliAnalysisTaskPIDqa::SetupITSqa()
   fListQAits->Add(hSig);\r
 \r
   //ITS Standalone tracks\r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),\r
                                Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                                vX->GetNrows()-1,vX->GetMatrixArray(),\r
@@ -957,7 +957,7 @@ void AliAnalysisTaskPIDqa::SetupITSqa()
   fListQAitsSA->Add(hSigSA);\r
   \r
   //ITS Pure Standalone tracks\r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),\r
                                    Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                                    vX->GetNrows()-1,vX->GetMatrixArray(),\r
@@ -982,7 +982,7 @@ void AliAnalysisTaskPIDqa::SetupTPCqa()
   \r
   TVectorD *vX=MakeLogBinning(200,.1,30);\r
   \r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),\r
                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
@@ -1028,7 +1028,7 @@ void AliAnalysisTaskPIDqa::SetupTOFqa()
   \r
   TVectorD *vX=MakeLogBinning(200,.1,30);\r
 \r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),\r
                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
@@ -1158,7 +1158,7 @@ void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
   TVectorD *vX=MakeLogBinning(200,.1,30);\r
 \r
   //TPC signals after TOF cut\r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),\r
                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
@@ -1167,7 +1167,7 @@ void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
   }\r
 \r
   //TOF signals after TPC cut\r
-  for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){\r
+  for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){\r
     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),\r
                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),\r
                               vX->GetNrows()-1,vX->GetMatrixArray(),\r
index e1d1575..2831f2b 100644 (file)
@@ -30,15 +30,23 @@ AliAnalysisTask *AddTaskTender(Bool_t useV0=kFALSE,
   tender->SetCheckEventSelection(checkEvtSelection);
   tender->SetDefaultCDBStorage("raw://");
   mgr->AddTask(tender);
-  if (checkEvtSelection) {
-    if (mgr->GetTasks()->First() != (TObject*)tender ) {
-      TString firstName(mgr->GetTasks()->First()->GetName());
-      Bool_t isSecond=(mgr->GetTasks()->At(1) == (TObject*)tender);
-      if (! (firstName=="PIDResponseTask" && isSecond )){
-        Fatal("AddTaskTender","When setting the tender to check the event selection, it has to be the first wagon, or the first after the PID Response!");
-        return NULL;
-      }
-    }
+  
+  Bool_t cachePID=kFALSE;
+
+  //check that that tender is the first task after the pid response
+  TString firstName(mgr->GetTasks()->First()->GetName());
+  Bool_t isSecond=(mgr->GetTasks()->At(1) == (TObject*)tender);
+
+  if (! (firstName=="PIDResponseTask" && isSecond )){
+    Fatal("AddTaskTender","When using the tender the first task needs to be the PIDResponse and the tender the second task!!!");
+    return NULL;
+  }
+  
+  AliAnalysisTaskPIDResponse *pidResp=(AliAnalysisTaskPIDResponse*)mgr->GetTasks()->First();
+  if (pidResp->GetCachePID()){
+    usePID=kTRUE;
+    pidResp->SetCachePID(kFALSE);
+    cachePID=kTRUE;
   }
   
   //========= Attach VZERO supply ======
@@ -105,11 +113,6 @@ AliAnalysisTask *AddTaskTender(Bool_t useV0=kFALSE,
     tender->AddSupply(trdSupply);
   }  
 
-  //========= Attach PID supply ======
-  if (usePID) {
-    tender->AddSupply(new AliPIDTenderSupply("PIDtender"));
-  }  
-
   //========= Attach Primary Vertex supply ======
   if (useVTX) {
     tender->AddSupply(new AliVtxTenderSupply("PriVtxtender"));
@@ -122,6 +125,13 @@ AliAnalysisTask *AddTaskTender(Bool_t useV0=kFALSE,
     tender->AddSupply(emcSupply);
   }  
 
+  //========= Attach PID supply ======
+  if (usePID) {
+    AliPIDTenderSupply *pidSupply=new AliPIDTenderSupply("PIDtender");
+    pidSupply->SetCachePID(cachePID);
+    tender->AddSupply(pidSupply);
+  }
+
   //================================================
   //              data containers
   //================================================
index 6b567c9..cadf213 100644 (file)
 
 #include "AliPIDTenderSupply.h"
 
+ClassImp(AliPIDTenderSupply)
+
 AliPIDTenderSupply::AliPIDTenderSupply() :
-  AliTenderSupply()
+  AliTenderSupply(),
+  fCachePID(kFALSE)
 {
   //
   // default ctor
@@ -39,7 +42,8 @@ AliPIDTenderSupply::AliPIDTenderSupply() :
 
 //_____________________________________________________
 AliPIDTenderSupply::AliPIDTenderSupply(const char *name, const AliTender *tender) :
-  AliTenderSupply(name,tender)
+  AliTenderSupply(name,tender),
+  fCachePID(kFALSE)
 {
   //
   // named ctor
@@ -58,6 +62,11 @@ void AliPIDTenderSupply::ProcessEvent()
 
   AliESDpid *pid=fTender->GetESDhandler()->GetESDpid();
   if (!pid) return;
+  // chache pid if requested
+  if (fCachePID) {
+    pid->FillTrackDetectorPID();
+  }
+  
   //
   // recalculate combined PID probabilities
   //
index a82dc49..59d0330 100644 (file)
@@ -24,13 +24,15 @@ public:
   
   virtual void              Init(){;}
   virtual void              ProcessEvent();
-  
+
+  void SetCachePID(Bool_t cachePID) { fCachePID=cachePID; }
 private:
+  Bool_t fCachePID;                    // Cache PID values in transient object
   
   AliPIDTenderSupply(const AliPIDTenderSupply&c);
   AliPIDTenderSupply& operator= (const AliPIDTenderSupply&c);
   
-  ClassDef(AliPIDTenderSupply, 1);  // PID tender task
+  ClassDef(AliPIDTenderSupply, 2);  // PID tender task
 };
 
 
index a57cdf8..7d94b0e 100644 (file)
@@ -90,7 +90,7 @@ AliCFTrackCutPid::AliCFTrackCutPid() :
   //
   //Default constructor 
   //
-  for(Int_t j=0; j< AliPID::kSPECIESN; j++) {
+  for(Int_t j=0; j< AliPID::kSPECIES; j++) {
     fPriors[j]=0.2;
   }
   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
@@ -127,7 +127,7 @@ AliCFTrackCutPid::AliCFTrackCutPid(const Char_t* name, const Char_t* title) :
   //
   //Constructor
   // 
-  for(Int_t j=0; j< AliPID::kSPECIESN; j++) {
+  for(Int_t j=0; j< AliPID::kSPECIES; j++) {
     fPriors[j]=0.2;
   }
   for(Int_t j=0; j< AliPID::kSPECIES; j++) {
@@ -172,7 +172,7 @@ AliCFTrackCutPid::AliCFTrackCutPid(const AliCFTrackCutPid& c) :
       fhProb[i][iP]=c.fhProb[i][iP];
     }
   }
-  for(Int_t j=0; j< AliPID::kSPECIESN; j++){
+  for(Int_t j=0; j< AliPID::kSPECIES; j++){
     fPriors[j]=c.fPriors[j];
   }
   for(Int_t j=0; j< AliPID::kSPECIES; j++){
@@ -215,7 +215,7 @@ AliCFTrackCutPid& AliCFTrackCutPid::operator=(const AliCFTrackCutPid& c)
       }  
     }
 
-    for(Int_t j=0; j< AliPID::kSPECIESN; j++){
+    for(Int_t j=0; j< AliPID::kSPECIES; j++){
       this->fPriors[j]=c.fPriors[j];
     }
     for(Int_t j=0; j< AliPID::kSPECIES; j++){
index 3e82fcc..2b4df8d 100644 (file)
@@ -113,7 +113,7 @@ class AliCFTrackCutPid : public AliCFCutBase
                                                             // fgParticleType if the probability is larger than this threshold,
                                                             // regardless it is the highest or not (!)
 
-  Double_t fPriors[AliPID::kSPECIESN];                       // a priori concentrations
+  Double_t fPriors[AliPID::kSPECIES];                       // a priori concentrations
   TF1 *fPriorsFunc[AliPID::kSPECIES];                       // momentum dependent priors
   Bool_t fDets[kNdets];                                     // boolean(s) corresponding to the chosen detector(s) 
   Bool_t fDetsInAnd[kNdets];                                // detector to be in AND for the combined PID
index 71cd663..ec26872 100644 (file)
@@ -32,7 +32,7 @@
 //
 //   pid->GetPIDFinal(idx) gives the probabilities
 //
-//   Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
+//   Double_t PIDFinal[AliPID::kSPECIESCN]  is the standard PID for :
 //
 //     kElectron :  fPIDFinal[0]
 //     kMuon     :  fPIDFinal[1]
@@ -177,7 +177,7 @@ void AliEMCALPID::InitParameters()
   fPIDWeight[1] = -1;
   fPIDWeight[2] = -1;
   
-  for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
+  for(Int_t i=0; i<AliPID::kSPECIESCN+1; i++)
     fPIDFinal[i]= 0;
   
   const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
index bae8bcc..6c1c28e 100644 (file)
@@ -31,7 +31,7 @@
 //
 //   pid->GetPIDFinal(idx) gives the probabilities
 //
-//   Double_t PIDFinal[AliPID::kSPECIESN]  is the standard PID for :
+//   Double_t PIDFinal[AliPID::kSPECIESCN]  is the standard PID for :
 //
 //     kElectron :  fPIDFinal[0]
 //     kMuon     :  fPIDFinal[1]
@@ -164,19 +164,27 @@ void AliEMCALPIDUtils::ComputePID(Double_t energy, Double_t lambda0)
     AliInfo(Form( "PIDWeight in loop = %f ||| %f ||| %f",  fPIDWeight[0] , fPIDWeight[1], fPIDWeight[2]) );
     AliInfo("********************************************************" );
   }
-  
-  fPIDFinal[0]  = fPIDWeight[0]/2; // photon
-  fPIDFinal[1]  = fPIDWeight[2]/8;
-  fPIDFinal[2]  = fPIDWeight[2]/8;
-  fPIDFinal[3]  = fPIDWeight[2]/8;
-  fPIDFinal[4]  = fPIDWeight[2]/8;
-  fPIDFinal[5]  = fPIDWeight[0]/2; // electron
-  fPIDFinal[6]  = fPIDWeight[1]  ; // Pi0
-  fPIDFinal[7]  = fPIDWeight[2]/8;
-  fPIDFinal[8]  = fPIDWeight[2]/8;
-  fPIDFinal[9]  = fPIDWeight[2]/8;
-  fPIDFinal[10] = fPIDWeight[2]/8;
 
+  //default particles
+  fPIDFinal[AliPID::kElectron]  = fPIDWeight[0]/2; // photon
+  fPIDFinal[AliPID::kMuon]      = fPIDWeight[2]/8;
+  fPIDFinal[AliPID::kPion]      = fPIDWeight[2]/8;
+  fPIDFinal[AliPID::kKaon]      = fPIDWeight[2]/8;
+  fPIDFinal[AliPID::kProton]    = fPIDWeight[2]/8;
+  //light nuclei
+  fPIDFinal[AliPID::kDeuteron]  = 0;
+  fPIDFinal[AliPID::kTriton]    = 0;
+  fPIDFinal[AliPID::kHe3]       = 0;
+  fPIDFinal[AliPID::kAlpha]     = 0;
+  //neutral particles
+  fPIDFinal[AliPID::kPhoton]    = fPIDWeight[0]/2; // electron
+  fPIDFinal[AliPID::kPi0]       = fPIDWeight[1]  ; // Pi0
+  fPIDFinal[AliPID::kNeutron]   = fPIDWeight[2]/8;
+  fPIDFinal[AliPID::kKaon0]     = fPIDWeight[2]/8;
+  fPIDFinal[AliPID::kEleCon]    = fPIDWeight[2]/8;
+  //
+  fPIDFinal[AliPID::kUnknown]   = fPIDWeight[2]/8;
+  
 }
 
 
@@ -389,7 +397,7 @@ void AliEMCALPIDUtils::InitParameters()
   fPIDWeight[1] = -1;
   fPIDWeight[2] = -1;
   
-  for(Int_t i=0; i<AliPID::kSPECIESN+1; i++)
+  for(Int_t i=0; i<AliPID::kSPECIESCN+1; i++)
     fPIDFinal[i]= 0;
   
   //   init the parameters here instead of from loading from recparam
index d2b648c..eaa8b32 100644 (file)
@@ -34,11 +34,11 @@ public:
   Double_t DistEnergy(const Double_t energy, const Int_t nature) ;
 
   Double_t GetPID(Int_t idx) const {if (idx>=0&&idx<3) return fPID[idx]; else return 0.;}
-  Double_t GetPIDFinal(Int_t idx) const {if (idx>=0&&idx<AliPID::kSPECIESN) return fPIDFinal[idx]; else return 0.;}
+  Double_t GetPIDFinal(Int_t idx) const {if (idx>=0&&idx<AliPID::kSPECIESCN) return fPIDFinal[idx]; else return 0.;}
   Double_t GetPIDWeight(Int_t idx) const {if (idx>=0&&idx<3) return fPIDWeight[idx]; else return 0.;}
   
   void    SetPID(Double_t val, Int_t idx) {if (idx>=0&&idx<3) fPID[idx] = val;}
-  void    SetPIDFinal(Double_t val, Int_t idx) {if (idx>=0&&idx<AliPID::kSPECIESN) fPIDFinal[idx] = val;}
+  void    SetPIDFinal(Double_t val, Int_t idx) {if (idx>=0&&idx<AliPID::kSPECIESCN) fPIDFinal[idx] = val;}
   void    SetPIDWeight(Double_t val, Int_t idx) {if (idx>=0&&idx<3) fPIDWeight[idx] = val;}
   void    SetPrintInfo(Bool_t yesno) {fPrintInfo = yesno;}
        
@@ -67,7 +67,7 @@ protected:
    
   Float_t fPID[3];
   
-  Float_t  fPIDFinal[AliPID::kSPECIESN+1]; // final PID format
+  Float_t  fPIDFinal[AliPID::kSPECIESCN+1]; // final PID format
   Float_t  fPIDWeight[3];                  // order: gamma, pi0, hadrons,
   Double_t fProbGamma;                   // probility to be a Gamma
   Double_t fProbPiZero;                          // probility to be a PiO
index dac3c43..56cbc67 100644 (file)
@@ -1528,8 +1528,8 @@ void AliEMCALRecoUtils::RecalculateClusterPID(AliVCluster * cluster)
   if ( cluster->GetM02() != 0)
     fPIDUtils->ComputePID(cluster->E(),cluster->GetM02());
   
-  Float_t pidlist[AliPID::kSPECIESN+1];
-  for(Int_t i = 0; i < AliPID::kSPECIESN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
+  Float_t pidlist[AliPID::kSPECIESCN+1];
+  for(Int_t i = 0; i < AliPID::kSPECIESCN+1; i++) pidlist[i] = fPIDUtils->GetPIDFinal(i);
         
   cluster->SetPID(pidlist);
 }
index ab25cc9..fefa58d 100644 (file)
@@ -131,7 +131,7 @@ struct AliHLTCaloClusterDataStruct
   // In case all the weights are non-positive they are replaced by
   // uniform probabilities
 
-  Int_t n = AliPID::kSPECIESN;
+  Int_t n = AliPID::kSPECIESCN;
 
   Float_t uniform = 1./(Float_t)n;
 
@@ -283,7 +283,7 @@ struct AliHLTCaloClusterDataStruct
   Float_t fTrackDz; //COMMENT
 
   /** PID */
-  Float_t fPID[AliPID::kSPECIESN];            //COMMENT
+  Float_t fPID[AliPID::kSPECIESCN];            //COMMENT
 
   /** Unique ID of the cluster*/
   Int_t fID;                                     //COMMENT
index 04267b6..66fee31 100644 (file)
@@ -350,7 +350,7 @@ AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize
        }
       else
        {
-         for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
+         for(Int_t k = 0; k < AliPID::kSPECIESCN; k++)
            {
              caloClusterPtr->fPID[k] = 0;
            }
index a0d18a5..d1302d6 100644 (file)
@@ -220,7 +220,7 @@ AliHLTPHOSClusterAnalyser::CreateClusters(UInt_t availableSize, UInt_t& totSize)
        }
       else
        {
-         for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
+         for(Int_t k = 0; k < AliPID::kSPECIESCN; k++)
            {
              caloClusterPtr->fPID[k] = 0;
            }
index 2bea0db..6ef6bfb 100644 (file)
@@ -52,7 +52,7 @@ AliPHOSFastRecParticle::AliPHOSFastRecParticle() :
 {
   // ctor
 
-  for(Int_t i=0; i<AliPID::kSPECIESN; i++) {
+  for(Int_t i=0; i<AliPID::kSPECIESCN; i++) {
     fPID[i] = -111.;
   }
   
@@ -86,7 +86,7 @@ AliPHOSFastRecParticle::AliPHOSFastRecParticle() :
   fPolarPhi    = rp.fPolarPhi;
   fParticlePDG = rp.fParticlePDG;
 
-  for(Int_t i=0; i<AliPID::kSPECIESN; i++) {
+  for(Int_t i=0; i<AliPID::kSPECIESCN; i++) {
     fPID[i] = rp.fPID[i];
   }
  
@@ -121,7 +121,7 @@ AliPHOSFastRecParticle::AliPHOSFastRecParticle() :
   fPolarPhi    = p.fPolarPhi;
   fParticlePDG = p.fParticlePDG; 
   
-  for(Int_t i=0; i<AliPID::kSPECIESN; i++) {
+  for(Int_t i=0; i<AliPID::kSPECIESCN; i++) {
     fPID[i] = p.fPID[i];
   }
   
index 1570dca..d9a55e1 100644 (file)
@@ -105,7 +105,7 @@ class AliPHOSFastRecParticle : public TParticle {
   Int_t fIndexInList ; // the index of this RecParticle in the list stored in TreeR (to be set by analysis)
   Float_t fTof ;       // time of fliht
   Int_t fType ;        // particle type obtained by "virtual" reconstruction
-  Float_t fPID[AliPID::kSPECIESN] ; // PID probability densities
+  Float_t fPID[AliPID::kSPECIESCN] ; // PID probability densities
 
  private:
 
index 86b1598..c333e1b 100644 (file)
@@ -473,7 +473,7 @@ void AliPHOSPIDv1::InitParameters()
   fERecWeight ->SetParameters(fERecWeightPar[0],fERecWeightPar[1] ,fERecWeightPar[2] ,fERecWeightPar[3]) ; 
 
 
-  for (Int_t i =0; i<  AliPID::kSPECIESN ; i++)
+  for (Int_t i =0; i<  AliPID::kSPECIESCN ; i++)
     fInitPID[i] = 1.;
  
 }
@@ -944,7 +944,7 @@ void  AliPHOSPIDv1::MakePID()
 {
   // construct the PID weight from a Bayesian Method
   
-  const Int_t kSPECIES = AliPID::kSPECIESN ;
+  const Int_t kSPECIES = AliPID::kSPECIESCN ;
  
   Int_t nparticles = fRecParticles->GetEntriesFast() ;
 
@@ -1632,10 +1632,10 @@ void AliPHOSPIDv1::GetVertex(void)
 //_______________________________________________________________________
 void AliPHOSPIDv1::SetInitPID(const Double_t *p) {
   // Sets values for the initial population of each particle type 
-  for (Int_t i=0; i<AliPID::kSPECIESN; i++) fInitPID[i] = p[i];
+  for (Int_t i=0; i<AliPID::kSPECIESCN; i++) fInitPID[i] = p[i];
 }
 //_______________________________________________________________________
 void AliPHOSPIDv1::GetInitPID(Double_t *p) const {
   // Gets values for the initial population of each particle type 
-  for (Int_t i=0; i<AliPID::kSPECIESN; i++) p[i] = fInitPID[i];
+  for (Int_t i=0; i<AliPID::kSPECIESCN; i++) p[i] = fInitPID[i];
 }
index b0e86ba..3e7ff28 100644 (file)
@@ -155,7 +155,7 @@ private:
   TVector3   fVtx ;                      //! Vertex position in current event
 
   //Initial pid population
-  Double_t fInitPID[AliPID::kSPECIESN] ; // Initial population to do bayesian PID
+  Double_t fInitPID[AliPID::kSPECIESCN] ; // Initial population to do bayesian PID
   // pid probability function parameters
   // ToF
   Double_t fTphoton[3] ;       // gaussian tof response for photon
index e9fd03e..f4ce52b 100644 (file)
@@ -46,7 +46,7 @@ AliPHOSRecParticle::AliPHOSRecParticle():
   fPos()
 {
   // ctor
-  const Int_t nSPECIES = AliPID::kSPECIESN;
+  const Int_t nSPECIES = AliPID::kSPECIESCN;
   for(Int_t i = 0; i<nSPECIES ; i++)
     fPID[i]=0.;
 }
@@ -82,7 +82,7 @@ AliPHOSRecParticle::AliPHOSRecParticle(const AliPHOSRecParticle & rp):
   fPolarTheta  = rp.fPolarTheta;
   fPolarPhi    = rp.fPolarPhi;
   fParticlePDG = rp.fParticlePDG; 
-  const Int_t nSPECIES = AliPID::kSPECIESN;
+  const Int_t nSPECIES = AliPID::kSPECIESCN;
   for(Int_t i = 0; i<nSPECIES ; i++)
     fPID[i]=rp.fPID[i];
 }
index 21cbe65..3ec4703 100755 (executable)
@@ -253,16 +253,16 @@ Int_t AliCaloPID::GetIdentifiedParticleType(const AliVCluster * cluster)
   
   if(fUseBayesianWeights)
   {
-    Double_t weights[AliPID::kSPECIESN];
+    Double_t weights[AliPID::kSPECIESCN];
     
     if(cluster->IsEMCAL() && fRecalculateBayesian)
     {          
       fEMCALPIDUtils->ComputePID(energy, lambda0);
-      for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
+      for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = fEMCALPIDUtils->GetPIDFinal(i);
     }
     else 
     {
-      for(Int_t i = 0; i < AliPID::kSPECIESN; i++) weights[i] = cluster->GetPID()[i];
+      for(Int_t i = 0; i < AliPID::kSPECIESCN; i++) weights[i] = cluster->GetPID()[i];
     }
 
     if(fDebug > 0)  PrintClusterPIDWeights(weights);
index 5223146..1f27dec 100644 (file)
@@ -1286,7 +1286,7 @@ Int_t AliHFEdca::GetCombinedPid(const AliESDtrack *const track)
  track->GetESDpid(prob);
 
  // setting priors!
- Double_t priors[AliPID::kSPECIESN];
+ Double_t priors[AliPID::kSPECIES];
  priors[0] = 0.01;
  priors[1] = 0.01;
  priors[2] = 0.85;
index 81a3857..321fcc4 100644 (file)
@@ -43,7 +43,7 @@ AliAODPid::AliAODPid():
     fTPCdEdxInfo(0)
 {
   // default constructor
-    for(Int_t i=0; i<kSPECIES; i++) fIntTime[i]   = 0; 
+    for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntTime[i]   = 0; 
     for(Int_t i=0; i<5; i++) fTOFpidResolution[i] = 0.;
     for(Int_t i=0; i<6; i++) {
       fTRDmomentum[i]      = 0.;
@@ -78,7 +78,7 @@ AliAODPid::AliAODPid(const AliAODPid& pid) :
 {
   // Copy constructor
   SetTRDsignal(fTRDnSlices, pid.fTRDslices);
-    for(Int_t i=0; i<kSPECIES; i++) fIntTime[i]=pid.fIntTime[i];
+    for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntTime[i]=pid.fIntTime[i];
 
     for(Int_t i=0; i<6; i++){ 
       fTRDmomentum[i]=pid.fTRDmomentum[i];
@@ -136,16 +136,16 @@ AliAODPid& AliAODPid::operator=(const AliAODPid& pid)
   return *this;
 }
 //_______________________________________________________________________________
-void AliAODPid::GetIntegratedTimes(Double_t timeint[kSPECIES]) const
+void AliAODPid::GetIntegratedTimes(Double_t timeint[AliPID::kSPECIES]) const
 {
  // Returns the array with integrated times for each particle hypothesis
-for(Int_t i=0; i<kSPECIES; i++) timeint[i]=fIntTime[i];
+ for(Int_t i=0; i<AliPID::kSPECIES; i++) timeint[i]=fIntTime[i];
 }
 //_______________________________________________________________________________
-void AliAODPid::SetIntegratedTimes(Double_t timeint[kSPECIES])
+void AliAODPid::SetIntegratedTimes(Double_t timeint[AliPID::kSPECIES])
 {
  // Returns the array with integrated times for each particle hypothesis
-for(Int_t i=0; i<kSPECIES; i++) fIntTime[i]=timeint[i];
+ for(Int_t i=0; i<AliPID::kSPECIES; i++) fIntTime[i]=timeint[i];
 }
 //______________________________________________________________________________
 void AliAODPid::SetTOFpidResolution(Double_t tofPIDres[5])
index b44c7ee..30c834b 100644 (file)
@@ -22,7 +22,7 @@ class AliAODPid : public TObject {
   AliAODPid(const AliAODPid& pid); 
   AliAODPid& operator=(const AliAODPid& pid);
   
-  enum{kSPECIES=5, kTRDnPlanes=6};
+  enum{kTRDnPlanes=6};
 
  //setters
   void      SetITSsignal(Double_t its)                         {fITSsignal=its;}
index b24d177..2b52a86 100644 (file)
@@ -26,6 +26,7 @@
 #include "AliExternalTrackParam.h"
 #include "AliVVertex.h"
 #include "AliAODTrack.h"
+#include "AliDetectorPID.h"
 #include "AliAODEvent.h"
 
 ClassImp(AliAODTrack)
@@ -53,6 +54,7 @@ AliAODTrack::AliAODTrack() :
   fCaloIndex(kEMCALNoMatch),
   fCovMatrix(NULL),
   fDetPid(NULL),
+  fDetectorPID(NULL),
   fProdVertex(NULL),
   fTrackPhiOnEMCal(-999),
   fTrackEtaOnEMCal(-999),
@@ -106,6 +108,7 @@ AliAODTrack::AliAODTrack(Short_t id,
   fCaloIndex(kEMCALNoMatch),
   fCovMatrix(NULL),
   fDetPid(NULL),
+  fDetectorPID(NULL),
   fProdVertex(prodVertex),
   fTrackPhiOnEMCal(-999),
   fTrackEtaOnEMCal(-999),
@@ -163,6 +166,7 @@ AliAODTrack::AliAODTrack(Short_t id,
   fCaloIndex(kEMCALNoMatch),
   fCovMatrix(NULL),
   fDetPid(NULL),
+  fDetectorPID(NULL),
   fProdVertex(prodVertex),
   fTrackPhiOnEMCal(-999),
   fTrackEtaOnEMCal(-999),
@@ -188,6 +192,7 @@ AliAODTrack::~AliAODTrack()
   // destructor
   delete fCovMatrix;
   delete fDetPid;
+  delete fDetectorPID;
 }
 
 
@@ -214,6 +219,7 @@ AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
   fCaloIndex(trk.fCaloIndex),
   fCovMatrix(NULL),
   fDetPid(NULL),
+  fDetectorPID(NULL),
   fProdVertex(trk.fProdVertex),
   fTrackPhiOnEMCal(trk.fTrackPhiOnEMCal),
   fTrackEtaOnEMCal(trk.fTrackEtaOnEMCal),
@@ -231,6 +237,7 @@ AliAODTrack::AliAODTrack(const AliAODTrack& trk) :
   if(trk.fCovMatrix) fCovMatrix=new AliAODRedCov<6>(*trk.fCovMatrix);
   if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
   SetPID(trk.fPID);
+  if (trk.fDetectorPID) fDetectorPID = new AliDetectorPID(*trk.fDetectorPID);
 }
 
 //______________________________________________________________________________
@@ -277,11 +284,15 @@ AliAODTrack& AliAODTrack::operator=(const AliAODTrack& trk)
     SetUsedForVtxFit(trk.GetUsedForVtxFit());
     SetUsedForPrimVtxFit(trk.GetUsedForPrimVtxFit());
 
+    //detector raw signals
     delete fDetPid;
     if(trk.fDetPid) fDetPid=new AliAODPid(*trk.fDetPid);
     else fDetPid=NULL;
 
-
+    //calibrated PID cache
+    delete fDetectorPID;
+    fDetectorPID=0x0;
+    if (trk.fDetectorPID) fDetectorPID = new AliDetectorPID(*trk.fDetectorPID);
   }
 
   return *this;
@@ -773,6 +784,16 @@ Int_t AliAODTrack::GetTOFBunchCrossing(Double_t b, Bool_t) const
   return bcid;
 }
 
+void AliAODTrack::SetDetectorPID(const AliDetectorPID *pid)
+{
+  //
+  // Set the detector PID
+  //
+  if (fDetectorPID) delete fDetectorPID;
+  fDetectorPID=pid;
+  
+}
+
 //_____________________________________________________________________________
 Bool_t AliAODTrack::GetXYZAt(Double_t x, Double_t b, Double_t *r) const
 {
index 7d4e44c..1a5487e 100644 (file)
@@ -20,6 +20,7 @@
  
 
 class AliVVertex;
+class AliDetectorPID;
 class AliTPCdEdxInfo;
 class AliAODEvent;
 
@@ -368,7 +369,9 @@ class AliAODTrack : public AliVTrack {
   void     SetProdVertex(TObject *vertex) { fProdVertex = vertex; }
   void     SetType(AODTrk_t ttype) { fType=ttype; }
 
-
+  // Trasient PID object, is owned by the track
+  virtual void  SetDetectorPID(const AliDetectorPID *pid);
+  virtual const AliDetectorPID* GetDetectorPID() const { return fDetectorPID; }
 
   // Dummy
   Int_t    PdgCode() const {return 0;}
@@ -413,7 +416,8 @@ class AliAODTrack : public AliVTrack {
 
   
   AliAODRedCov<6> *fCovMatrix;      // covariance matrix (x, y, z, px, py, pz)
-  AliAODPid    *fDetPid;            // more detailed or detector specific pid information
+  AliAODPid    *fDetPid;            // more detailed or detector specific raw pid information
+  mutable const AliDetectorPID* fDetectorPID; //! transient object to cache calibrated PID information
   TRef          fProdVertex;        // vertex of origin
 
   Double_t      fTrackPhiOnEMCal;   // phi of track after being propagated to 430cm
index 34c022f..d3c1d40 100644 (file)
@@ -34,6 +34,8 @@
 #include "AliAODMCHeader.h"
 #include "AliAODMCParticle.h"
 
+#include <AliDetectorPID.h>
+
 ClassImp(AliAODpidUtil)
 
   Int_t AliAODpidUtil::MakePID(const AliAODTrack *track,Double_t *p) const {
@@ -239,8 +241,8 @@ void AliAODpidUtil::MakeTOFPID(const AliAODTrack *track, Double_t *p) const
   if ((track->GetStatus()&AliESDtrack::kTIME )==0)     return;
   if ((track->GetStatus()&AliESDtrack::kTOFpid )==0)   return;
 
-  Double_t time[AliPID::kSPECIESN];
-  Double_t sigma[AliPID::kSPECIESN];
+  Double_t time[AliPID::kSPECIES];
+  Double_t sigma[AliPID::kSPECIES];
   AliAODPid *pidObj = track->GetDetPid();
   pidObj->GetIntegratedTimes(time);
 
@@ -296,3 +298,45 @@ void AliAODpidUtil::MakeTRDPID(const AliAODTrack *track,Double_t *p) const
   return;
 }
 
+//_________________________________________________________________________
+Float_t AliAODpidUtil::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Number of sigma implementation for the TOF
+  //
+  
+  AliAODTrack *track=(AliAODTrack*)vtrack;
+
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetNumberOfSigmas(kTOF, type);
+  }  
+  
+  Bool_t oldAod=kTRUE;
+  Double_t sigTOF;
+  AliAODPid *pidObj = track->GetDetPid();
+  if (!pidObj) return -999.;
+  Double_t tofTime=pidObj->GetTOFsignal();
+  Double_t expTime=fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type);
+  Double_t sigmaTOFPid[AliPID::kSPECIES];
+  pidObj->GetTOFpidResolution(sigmaTOFPid);
+  AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
+  if (event) {  // protection if the user didn't call GetTrack, which sets the internal pointer
+    AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
+    if (tofH && (TMath::Abs(sigmaTOFPid[0]) <= 1.E-16) ) { // new AOD
+        sigTOF=fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type)); //fTOFResponse is set in InitialiseEvent
+        tofTime -= fTOFResponse.GetStartTime(vtrack->P());
+        oldAod=kFALSE;
+    }
+  } else {
+    AliError("pointer to AliAODEvent not found, please call GetTrack to set it");
+    return -996.;
+  }
+  if (oldAod) { // old AOD
+    if (type <= AliPID::kProton) {
+      sigTOF=sigmaTOFPid[type];
+    } else return -998.;  // light nuclei cannot be supported on old AOD because we don't have timeZero resolution
+  }
+  if (sigTOF>0) return (tofTime - expTime)/sigTOF;
+  else return -997.;
+}
index 7f3718a..c962d9b 100644 (file)
@@ -48,36 +48,6 @@ private:
   ClassDef(AliAODpidUtil,3)  // PID calculation class
 };
 
-inline Float_t AliAODpidUtil::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const {
-  AliAODTrack *track=(AliAODTrack*)vtrack;
-  Bool_t oldAod=kTRUE;
-  Double_t sigTOF;
-  AliAODPid *pidObj = track->GetDetPid();
-  if (!pidObj) return -999.;
-  Double_t tofTime=pidObj->GetTOFsignal();
-  Double_t expTime=fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type);
-  Double_t sigmaTOFPid[AliPID::kSPECIES];
-  pidObj->GetTOFpidResolution(sigmaTOFPid);
-  AliAODEvent *event=(AliAODEvent*)track->GetAODEvent();
-  if (event) {  // protection if the user didn't call GetTrack, which sets the internal pointer
-    AliTOFHeader* tofH=(AliTOFHeader*)event->GetTOFHeader();
-    if (tofH && (TMath::Abs(sigmaTOFPid[0]) <= 1.E-16) ) { // new AOD
-       sigTOF=fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type)); //fTOFResponse is set in InitialiseEvent
-       tofTime -= fTOFResponse.GetStartTime(vtrack->P());
-       oldAod=kFALSE;
-    } 
-  } else {
-    AliError("pointer to AliAODEvent not found, please call GetTrack to set it");
-    return -996.;
-  }
-  if (oldAod) { // old AOD
-    if (type <= AliPID::kProton) {
-      sigTOF=sigmaTOFPid[type];
-    } else return -998.;  // light nuclei cannot be supported on old AOD because we don't have timeZero resolution
-  }
-  if (sigTOF>0) return (tofTime - expTime)/sigTOF;
-  else return -997.;
-}
 
 #endif
 
index ef4d53d..c3187b3 100644 (file)
@@ -72,6 +72,8 @@ set ( SRCS
     STEERBase/AliTRDPIDResponse.cxx 
     STEERBase/AliEMCALPIDResponse.cxx 
     STEERBase/AliPIDCombined.cxx
+    STEERBase/AliPIDValues.cxx
+    STEERBase/AliDetectorPID.cxx
     STEERBase/AliDAQ.cxx 
     STEERBase/AliRefArray.cxx
     STEERBase/AliOADBContainer.cxx     
index 4aaa0bc..e05f7d4 100644 (file)
@@ -53,7 +53,7 @@ AliESDCaloCluster::AliESDCaloCluster() :
   // The default ESD constructor 
   //
   fGlobalPos[0] = fGlobalPos[1] = fGlobalPos[2] = 0.;
-  for(Int_t i=0; i<AliPID::kSPECIESN; i++) fPID[i] = 0.;
+  for(Int_t i=0; i<AliPID::kSPECIESCN; i++) fPID[i] = 0.;
 }
 
 //_______________________________________________________________________
@@ -85,7 +85,7 @@ AliESDCaloCluster::AliESDCaloCluster(const AliESDCaloCluster& clus) :
   fGlobalPos[1] = clus.fGlobalPos[1];
   fGlobalPos[2] = clus.fGlobalPos[2];
 
-  for(Int_t i=0; i<AliPID::kSPECIESN; i++) fPID[i] = clus.fPID[i];
+  for(Int_t i=0; i<AliPID::kSPECIESCN; i++) fPID[i] = clus.fPID[i];
 
   if (clus.fNCells > 0) {
 
@@ -125,7 +125,7 @@ AliESDCaloCluster &AliESDCaloCluster::operator=(const AliESDCaloCluster& source)
   fTrackDx= source.fTrackDx ;
   fTrackDz= source.fTrackDz ;
   fDistToBadChannel = source.fDistToBadChannel ;
-  for(Int_t i=0; i<AliPID::kSPECIESN; i++) fPID[i] = source.fPID[i];
+  for(Int_t i=0; i<AliPID::kSPECIESCN; i++) fPID[i] = source.fPID[i];
   fID = source.fID;
 
   fNCells= source.fNCells;
@@ -233,7 +233,7 @@ void AliESDCaloCluster::SetPID(const Float_t *p) {
   // In case all the weights are non-positive they are replaced by
   // uniform probabilities
 
-  Int_t n = AliPID::kSPECIESN;
+  Int_t n = AliPID::kSPECIESCN;
 
   Float_t uniform = 1./(Float_t)n;
 
index f7d21b1..126df93 100644 (file)
@@ -61,9 +61,9 @@ class AliESDCaloCluster : public AliVCluster
   Double_t Chi2() const         { return fChi2; }
   
   const Double_t *GetPID() const { return fPID; }
-  //for(Int_t i=0; i<AliPID::kSPECIESN; ++i) pid[i]=fPID[i];}
+  //for(Int_t i=0; i<AliPID::kSPECIESCN; ++i) pid[i]=fPID[i];}
   void SetPID  (const Float_t *pid) ;
-  void SetPIDAt(Float_t p, Int_t i) {if(i>=0 && i<AliPID::kSPECIESN) fPID[i] = p ; 
+  void SetPIDAt(Float_t p, Int_t i) {if(i>=0 && i<AliPID::kSPECIESCN) fPID[i] = p ; 
     else AliInfo(Form("Bad index for PID array, i = %d \n",i));}
   
   void     SetM20(Double_t m20) { fM20 = m20; }
@@ -156,7 +156,7 @@ class AliESDCaloCluster : public AliVCluster
   Double32_t   fTrackDz ;         // Distance to closest track in z
   
   Double32_t   fDistToBadChannel; // Distance to nearest bad channel
-  Double32_t   fPID[AliPID::kSPECIESN]; //[0,1,8]"detector response  probabilities" (for the PID)
+  Double32_t   fPID[AliPID::kSPECIESCN]; //[0,1,8]"detector response  probabilities" (for the PID)
   Int_t        fID;                // Unique Id of the cluster
   UChar_t      fNExMax ;           // number of (Ex-)maxima before unfolding  
   Char_t       fClusterType;       // Flag for different cluster type/versions
index 79ceba1..4c027b8 100644 (file)
@@ -35,6 +35,7 @@
 #include "AliESDtrack.h"
 #include "AliMCEvent.h"
 
+#include <AliDetectorPID.h>
 
 ClassImp(AliESDpid)
 
@@ -255,7 +256,7 @@ void AliESDpid::MakeTOFPID(AliESDtrack *track, Float_t /*timeZeroTOF*/) const
   Int_t ibin = fTOFResponse.GetMomBin(track->GetP());
   Float_t timezero = fTOFResponse.GetT0bin(ibin);
 
-  Double_t time[AliPID::kSPECIESN];
+  Double_t time[AliPID::kSPECIES];
   track->GetIntegratedTimes(time);
 
   Double_t sigma[AliPID::kSPECIES];
@@ -406,3 +407,19 @@ Bool_t AliESDpid::CheckTOFMatching(AliESDtrack *track) const{
     
     return status;
 }
+
+Float_t AliESDpid::NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type, const Float_t /*timeZeroTOF*/) const
+{
+  //
+  // Number of sigma implementation for the TOF
+  //
+  
+  AliVTrack *vtrack=(AliVTrack*)track;
+  // look for cached value first
+  if (vtrack->GetDetectorPID()){
+    return vtrack->GetDetectorPID()->GetNumberOfSigmas(kTOF, type);
+  }
+  
+  Double_t expTime = fTOFResponse.GetExpectedSignal(vtrack,type);
+  return (vtrack->GetTOFsignal() - fTOFResponse.GetStartTime(vtrack->P()) - expTime)/fTOFResponse.GetExpectedSigma(vtrack->P(),expTime,AliPID::ParticleMassZ(type));
+}
index 4d272fd..4d9a3a4 100644 (file)
@@ -40,8 +40,8 @@ AliESDpid(const AliESDpid&a): AliPIDResponse(a), fRangeTOFMismatch(a.fRangeTOFMi
   void MakeTRDPID(AliESDtrack *track) const;
   void CombinePID(AliESDtrack *track) const;
   
-  virtual Float_t NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type, const Float_t timeZeroTOF) const;
-  virtual Float_t NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type) const {return NumberOfSigmasTOF(vtrack,type,0); }
+  virtual Float_t NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type, const Float_t timeZeroTOF) const;
+  virtual Float_t NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type) const {return NumberOfSigmasTOF(track,type,0); }
   
   void SetNMaxSigmaTOFTPCMismatch(Float_t range) {fRangeTOFMismatch=range;}
   Float_t GetNMaxSigmaTOFTPCMismatch() const {return fRangeTOFMismatch;}
@@ -59,15 +59,6 @@ private:
 };
 
 
-inline Float_t AliESDpid::NumberOfSigmasTOF(const AliVParticle *vtrack, AliPID::EParticleType type, const Float_t /*timeZeroTOF*/) const {
-  AliESDtrack *track=(AliESDtrack*)vtrack;
-  //  Double_t times[AliPID::kSPECIES];
-  // track->GetIntegratedTimes(times);
-  Double_t expTime = fTOFResponse.GetExpectedSignal((AliVTrack*)vtrack,type);
-  //  return (track->GetTOFsignal() - fTOFResponse.GetStartTime(track->GetP()) - times[type])/fTOFResponse.GetExpectedSigma(track->GetP(),times[type],AliPID::ParticleMass(type));
-  return (track->GetTOFsignal() - fTOFResponse.GetStartTime(track->GetP()) - expTime)/fTOFResponse.GetExpectedSigma(track->GetP(),expTime,AliPID::ParticleMassZ(type));
-}
-
 #endif
 
 
index 0a148c2..68a77bd 100644 (file)
 #include "TPolyMarker3D.h"
 #include "AliTrackerBase.h"
 #include "AliTPCdEdxInfo.h"
+#include "AliDetectorPID.h"
 
 ClassImp(AliESDtrack)
 
@@ -239,6 +240,7 @@ AliESDtrack::AliESDtrack() :
   fCacheNCrossedRows(-10),
   fCacheChi2TPCConstrainedVsGlobal(-10),
   fCacheChi2TPCConstrainedVsGlobalVertex(0),
+  fDetectorPID(0x0),
   fTrackPhiOnEMCal(-999),
   fTrackEtaOnEMCal(-999)
 {
@@ -354,6 +356,7 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
   fCacheNCrossedRows(track.fCacheNCrossedRows),
   fCacheChi2TPCConstrainedVsGlobal(track.fCacheChi2TPCConstrainedVsGlobal),
   fCacheChi2TPCConstrainedVsGlobalVertex(track.fCacheChi2TPCConstrainedVsGlobalVertex),
+  fDetectorPID(0x0),
   fTrackPhiOnEMCal(track.fTrackPhiOnEMCal),
   fTrackEtaOnEMCal(track.fTrackEtaOnEMCal)
 {
@@ -381,6 +384,7 @@ AliESDtrack::AliESDtrack(const AliESDtrack& track):
     for (Int_t i=0; i<fTRDnSlices; i++) fTRDslices[i]=track.fTRDslices[i];
   }
 
+  if (track.fDetectorPID) fDetectorPID = new AliDetectorPID(*track.fDetectorPID);
 
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTRDr[i]=track.fTRDr[i]; 
   for (Int_t i=0;i<AliPID::kSPECIES;i++) fTOFr[i]=track.fTOFr[i];
@@ -481,6 +485,7 @@ AliESDtrack::AliESDtrack(const AliVTrack *track) :
   fCacheNCrossedRows(-10),
   fCacheChi2TPCConstrainedVsGlobal(-10),
   fCacheChi2TPCConstrainedVsGlobalVertex(0),
+  fDetectorPID(0x0),
   fTrackPhiOnEMCal(-999),
   fTrackEtaOnEMCal(-999)
 {
@@ -628,6 +633,7 @@ AliESDtrack::AliESDtrack(TParticle * part) :
   fCacheNCrossedRows(-10),
   fCacheChi2TPCConstrainedVsGlobal(-10),
   fCacheChi2TPCConstrainedVsGlobalVertex(0),
+  fDetectorPID(0x0),
   fTrackPhiOnEMCal(-999),
   fTrackEtaOnEMCal(-999)
 {
@@ -769,6 +775,7 @@ AliESDtrack::~AliESDtrack(){
   fCacheChi2TPCConstrainedVsGlobal = -10.;
   if(fCacheChi2TPCConstrainedVsGlobalVertex) fCacheChi2TPCConstrainedVsGlobalVertex = 0;
 
+  delete fDetectorPID;
 }
 
 AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
@@ -979,6 +986,10 @@ AliESDtrack &AliESDtrack::operator=(const AliESDtrack &source){
   fCacheChi2TPCConstrainedVsGlobal = source.fCacheChi2TPCConstrainedVsGlobal;
   fCacheChi2TPCConstrainedVsGlobalVertex = source.fCacheChi2TPCConstrainedVsGlobalVertex;
 
+  delete fDetectorPID;
+  fDetectorPID=0x0;
+  if (source.fDetectorPID) fDetectorPID = new AliDetectorPID(*source.fDetectorPID);
+  
   fTrackPhiOnEMCal= source.fTrackPhiOnEMCal;
   fTrackEtaOnEMCal= source.fTrackEtaOnEMCal;
 
@@ -2414,7 +2425,7 @@ void AliESDtrack::Print(Option_t *) const {
   // Prints info on the track
   AliExternalTrackParam::Print();
   printf("ESD track info\n") ; 
-  Double_t p[AliPID::kSPECIESN] ; 
+  Double_t p[AliPID::kSPECIES] ;
   Int_t index = 0 ; 
   if( IsOn(kITSpid) ){
     printf("From ITS: ") ; 
@@ -2632,3 +2643,13 @@ Double_t AliESDtrack::GetChi2TPCConstrainedVsGlobal(const AliESDVertex* vtx) con
   fCacheChi2TPCConstrainedVsGlobal = chi2(0,0);
   return fCacheChi2TPCConstrainedVsGlobal;
 }
+
+void AliESDtrack::SetDetectorPID(const AliDetectorPID *pid)
+{
+  //
+  // Set the detector PID
+  //
+  if (fDetectorPID) delete fDetectorPID;
+  fDetectorPID=pid;
+  
+}
index 83218a6..3f40cb0 100644 (file)
@@ -40,7 +40,7 @@ class AliESDEvent;
 class AliKalmanTrack;
 class AliTrackPointArray;
 class TPolyMarker3D;
-
+class AliDetectorPID;
 
 class AliESDtrack : public AliExternalTrackParam {
 public:
@@ -395,7 +395,12 @@ public:
   }
   virtual void Print(Option_t * opt) const ;
   const AliESDEvent* GetESDEvent() const {return fESDEvent;}
-  void         SetESDEvent(const AliESDEvent* evt) {fESDEvent = evt;}  
+  void         SetESDEvent(const AliESDEvent* evt) {fESDEvent = evt;}
+
+  // Trasient PID object, is owned by the track
+  virtual void  SetDetectorPID(const AliDetectorPID *pid);
+  virtual const AliDetectorPID* GetDetectorPID() const { return fDetectorPID; }
+  
   //
   // visualization (M. Ivanov)
   //
@@ -534,6 +539,8 @@ protected:
   mutable Float_t fCacheChi2TPCConstrainedVsGlobal; //! Cache for the chi2 of constrained TPC vs global track
   mutable const AliESDVertex* fCacheChi2TPCConstrainedVsGlobalVertex; //! Vertex for which the cache is valid
 
+  mutable const AliDetectorPID* fDetectorPID; //! transient object to cache PID information
+
   Double_t fTrackPhiOnEMCal;   // phi of track after being propagated to 430cm
   Double_t fTrackEtaOnEMCal;   // eta of track after being propagated to 430cm
   
index 4eb9b86..7efcb58 100644 (file)
@@ -45,7 +45,7 @@
 //             -->  Probability from Gaussian distribution              //
 //                  (proper normalization to each other?)               //
 //                                                                      //
-// NO Parametrization (outside pT range): --> return 999                //
+// NO Parametrization (outside pT range): --> return kFALSE             //
 //////////////////////////////////////////////////////////////////////////
 
 #include <TF1.h>
@@ -179,38 +179,37 @@ Double_t  AliEMCALPIDResponse::GetNumberOfSigmas( Float_t pt,  Float_t eop, AliP
 
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
+Bool_t AliEMCALPIDResponse::ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const {
   //
   //
   Double_t fRange  = 5.0;   // hardcoded (???)
   Double_t nsigma  = 0.0;
-  Double_t proba   = 999.;
   
 
   // Check the charge
   if( charge != -1 && charge != 1){
-    return proba;
+    return kFALSE;
   }
 
  
   // default value (will be returned, if pt below threshold)
-  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
-    pEMCAL[species] = 999.;
+  for (Int_t species = 0; species < nSpecies; species++) {
+    pEMCAL[species] = 1./nSpecies;
   }
 
   // set E/p range
   if(eop < 0.05) eop = 0.05;
   if(eop > 2.00) eop = 2.00;
   
-  for (Int_t species = 0; species < AliPID::kSPECIES; species++) {
+  for (Int_t species = 0; species < nSpecies; species++) {
     
     AliPID::EParticleType type = AliPID::EParticleType(species);
 
     // Get the parameters for this particle type and pt
     const TVectorD *params = GetParams(species, pt, charge);
     
-    // IF not in momentum range, NULL is returned --> return default value
-    if(!params) return proba;
+    // IF not in momentum/species (only for kSPECIES so far) range, NULL is returned --> return kFALSE
+    if(!params) return kFALSE;
 
     Double_t sigma    = (*params)[3];   // sigma value of Gausiian parametrization
     Double_t probLow  = (*params)[6];   // probability to be below eopMin
@@ -246,13 +245,9 @@ Double_t AliEMCALPIDResponse::ComputeEMCALProbability(Float_t pt, Float_t eop, I
        pEMCAL[species]*=GetExpectedNorm(pt,type,charge);
       }
     }
-    
-    // return the electron probability
-    proba = pEMCAL[AliPID::kElectron];  
-
   }
 
-  return proba;
+  return kTRUE;
 
 }
 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
index 52cf1b8..57c17fb 100644 (file)
@@ -37,8 +37,8 @@ public :
     void   SetPIDParams(const TObjArray * params) { fkPIDParams = params; }
     
 
-    // EMCAL probability -> should go to another place?
-    Double_t ComputeEMCALProbability( Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const;
+    // EMCAL probability
+    Bool_t ComputeEMCALProbability(Int_t nSpecies, Float_t pt, Float_t eop, Int_t charge, Double_t *pEMCAL) const;
 
 protected:
   
index c4b63d6..3556220 100644 (file)
@@ -36,7 +36,8 @@ public:
  Double_t GetResolution(Double_t bethe, Int_t nPtsForPid=4, Bool_t isSA=kFALSE) const;
  void GetITSProbabilities(Float_t mom, Double_t qclu[4], Double_t condprobfun[AliPID::kSPECIES],Bool_t isMC=kFALSE) const;
  Float_t GetNumberOfSigmas(Float_t mom, Float_t signal, AliPID::EParticleType type, Int_t nPtsForPid=4, Bool_t isSA=kFALSE) const {
-   Float_t bethe = Bethe(mom,AliPID::ParticleMass(type),isSA);
+   const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(type),2.);
+   Float_t bethe = Bethe(mom,AliPID::ParticleMassZ(type),isSA)*chargeFactor;
    return (signal - bethe)/GetResolution(bethe,nPtsForPid,isSA);
  }
  Int_t GetParticleIdFromdEdxVsP(Float_t mom, Float_t signal, Bool_t isSA=kFALSE) const;
index 7c5b261..96d2030 100644 (file)
 
 ClassImp(AliPID)
 
-const char* AliPID::fgkParticleName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
+const char* AliPID::fgkParticleName[AliPID::kSPECIESCN+1] = {
   "electron",
   "muon",
   "pion",
-  "kaon",
+  "kaon",  
   "proton",
+  
+  "deuteron",
+  "triton",
+  "helium-3",
+  "alpha",
+  
   "photon",
   "pi0",
   "neutron",
   "kaon0",
   "eleCon",
-  "deuteron",
-  "triton",
-  "helium-3",
-  "alpha",
+  
   "unknown"
 };
 
-const char* AliPID::fgkParticleShortName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
+const char* AliPID::fgkParticleShortName[AliPID::kSPECIESCN+1] = {
   "e",
   "mu",
   "pi",
   "K",
   "p",
+
+  "d",
+  "t",
+  "he3",
+  "alpha",
+
   "photon",
   "pi0",
   "n",
   "K0",
   "eleCon",
-  "d",
-  "t",
-  "he3",
-  "alpha",
+  
   "unknown"
 };
 
-const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
+const char* AliPID::fgkParticleLatexName[AliPID::kSPECIESCN+1] = {
   "e",
   "#mu",
   "#pi",
   "K",
   "p",
+
+  "d",
+  "t",
+  "he3",
+  "#alpha",
+
   "#gamma",
   "#pi_{0}",
   "n",
   "K_{0}",
   "eleCon",
-  "d",
-  "t",
-  "he3",
-  "#alpha",
+
   "unknown"
 };
 
-const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
+const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESCN+1] = {
   ::kElectron, 
   ::kMuonMinus, 
   ::kPiPlus, 
   ::kKPlus, 
   ::kProton,
+  
+  1000010020,
+  1000010030,
+  1000020030,
+  1000020040,
+
   ::kGamma,
   ::kPi0,
   ::kNeutron,
   ::kK0,
   ::kElectron,
-  1000010020,
-  1000010030,
-  1000020030,
-  1000020040,
   0
 };
 
-/*const*/ Float_t AliPID::fgkParticleMass[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
+/*const*/ Float_t AliPID::fgkParticleMass[AliPID::kSPECIESCN+1] = {
   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
   /*
   M(kElectron),  // electron
@@ -142,7 +153,7 @@ const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
   0.00000        // unknown
   */
 };
-/*const*/ Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
+/*const*/ Float_t AliPID::fgkParticleMassZ[AliPID::kSPECIESCN+1] = {
   0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
   /*
   M(kElectron),  // electron
@@ -163,7 +174,7 @@ const Int_t AliPID::fgkParticleCode[AliPID::kSPECIESN+AliPID::kSPECIESLN+1] = {
   */
 };
 
-Double_t AliPID::fgPrior[kSPECIESN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+Double_t AliPID::fgPrior[kSPECIESCN] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
 
 
 //_______________________________________________________________________
@@ -176,8 +187,8 @@ AliPID::AliPID() :
   //
   Init();
   // set default values (= equal probabilities)
-  for (Int_t i = 0; i < kSPECIESN; i++)
-    fProbDensity[i] = 1./kSPECIESN;
+  for (Int_t i = 0; i < kSPECIESCN; i++)
+    fProbDensity[i] = 1./kSPECIESCN;
 }
 
 //_______________________________________________________________________
@@ -190,10 +201,10 @@ AliPID::AliPID(const Double_t* probDensity, Bool_t charged) :
   //
   Init();
   // set given probability densities
-  for (Int_t i = 0; i < kSPECIES; i++) 
+  for (Int_t i = 0; i < kSPECIESC; i++)
     fProbDensity[i] = probDensity[i];
 
-  for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
+  for (Int_t i = kSPECIESC; i < kSPECIESCN; i++)
     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 }
 
@@ -207,10 +218,10 @@ AliPID::AliPID(const Float_t* probDensity, Bool_t charged) :
   //
   Init();
   // set given probability densities
-  for (Int_t i = 0; i < kSPECIES; i++) 
+  for (Int_t i = 0; i < kSPECIESC; i++) 
     fProbDensity[i] = probDensity[i];
 
-  for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
+  for (Int_t i = kSPECIESC; i < kSPECIESCN; i++) 
     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 }
 
@@ -223,7 +234,7 @@ AliPID::AliPID(const AliPID& pid) :
   // copy constructor
   //
   // We do not call init here, MUST already be done
-  for (Int_t i = 0; i < kSPECIESN; i++) 
+  for (Int_t i = 0; i < kSPECIESCN; i++) 
     fProbDensity[i] = pid.fProbDensity[i];
 }
 
@@ -233,10 +244,10 @@ void AliPID::SetProbabilities(const Double_t* probDensity, Bool_t charged)
   //
   // Set the probability densities
   //
-  for (Int_t i = 0; i < kSPECIES; i++) 
+  for (Int_t i = 0; i < kSPECIESC; i++) 
     fProbDensity[i] = probDensity[i];
 
-  for (Int_t i = kSPECIES; i < kSPECIESN; i++) 
+  for (Int_t i = kSPECIESC; i < kSPECIESCN; i++) 
     fProbDensity[i] = ((charged) ? 0 : probDensity[i]);
 }
 
@@ -247,7 +258,7 @@ AliPID& AliPID::operator = (const AliPID& pid)
 
   if(this != &pid) {
     fCharged = pid.fCharged;
-    for (Int_t i = 0; i < kSPECIESN; i++) {
+    for (Int_t i = 0; i < kSPECIESCN; i++) {
       fProbDensity[i] = pid.fProbDensity[i];
     }
   }
@@ -263,7 +274,7 @@ void AliPID::Init()
   // Initialise only once... 
   if(!fgkParticleMass[0]) {
     AliPDG::AddParticlesToPdgDataBase();
-    for (Int_t i = 0; i < kSPECIESN + kSPECIESLN; i++) {
+    for (Int_t i = 0; i < kSPECIESC; i++) {
       fgkParticleMass[i] = M(i);
       if (i == kHe3 || i == kAlpha) fgkParticleMassZ[i] = M(i)/2.;
       else fgkParticleMassZ[i]=M(i);
@@ -280,7 +291,7 @@ Double_t AliPID::GetProbability(EParticleType iType,
   // assuming the a priori probabilities "prior"
   //
   Double_t sum = 0.;
-  Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
+  Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
   for (Int_t i = 0; i < nSpecies; i++) {
     sum += fProbDensity[i] * prior[i];
   }
@@ -308,7 +319,7 @@ void AliPID::GetProbabilities(Double_t* probabilities,
 // assuming the a priori probabilities "prior"
 
   Double_t sum = 0.;
-  Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
+  Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
   for (Int_t i = 0; i < nSpecies; i++) {
     sum += fProbDensity[i] * prior[i];
   }
@@ -339,7 +350,7 @@ AliPID::EParticleType AliPID::GetMostProbable(const Double_t* prior) const
 
   Double_t max = 0.;
   EParticleType id = kPion;
-  Int_t nSpecies = ((fCharged) ? kSPECIES : kSPECIESN);
+  Int_t nSpecies = ((fCharged) ? kSPECIESC : kSPECIESCN);
   for (Int_t i = 0; i < nSpecies; i++) {
     Double_t prob = fProbDensity[i] * prior[i];
     if (prob > max) {
@@ -369,8 +380,8 @@ void AliPID::SetPriors(const Double_t* prior, Bool_t charged)
 // use the given priors as global a priori probabilities
 
   Double_t sum = 0;
-  for (Int_t i = 0; i < kSPECIESN; i++) {
-    if (charged && (i >= kSPECIES)) {
+  for (Int_t i = 0; i < kSPECIESCN; i++) {
+    if (charged && (i >= kSPECIESC)) {
       fgPrior[i] = 0;      
     } else {
       if (prior[i] < 0) {
@@ -409,7 +420,7 @@ AliPID& AliPID::operator *= (const AliPID& pid)
 {
 // combine this probability densities with the one of "pid"
 
-  for (Int_t i = 0; i < kSPECIESN; i++) {
+  for (Int_t i = 0; i < kSPECIESCN; i++) {
     fProbDensity[i] *= pid.fProbDensity[i];
   }
   return *this;
index e254412..c4b6875 100644 (file)
 
 
 #include <TObject.h>
-
+#include <TMath.h>
 
 class AliPID : public TObject {
  public:
   enum {
-    kSPECIES = 5,    // Number of particle species recognized by the PID
-    kSPECIESN = 10,  // Number of charged+neutral particle species recognized by the PHOS/EMCAL PID
-    kSPECIESLN = 4   // Number of light nuclei: deuteron, triton, helium-3 and alpha
+    kSPECIES = 5,     // Number of default particle species recognized by the PID
+    kSPECIESC = 9,    // Number of default particles + light nuclei recognized by the PID
+    kSPECIESCN = 14,  // Number of charged+neutral particle species recognized by the PHOS/EMCAL PID
   };
   enum EParticleType {
     kElectron = 0, 
@@ -29,17 +29,21 @@ class AliPID : public TObject {
     kPion = 2, 
     kKaon = 3, 
     kProton = 4, 
-    kPhoton = 5, 
-    kPi0 = 6, 
-    kNeutron = 7, 
-    kKaon0 = 8, 
-    kEleCon = 9,
-    kDeuteron = 10,
-    kTriton = 11,
-    kHe3 = 12,
-    kAlpha = 13,
+
+    kDeuteron = 5,
+    kTriton = 6,
+    kHe3 = 7,
+    kAlpha = 8,
+    
+    kPhoton = 9,
+    kPi0 = 10, 
+    kNeutron = 11, 
+    kKaon0 = 12, 
+    kEleCon = 13,
+    
     kUnknown = 14
   };
+  
   static Float_t       ParticleMass(Int_t iType) {
      if(!fgkParticleMass[0]) Init(); 
      return fgkParticleMass[iType];
@@ -48,6 +52,11 @@ class AliPID : public TObject {
      if(!fgkParticleMass[0]) Init(); 
      return fgkParticleMassZ[iType];
   }
+  static Int_t         ParticleCharge(Int_t iType){
+    if(!fgkParticleMass[0]) Init();
+    if (iType<0||iType>=kSPECIESC) return 0;
+    return TMath::Nint(fgkParticleMass[iType]/fgkParticleMassZ[iType]);
+  }
   static const char*   ParticleName(Int_t iType) 
     {return fgkParticleName[iType];};
   static const char*   ParticleShortName(Int_t iType) 
@@ -85,18 +94,18 @@ class AliPID : public TObject {
 
   static void          Init();
 
-  Bool_t               fCharged;                   // flag for charged/neutral
-  Double_t             fProbDensity[kSPECIESN];    // probability densities
-  static Double_t      fgPrior[kSPECIESN];         // a priori probabilities
+  Bool_t               fCharged;                           // flag for charged/neutral
+  Double_t             fProbDensity[kSPECIESCN];           // probability densities
+  static Double_t      fgPrior[kSPECIESCN];                // a priori probabilities
 
-  static /*const*/ Float_t fgkParticleMass[kSPECIESN+kSPECIESLN+1]; // particle masses
-  static /*const*/ Float_t fgkParticleMassZ[kSPECIESN+kSPECIESLN+1]; // particle masses/charge
-  static const char*   fgkParticleName[kSPECIESN+kSPECIESLN+1]; // particle names
-  static const char*   fgkParticleShortName[kSPECIESN+kSPECIESLN+1]; // particle names
-  static const char*   fgkParticleLatexName[kSPECIESN+kSPECIESLN+1]; // particle names
-  static const Int_t   fgkParticleCode[kSPECIESN+kSPECIESLN+1]; // particle codes
+  static /*const*/ Float_t fgkParticleMass[kSPECIESCN+1];  // particle masses
+  static /*const*/ Float_t fgkParticleMassZ[kSPECIESCN+1]; // particle masses/charge
+  static const char*   fgkParticleName[kSPECIESCN+1];      // particle names
+  static const char*   fgkParticleShortName[kSPECIESCN+1]; // particle names
+  static const char*   fgkParticleLatexName[kSPECIESCN+1]; // particle names
+  static const Int_t   fgkParticleCode[kSPECIESCN+1];      // particle codes
 
-  ClassDef(AliPID, 3)    // particle id probability densities
+  ClassDef(AliPID, 3)                                      // particle id probability densities
 };
 
 
index 15a6556..248a7f8 100644 (file)
@@ -38,7 +38,7 @@
 #include "AliOADBContainer.h"
 
 // initialize static members
-TH2F* AliPIDCombined::fDefaultPriorsTPC[5];
+TH2F* AliPIDCombined::fDefaultPriorsTPC[]={0x0};
 
 ClassImp(AliPIDCombined);
 
@@ -53,7 +53,7 @@ AliPIDCombined::AliPIDCombined():
   //
   // default constructor
   //   
-  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
+  for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
   AliLog::SetClassDebugLevel("AliPIDCombined",10);
 }
 
@@ -69,7 +69,7 @@ AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
   //
   // default constructor with name and title
   //
-  for (Int_t i=0;i<AliPID::kSPECIES+AliPID::kSPECIESLN;i++) fPriorsDistributions[i]=NULL;
+  for (Int_t i=0;i<AliPID::kSPECIESC;i++) fPriorsDistributions[i]=NULL;
   AliLog::SetClassDebugLevel("AliPIDCombined",10);
 
 }
@@ -77,7 +77,7 @@ AliPIDCombined::AliPIDCombined(const TString& name,const TString& title):
 //-------------------------------------------------------------------------------------------------    
 AliPIDCombined::~AliPIDCombined() {
 
-  for(Int_t i=0;i < AliPID::kSPECIES+AliPID::kSPECIESLN;i++){
+  for(Int_t i=0;i < AliPID::kSPECIESC;i++){
     if(fPriorsDistributions[i])
       delete fPriorsDistributions[i];
   }
@@ -85,22 +85,12 @@ AliPIDCombined::~AliPIDCombined() {
 
 //-------------------------------------------------------------------------------------------------    
 void AliPIDCombined::SetPriorDistribution(AliPID::EParticleType type,TH1F *prior) {
-  if ( (type < 0) || ( type >= (AliPID::kSPECIESN+AliPID::kSPECIESLN) ) ){
+  if ( (type < 0) || ( type >= ((AliPID::EParticleType)AliPID::kSPECIESC) ) ){
     AliError(Form("Invalid EParticleType setting prior (offending type: %d)",type));
     return;
   }
   if(prior) {
     Int_t i = (Int_t)type;
-    if (i >= AliPID::kSPECIES ) {
-      if (i < (Int_t)AliPID::kDeuteron) {
-       AliError(Form("Invalid EParticleType setting prior. Type: %d (neutral) not supported",i));
-       return;
-      } else i -= (AliPID::kSPECIESN-AliPID::kSPECIES);
-    }
-    if (i>(AliPID::kSPECIES+AliPID::kSPECIESLN)) {             // coverity fix (to put it mildly....)
-      AliError(Form("Unexpected EParticleType setting prior. Type: %d (neutral) not supported",i));
-      return;
-    }
     if (fPriorsDistributions[i]) {
       delete fPriorsDistributions[i]; 
     }
@@ -371,20 +361,29 @@ void AliPIDCombined::SetDefaultTPCPriors(){
   fEnablePriors=kTRUE;
   fUseDefaultTPCPriors = kTRUE;
 
+  //check if priors are already initialized
+  if (fDefaultPriorsTPC[0]) return;
+  
   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()));
+  if(!foadb || !foadb->IsOpen()) {
+    delete foadb;
+    AliFatal(Form("Cannot open OADB file %s", oadbfilename.Data()));
+    return;
+  }
   
   AliOADBContainer * priorsContainer = (AliOADBContainer*) foadb->Get("priorsTPC");
   if (!priorsContainer) AliFatal("Cannot fetch OADB container for Priors");
   
-  const char *namespecies[5] = {"El","Mu","Pi","Ka","Pr"};
+  const char *namespecies[AliPID::kSPECIES] = {"El","Mu","Pi","Ka","Pr"};
   char name[100];
 
-  for(Int_t i=0;i < 5;i++){
+  for(Int_t i=0;i < AliPID::kSPECIES;i++){
     snprintf(name,99,"hDefault%sPriors",namespecies[i]);
     fDefaultPriorsTPC[i] = (TH2F*) priorsContainer->GetDefaultObject(name);
     if (!fDefaultPriorsTPC[i]) AliFatal(Form("Cannot find priors for %s", namespecies[i]));
   }
+
+  delete foadb;
 }
index ce1a52b..d7704b2 100644 (file)
@@ -60,11 +60,11 @@ private:
   Int_t fRejectMismatchMask; // Detectors set return flat prob. if mismatch detected \r
   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
+  TH1F *fPriorsDistributions[AliPID::kSPECIESC]; // priors\r
   Bool_t fUseDefaultTPCPriors; // switch to use Defaul TPC Priors\r
-  static TH2F *fDefaultPriorsTPC[5]; // Default priors for TPC tracks\r
+  static TH2F *fDefaultPriorsTPC[AliPID::kSPECIES]; // Default priors for TPC tracks\r
 \r
-  ClassDef(AliPIDCombined,2);\r
+  ClassDef(AliPIDCombined, 3);   // Combined PID using priors\r
 };\r
 \r
 #endif\r
index 491b6b0..08f33d2 100644 (file)
@@ -29,6 +29,7 @@
 #include <TF1.h>
 #include <TSpline.h>
 #include <TFile.h>
+#include <TArrayI.h>
 #include <TArrayF.h>
 
 #include <AliVEvent.h>
@@ -40,6 +41,7 @@
 #include <AliTOFPIDParams.h>
 
 #include "AliPIDResponse.h"
+#include "AliDetectorPID.h"
 
 #include "AliCentrality.h"
 
@@ -56,6 +58,7 @@ fRange(5.),
 fITSPIDmethod(kITSTruncMean),
 fIsMC(isMC),
 fOADBPath(),
+fCustomTPCpidResponse(),
 fBeamType("PP"),
 fLHCperiod(),
 fMCperiodTPC(),
@@ -91,9 +94,9 @@ AliPIDResponse::~AliPIDResponse()
   //
   // dtor
   //
-    delete fArrPidResponseMaster;
-    delete fTRDPIDResponseObject;
-  if (fTOFPIDParams) delete fTOFPIDParams;
+  delete fArrPidResponseMaster;
+  delete fTRDPIDResponseObject;
+  delete fTOFPIDParams;
 }
 
 //______________________________________________________________________________
@@ -108,6 +111,7 @@ fRange(other.fRange),
 fITSPIDmethod(other.fITSPIDmethod),
 fIsMC(other.fIsMC),
 fOADBPath(other.fOADBPath),
+fCustomTPCpidResponse(other.fCustomTPCpidResponse),
 fBeamType("PP"),
 fLHCperiod(),
 fMCperiodTPC(),
@@ -150,6 +154,7 @@ AliPIDResponse& AliPIDResponse::operator=(const AliPIDResponse &other)
     fRange=other.fRange;
     fITSPIDmethod=other.fITSPIDmethod;
     fOADBPath=other.fOADBPath;
+    fCustomTPCpidResponse=other.fCustomTPCpidResponse;
     fIsMC=other.fIsMC;
     fBeamType="PP";
     fLHCperiod="";
@@ -183,18 +188,94 @@ Float_t AliPIDResponse::NumberOfSigmas(EDetCode detCode, const AliVParticle *tra
     case kDetITS: return NumberOfSigmasITS(track, type); break;
     case kDetTPC: return NumberOfSigmasTPC(track, type); break;
     case kDetTOF: return NumberOfSigmasTOF(track, type); break;
-//     case kDetTRD: return ComputeTRDProbability(track, type); break;
-//     case kDetPHOS: return ComputePHOSProbability(track, type); break;
-//     case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
-//     case kDetHMPID: return ComputeHMPIDProbability(track, type); break;
+    case kDetEMCAL: return NumberOfSigmasEMCAL(track, type); break;
     default: return -999.;
   }
 
 }
 
 //______________________________________________________________________________
-Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const {
+Float_t AliPIDResponse::NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const
+{
+  //
+  // NumberOfSigmas for 'detCode'
+  //
+  return NumberOfSigmas((EDetCode)(1<<detCode), track, type);
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the ITS
+  //
+  
+  AliVTrack *track=(AliVTrack*)vtrack;
+  
+  // look for cached value first
+  // only the non SA tracks are cached
+  if ( track->GetDetectorPID() ){
+    return track->GetDetectorPID()->GetNumberOfSigmas(kITS, type);
+  }
+  
+  Float_t dEdx=track->GetITSsignal();
+  if (dEdx<=0) return -999.;
+  
+  UChar_t clumap=track->GetITSClusterMap();
+  Int_t nPointsForPid=0;
+  for(Int_t i=2; i<6; i++){
+    if(clumap&(1<<i)) ++nPointsForPid;
+  }
+  Float_t mom=track->P();
+
+  //check for ITS standalone tracks
+  Bool_t isSA=kTRUE;
+  if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
+  
+  //TODO: in case of the electron, use the SA parametrisation,
+  //      this needs to be changed if ITS provides a parametrisation
+  //      for electrons also for ITS+TPC tracks
+  return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA || (type==AliPID::kElectron));
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the TPC
+  //
+  
+  AliVTrack *track=(AliVTrack*)vtrack;
+  
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetNumberOfSigmas(kTPC, type);
+  }
+  
+  Double_t mom  = track->GetTPCmomentum();
+  Double_t sig  = track->GetTPCsignal();
+  UInt_t   sigN = track->GetTPCsignalN();
+  
+  Double_t nSigma = -999.;
+  if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
+  
+  return nSigma;
+}
+
+//______________________________________________________________________________
+Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type) const
+{
+  //
+  // Calculate the number of sigmas in the EMCAL
+  //
+  
+  AliVTrack *track=(AliVTrack*)vtrack;
 
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetNumberOfSigmas(kEMCAL, type);
+  }
+  
   AliVCluster *matchedClus = NULL;
 
   Double_t mom     = -1.; 
@@ -235,8 +316,10 @@ Float_t AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EPar
 }
 
 //______________________________________________________________________________
-Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
+Float_t  AliPIDResponse::NumberOfSigmasEMCAL(const AliVParticle *vtrack, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const {
 
+  AliVTrack *track=(AliVTrack*)vtrack;
+  
   AliVCluster *matchedClus = NULL;
 
   Double_t mom     = -1.; 
@@ -297,8 +380,8 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode d
   switch (detCode){
     case kDetITS: return ComputeITSProbability(track, nSpecies, p); break;
     case kDetTPC: return ComputeTPCProbability(track, nSpecies, p); break;
-    case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
     case kDetTRD: return ComputeTRDProbability(track, nSpecies, p); break;
+    case kDetTOF: return ComputeTOFProbability(track, nSpecies, p); break;
     case kDetPHOS: return ComputePHOSProbability(track, nSpecies, p); break;
     case kDetEMCAL: return ComputeEMCALProbability(track, nSpecies, p); break;
     case kDetHMPID: return ComputeHMPIDProbability(track, nSpecies, p); break;
@@ -307,24 +390,41 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetCode d
 }
 
 //______________________________________________________________________________
+AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePIDProbability  (EDetector detCode,  const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
+{
+  //
+  // Compute PID response of 'detCode'
+  //
+
+  return ComputePIDProbability((EDetCode)(1<<detCode),track,nSpecies,p);
+}
+
+//______________________________________________________________________________
 AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const
 {
   //
   // Compute PID response for the ITS
   //
 
+  // look for cached value first
+  // only the non SA tracks are cached
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies);
+  }
+
   // set flat distribution (no decision)
   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 
   if ((track->GetStatus()&AliVTrack::kITSin)==0 &&
     (track->GetStatus()&AliVTrack::kITSout)==0) return kDetNoSignal;
+
+  //check for ITS standalone tracks
+  Bool_t isSA=kTRUE;
+  if( track->GetStatus() & AliVTrack::kTPCin ) isSA=kFALSE;
   
   Double_t mom=track->P();
   Double_t dedx=track->GetITSsignal();
-  Bool_t isSA=kTRUE;
   Double_t momITS=mom;
-  ULong_t trStatus=track->GetStatus();
-  if(trStatus&AliVTrack::kTPCin) isSA=kFALSE;
   UChar_t clumap=track->GetITSClusterMap();
   Int_t nPointsForPid=0;
   for(Int_t i=2; i<6; i++){
@@ -338,9 +438,13 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeITSProbability  (const AliV
 
   Bool_t mismatch=kTRUE/*, heavy=kTRUE*/;
   for (Int_t j=0; j<AliPID::kSPECIES; j++) {
-    Double_t mass=AliPID::ParticleMass(j);//GeV/c^2
-    Double_t bethe=fITSResponse.Bethe(momITS,mass);
-    Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA);
+    Double_t mass=AliPID::ParticleMassZ(j);//GeV/c^2
+    const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(j),2.);
+    Double_t bethe=fITSResponse.Bethe(momITS,mass)*chargeFactor;
+    //TODO: in case of the electron, use the SA parametrisation,
+    //      this needs to be changed if ITS provides a parametrisation
+    //      for electrons also for ITS+TPC tracks
+    Double_t sigma=fITSResponse.GetResolution(bethe,nPointsForPid,isSA || (j==(Int_t)AliPID::kElectron));
     if (TMath::Abs(dedx-bethe) > fRange*sigma) {
       p[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
     } else {
@@ -367,7 +471,12 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliV
   //
   // Compute PID response for the TPC
   //
-
+  
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetRawProbability(kTPC, p, nSpecies);
+  }
+  
   // set flat distribution (no decision)
   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 
@@ -381,7 +490,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliV
 
   if(fTuneMConData) dedx = this->GetTPCsignalTunedOnData(track);
 
-  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+  for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
     AliPID::EParticleType type=AliPID::EParticleType(j);
     Double_t bethe=fTPCResponse.GetExpectedSignal(mom,type);
     Double_t sigma=fTPCResponse.GetExpectedSigma(mom,track->GetTPCsignalN(),type);
@@ -391,12 +500,6 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTPCProbability  (const AliV
       p[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
       mismatch=kFALSE;
     }
-
-    // TODO: Light nuclei, also in TPC pid response
-    
-    // Check for particles heavier than (AliPID::kSPECIES - 1)
-//     if (dedx < (bethe + fRange*sigma)) heavy=kFALSE;
-
   }
 
   if (mismatch){
@@ -412,7 +515,12 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliV
   //
   // Compute PID response for the
   //
-
+  
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetRawProbability(kTOF, p, nSpecies);
+  }
+  
   Double_t meanCorrFactor = 0.11/fTOFtail; // Correction factor on the mean because of the tail (should be ~ 0.1 with tail = 1.1)
 
   // set flat distribution (no decision)
@@ -421,20 +529,11 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliV
   if ((track->GetStatus()&AliVTrack::kTOFout)==0) return kDetNoSignal;
   if ((track->GetStatus()&AliVTrack::kTIME)==0) return kDetNoSignal;
   
-  Double_t time[AliPID::kSPECIESN];
-  track->GetIntegratedTimes(time);
-  
-  //  Double_t sigma[AliPID::kSPECIES];
-  // for (Int_t iPart = 0; iPart < AliPID::kSPECIES; iPart++) {
-  //    sigma[iPart] = fTOFResponse.GetExpectedSigma(track->P(),time[iPart],AliPID::ParticleMass(iPart));
-  //  }
-  
   Bool_t mismatch = kTRUE/*, heavy = kTRUE*/;
-  for (Int_t j=0; j<AliPID::kSPECIES; j++) {
+  for (Int_t j=0; j<AliPID::kSPECIESC; j++) {
     AliPID::EParticleType type=AliPID::EParticleType(j);
     Double_t nsigmas=NumberOfSigmasTOF(track,type) + meanCorrFactor;
 
-    //    Double_t sig = sigma[j];
     Double_t expTime = fTOFResponse.GetExpectedSignal(track,type);
     Double_t sig = fTOFResponse.GetExpectedSigma(track->P(),expTime,AliPID::ParticleMassZ(type));
     if (TMath::Abs(nsigmas) > (fRange+2)) {
@@ -449,13 +548,6 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTOFProbability  (const AliV
        p[j] = TMath::Exp(-(nsigmas - fTOFtail*0.5)*fTOFtail)/sig;
     }
 
-    /* OLD Gaussian shape
-    if (TMath::Abs(nsigmas) > (fRange+2)) {
-      p[j] = TMath::Exp(-0.5*(fRange+2)*(fRange+2))/sig;
-    } else
-      p[j] = TMath::Exp(-0.5*nsigmas*nsigmas)/sig;
-    */
-
     if (TMath::Abs(nsigmas)<5.){
       Double_t nsigmasTPC=NumberOfSigmasTPC(track,type);
       if (TMath::Abs(nsigmasTPC)<5.) mismatch=kFALSE;
@@ -476,13 +568,18 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeTRDProbability  (const AliV
   //
   // Compute PID response for the
   //
-
+  
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetRawProbability(kTRD, p, nSpecies);
+  }
+  
   // set flat distribution (no decision)
   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
   if((track->GetStatus()&AliVTrack::kTRDout)==0) return kDetNoSignal;
 
-  Float_t mom[6];
-  Double_t dedx[48];  // Allocate space for the maximum number of TRD slices
+  Float_t mom[6]={0.};
+  Double_t dedx[48]={0.};  // Allocate space for the maximum number of TRD slices
   Int_t nslices = fTRDslicesForPID[1] - fTRDslicesForPID[0] + 1;
   AliDebug(1, Form("First Slice: %d, Last Slice: %d, Number of slices: %d",  fTRDslicesForPID[0], fTRDslicesForPID[1], nslices));
   for(UInt_t ilayer = 0; ilayer < 6; ilayer++){
@@ -500,6 +597,13 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const Al
   //
   // Compute PID response for the EMCAL
   //
+  
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetRawProbability(kEMCAL, p, nSpecies);
+  }
+
+  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
 
   AliVCluster *matchedClus = NULL;
 
@@ -515,7 +619,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const Al
   nMatchClus = track->GetEMCALcluster();
 
   if(nMatchClus > -1){
-
+    
     mom    = track->P();
     pt     = track->Pt();
     charge = track->Charge();
@@ -523,27 +627,26 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeEMCALProbability  (const Al
     matchedClus = (AliVCluster*)fCurrentEvent->GetCaloCluster(nMatchClus);
     
     if(matchedClus){    
-
-    // matched cluster is EMCAL
-    if(matchedClus->IsEMCAL()){
-
+      
+      // matched cluster is EMCAL
+      if(matchedClus->IsEMCAL()){
+      
       fClsE       = matchedClus->E();
       EovP        = fClsE/mom;
       
       
       // compute the probabilities 
-      if( 999 != fEMCALResponse.ComputeEMCALProbability(pt,EovP,charge,p)){    
-
-       // in case everything is OK
-       return kDetPidOk;
-       
+        if(fEMCALResponse.ComputeEMCALProbability(nSpecies,pt,EovP,charge,p)){ 
+          
+          // in case everything is OK
+             return kDetPidOk;
+        }
       }
     }
   }
-  }
   
   // in all other cases set flat distribution (no decision)
-  for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
+  for (Int_t j=0; j<nSpecies; j++) p[j] = 1./nSpecies;
   return kDetNoSignal;
   
 }
@@ -553,7 +656,12 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputePHOSProbability (const AliV
   //
   // Compute PID response for the PHOS
   //
-
+  
+  // look for cached value first
+//   if (track->GetDetectorPID()){
+//     return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies);
+//   }
+  
   // set flat distribution (no decision)
   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
   return kDetNoSignal;
@@ -565,6 +673,12 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliV
   // Compute PID response for the HMPID
   //
 
+
+  // look for cached value first
+  if (track->GetDetectorPID()){
+    return track->GetDetectorPID()->GetRawProbability(kHMPID, p, nSpecies);
+  }
+  
   // set flat distribution (no decision)
   for (Int_t j=0; j<nSpecies; j++) p[j]=1./nSpecies;
   if((track->GetStatus()&AliVTrack::kHMPIDpid)==0) return kDetNoSignal;
@@ -575,7 +689,7 @@ AliPIDResponse::EDetPidStatus AliPIDResponse::ComputeHMPIDProbability(const AliV
 }
 
 //______________________________________________________________________________
-void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
+void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run)
 {
   //
   // Apply settings for the current event
@@ -585,7 +699,8 @@ void AliPIDResponse::InitialiseEvent(AliVEvent *event, Int_t pass)
   fCurrentEvent=0x0;
   if (!event) return;
   fCurrentEvent=event;
-  fRun=event->GetRunNumber();
+  if (run>0) fRun=run;
+  else fRun=event->GetRunNumber();
   
   if (fRun!=fOldRun){
     ExecNewRun();
@@ -689,6 +804,10 @@ void AliPIDResponse::SetRecoInfo()
   //TODO: periods 11B, 11C are not yet treated assume 11d for the moment
   else if (fRun>=148531&&fRun<=155384) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
   else if (fRun>=156477&&fRun<=159635) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+  // also for 11e,f use 11d
+  else if (fRun>=160676&&fRun<=162740) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+  else if (fRun>=162933&&fRun<=165746) { fLHCperiod="LHC11D"; fMCperiodTPC="LHC10F6A"; }
+  
   else if (fRun>=166529) {
     fLHCperiod="LHC11H";
     fMCperiodTPC="LHC11A10";
@@ -700,7 +819,7 @@ void AliPIDResponse::SetRecoInfo()
   //exception new pp MC productions from 2011
   if (fBeamType=="PP" && reg.MatchB(fCurrentFile)) fMCperiodTPC="LHC11B2";
   // exception for 11f1
-  if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";                                                                                                                
+  if (fCurrentFile.Contains("LHC11f1/")) fMCperiodTPC="LHC11F1";
 }
 
 //______________________________________________________________________________
@@ -726,6 +845,7 @@ void AliPIDResponse::SetTPCPidResponseMaster()
   fArrPidResponseMaster=0x0;
   
   TString fileName(Form("%s/COMMON/PID/data/TPCPIDResponse.root", fOADBPath.Data()));
+  if (!fCustomTPCpidResponse.IsNull()) fileName=fCustomTPCpidResponse;
   
   TFile *f=TFile::Open(fileName.Data());
   if (f && f->IsOpen() && !f->IsZombie()){
@@ -767,7 +887,7 @@ void AliPIDResponse::SetTPCParametrisation()
   //
   //reset old splines
   //
-  for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
+  for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
     fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,0x0);
   }
 
@@ -802,16 +922,18 @@ void AliPIDResponse::SetTPCParametrisation()
       TString particleName=arr->At(1)->GetName();
       delete arr;
       if (particleName.IsNull()) continue;
-      if (particleName=="ALL") grAll=responseFunction;
+      if (!grAll&&particleName=="ALL") grAll=responseFunction;
       else {
         //find particle id
-        for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
+        for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
           TString particle=AliPID::ParticleName(ispec);
           particle.ToUpper();
           if ( particle == particleName ){
             fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,responseFunction);
             fTPCResponse.SetUseDatabase(kTRUE);
             AliInfo(Form("Adding graph: %d - %s",ispec,responseFunction->GetName()));
+            // overwrite default with proton spline (for light nuclei)
+            if (ispec==AliPID::kProton) grAll=responseFunction;
             found=kTRUE;
             break;
           }
@@ -821,7 +943,7 @@ void AliPIDResponse::SetTPCParametrisation()
     
     //set default response function to all particles which don't have a specific one
     if (grAll){
-      for (Int_t ispec=0; ispec<AliPID::kSPECIES; ++ispec){
+      for (Int_t ispec=0; ispec<AliPID::kSPECIESC; ++ispec){
         if (!fTPCResponse.GetResponseFunction((AliPID::EParticleType)ispec)){
           fTPCResponse.SetResponseFunction((AliPID::EParticleType)ispec,grAll);
           AliInfo(Form("Adding graph: %d - %s",ispec,grAll->GetName()));
@@ -920,26 +1042,35 @@ void AliPIDResponse::SetTOFPidResponseMaster()
   //
   // Load the TOF pid params from the OADB
   //
+
+  if (fTOFPIDParams) delete fTOFPIDParams;
+  fTOFPIDParams=0x0;
+
   TFile *oadbf = new TFile(Form("%s/COMMON/PID/data/TOFPIDParams.root",fOADBPath.Data()));
-  if (oadbf->IsOpen()) {
+  if (oadbf && oadbf->IsOpen()) {
     AliInfo(Form("Loading TOF Params from %s/COMMON/PID/data/TOFPIDParams.root", fOADBPath.Data()));
     AliOADBContainer *oadbc = (AliOADBContainer *)oadbf->Get("TOFoadb");
-    if (fTOFPIDParams) delete fTOFPIDParams;
-    fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
+    if (oadbc) fTOFPIDParams = dynamic_cast<AliTOFPIDParams *>(oadbc->GetObject(fRun,"TOFparams"));
     oadbf->Close();
     delete oadbc;
-  } else {
-    AliError(Form("TOFPIDParams.root not found in %s/COMMON/PID/data !!",fOADBPath.Data()));
   }
   delete oadbf;
 
-  }
+  if (!fTOFPIDParams) AliFatal("TOFPIDParams could not be retrieved");
+}
 
 //______________________________________________________________________________
 void AliPIDResponse::InitializeTOFResponse(){
   //
   // Set PID Params to the TOF PID response
-  // 
+  //
+
+  AliInfo("TOF PID Params loaded from OADB");
+  AliInfo(Form("  TOF resolution %5.2f [ps]",fTOFPIDParams->GetTOFresolution()));
+  AliInfo(Form("  StartTime method %d",fTOFPIDParams->GetStartTimeMethod()));
+  AliInfo(Form("  TOF res. mom. params: %5.2f %5.2f %5.2f %5.2f",
+               fTOFPIDParams->GetSigParams(0),fTOFPIDParams->GetSigParams(1),fTOFPIDParams->GetSigParams(2),fTOFPIDParams->GetSigParams(3)));
+  
   for (Int_t i=0;i<4;i++) {
     fTOFResponse.SetTrackParameter(i,fTOFPIDParams->GetSigParams(i));
   }
@@ -1016,6 +1147,40 @@ void AliPIDResponse::InitializeEMCALResponse(){
   fEMCALResponse.SetPIDParams(fEMCALPIDParams);
 
 }
+
+//_____________________________________________________
+void AliPIDResponse::FillTrackDetectorPID()
+{
+  //
+  // create detector PID information and setup the transient pointer in the track
+  //
+
+  if (!fCurrentEvent) return;
+  
+  //TODO: which particles to include? See also the loops below...
+  Double_t values[AliPID::kSPECIESC]={0};
+  
+  for (Int_t itrack=0; itrack<fCurrentEvent->GetNumberOfTracks(); ++itrack){
+    AliVTrack *track=dynamic_cast<AliVTrack*>(fCurrentEvent->GetTrack(itrack));
+    if (!track) continue;
+
+    AliDetectorPID *detPID=new AliDetectorPID;
+    for (Int_t idet=0; idet<kNdetectors; ++idet){
+
+      //nsigmas
+      for (Int_t ipart=0; ipart<AliPID::kSPECIESC; ++ipart)
+        values[ipart]=NumberOfSigmas((EDetector)idet,track,(AliPID::EParticleType)ipart);
+      detPID->SetNumberOfSigmas((EDetector)idet, values, (Int_t)AliPID::kSPECIESC);
+
+      //probabilities
+      EDetPidStatus status=ComputePIDProbability((EDetector)idet,track,AliPID::kSPECIESC,values);
+      detPID->SetRawProbability((EDetector)idet, values, (Int_t)AliPID::kSPECIESC, status);
+    }
+
+    track->SetDetectorPID(detPID);
+  }
+}
+
 //_________________________________________________________________________
 void AliPIDResponse::SetTOFResponse(AliVEvent *vevent,EStartTimeType_t option){
   //
index df58190..307d5b4 100644 (file)
 
 #include "TNamed.h"
 
-class AliVEvent;
 class TF1;
-class AliTRDPIDResponseObject; 
+class TObjArray;
+
+class AliVEvent;
+class AliTRDPIDResponseObject;
+class AliTOFPIDParams;
 
 class AliPIDResponse : public TNamed {
 public:
   AliPIDResponse(Bool_t isMC=kFALSE);
   virtual ~AliPIDResponse();
 
+  enum EDetector {
+    kITS=0,
+    kTPC=1,
+    kTRD=2,
+    kTOF=3,
+    kHMPID=4,
+    kEMCAL=5,
+    kPHOS=6,
+    kNdetectors=7
+  };
+  
   enum EDetCode {
     kDetITS = 0x1,
     kDetTPC = 0x2,
@@ -58,18 +72,21 @@ public:
   AliTRDPIDResponse &GetTRDResponse() {return fTRDResponse;}
   AliEMCALPIDResponse &GetEMCALResponse() {return fEMCALResponse;}
 
-  Float_t NumberOfSigmas(EDetCode detCode, const AliVParticle *track, AliPID::EParticleType type) const;
+  Float_t NumberOfSigmas(EDetector detCode, const AliVParticle *track, AliPID::EParticleType type) const;
+  Float_t NumberOfSigmas(EDetCode  detCode, const AliVParticle *track, AliPID::EParticleType type) const;
   
-  virtual Float_t NumberOfSigmasITS(const AliVParticle *track, AliPID::EParticleType type) const;
-  virtual Float_t NumberOfSigmasTPC(const AliVParticle *track, AliPID::EParticleType type) const;
-  virtual Float_t NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type) const;
-  virtual Float_t NumberOfSigmasEMCAL(const AliVTrack *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const;
-  virtual Float_t NumberOfSigmasTOF(const AliVParticle *track, AliPID::EParticleType type) const = 0;
+  virtual Float_t NumberOfSigmasITS  (const AliVParticle *track, AliPID::EParticleType type) const;
+  virtual Float_t NumberOfSigmasTPC  (const AliVParticle *track, AliPID::EParticleType type) const;
+  virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type, Double_t &eop, Double_t showershape[4]) const;
+  virtual Float_t NumberOfSigmasTOF  (const AliVParticle *track, AliPID::EParticleType type) const = 0;
+  virtual Float_t NumberOfSigmasEMCAL(const AliVParticle *track, AliPID::EParticleType type) const;
+
   virtual Bool_t IdentifiedAsElectronTRD(const AliVTrack *track, Double_t efficiencyLevel) const;
 
-  EDetPidStatus ComputePIDProbability  (EDetCode detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputePIDProbability  (EDetector detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
+  EDetPidStatus ComputePIDProbability  (EDetCode  detCode, const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
   
-  EDetPidStatus ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t  p[]) const;
+  EDetPidStatus ComputeITSProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
   EDetPidStatus ComputeTPCProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
   EDetPidStatus ComputeTOFProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
   EDetPidStatus ComputeTRDProbability  (const AliVTrack *track, Int_t nSpecies, Double_t p[]) const;
@@ -84,9 +101,16 @@ public:
   
   void SetOADBPath(const char* path) {fOADBPath=path;}
   const char *GetOADBPath() const {return fOADBPath.Data();}
-  void InitialiseEvent(AliVEvent *event, Int_t pass);
+
+  void SetCustomTPCpidResponse(const char* tpcpid) { fCustomTPCpidResponse = tpcpid; }
+  const char* GetCustomTPCpidResponse() const { return fCustomTPCpidResponse.Data(); }
+  
+  void InitialiseEvent(AliVEvent *event, Int_t pass, Int_t run=-1);
   void SetCurrentFile(const char* file) { fCurrentFile=file; }
 
+  // cache PID in the track
+  void FillTrackDetectorPID();
+
   AliVEvent * GetCurrentEvent() const {return fCurrentEvent;}
 
   // User settings for the MC period and reco pass
@@ -122,6 +146,7 @@ private:
   Bool_t fIsMC;                        //  If we run on MC data
 
   TString fOADBPath;                   // OADB path to use
+  TString fCustomTPCpidResponse;       // Custom TPC Pid Response file for debugging purposes
   
   TString fBeamType;                   //! beam type (PP) or (PBPB)
   TString fLHCperiod;                  //! LHC period
@@ -179,37 +204,7 @@ private:
   //
   void SetRecoInfo();
   
-  ClassDef(AliPIDResponse, 8);  //PID response handling
+  ClassDef(AliPIDResponse, 9);  //PID response handling
 };
 
-inline Float_t AliPIDResponse::NumberOfSigmasTPC(const AliVParticle *vtrack, AliPID::EParticleType type) const {
-  AliVTrack *track=(AliVTrack*)vtrack;
-  Double_t mom  = track->GetTPCmomentum();
-  Double_t sig  = track->GetTPCsignal();
-  UInt_t   sigN = track->GetTPCsignalN();
-
-  if(fTuneMConData) sig = GetTPCsignalTunedOnData(track);
-
-  Double_t nSigma = -999.;
-  if (sigN>0) nSigma=fTPCResponse.GetNumberOfSigmas(mom,sig,sigN,type);
-
-  return nSigma;
-}
-
-inline Float_t AliPIDResponse::NumberOfSigmasITS(const AliVParticle *vtrack, AliPID::EParticleType type) const {
-  AliVTrack *track=(AliVTrack*)vtrack;
-  Float_t dEdx=track->GetITSsignal();
-  if (dEdx<=0) return -999.;
-  
-  UChar_t clumap=track->GetITSClusterMap();
-  Int_t nPointsForPid=0;
-  for(Int_t i=2; i<6; i++){
-    if(clumap&(1<<i)) ++nPointsForPid;
-  }
-  Float_t mom=track->P();
-  Bool_t isSA=kTRUE;
-  if(track->GetTPCNcls()>0) isSA=kFALSE;
-  return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA);
-}
-
 #endif
index 93811c2..b0bf074 100644 (file)
@@ -25,6 +25,7 @@
 #include <TGraph.h>
 #include <TObjArray.h>
 #include <TSpline.h>
+#include <TMath.h>
 
 #include "AliExternalTrackParam.h"
 
@@ -129,8 +130,12 @@ Double_t AliTPCPIDResponse::GetExpectedSignal(const Float_t mom,
   if (!fUseDatabase||fResponseFunctions.GetEntriesFast()>AliPID::kUnknown) return Bethe(mom/mass);
   //
   TSpline3 * responseFunction = (TSpline3 *) fResponseFunctions.UncheckedAt(n);
+
   if (!responseFunction) return Bethe(mom/mass);
-  return fMIP*responseFunction->Eval(mom/mass);
+  //charge factor. BB goes with z^2, however in reality it is slightly larger (calibration, threshold effects, ...)
+  // !!! Splines for light nuclei need to be normalised to this factor !!!
+  const Double_t chargeFactor = TMath::Power(AliPID::ParticleCharge(n),2.3);
+  return fMIP*responseFunction->Eval(mom/mass)*chargeFactor;
 
 }
 
index ae1af45..bb40010 100644 (file)
@@ -16,6 +16,7 @@
 class AliVVertex;\r
 class AliExternalTrackParam;\r
 class AliTPCdEdxInfo;\r
+class AliDetectorPID;\r
  \r
 class AliVTrack: public AliVParticle {\r
 \r
@@ -60,6 +61,8 @@ public:
   virtual Double_t GetTRDslice(Int_t /*plane*/, Int_t /*slice*/) const { return -1.; }\r
   virtual Int_t    GetNumberOfTRDslices() const { return 0; }\r
   virtual UChar_t  GetTRDntrackletsPID() const { return 0;}\r
+  virtual void     SetDetectorPID(const AliDetectorPID */*pid*/) {;}\r
+  virtual const    AliDetectorPID* GetDetectorPID() const { return 0x0; }\r
   virtual Double_t GetTRDchi2()          const { return -1;}\r
   \r
   virtual Int_t GetEMCALcluster()     const {return -1;}\r
index e9906bd..e4240fa 100644 (file)
@@ -99,6 +99,8 @@
 #pragma link C++ class AliTRDPIDResponse+;
 #pragma link C++ class AliEMCALPIDResponse+;
 #pragma link C++ class AliPIDCombined+;
+#pragma link C++ class AliPIDValues+;
+#pragma link C++ class AliDetectorPID+;
 #pragma link C++ class AliTOFHeader+;
 
 #pragma link C++ class AliDAQ+;