From 00a38d0729217a72e2a4c3b7dae1f90aa2cbef58 Mon Sep 17 00:00:00 2001 From: morsch Date: Thu, 23 Aug 2012 16:08:32 +0000 Subject: [PATCH] - Reshuffling of the particle codes in AliPID. Now the light nuclei are between the 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 --- ANALYSIS/AliAnalysisTaskESDfilter.cxx | 2 +- ANALYSIS/AliAnalysisTaskPIDResponse.cxx | 10 +- ANALYSIS/AliAnalysisTaskPIDResponse.h | 7 +- ANALYSIS/AliAnalysisTaskPIDqa.cxx | 30 +- ANALYSIS/TenderSupplies/AddTaskTender.C | 38 ++- .../TenderSupplies/AliPIDTenderSupply.cxx | 13 +- ANALYSIS/TenderSupplies/AliPIDTenderSupply.h | 6 +- CORRFW/AliCFTrackCutPid.cxx | 8 +- CORRFW/AliCFTrackCutPid.h | 2 +- EMCAL/AliEMCALPID.cxx | 4 +- EMCAL/AliEMCALPIDUtils.cxx | 36 +- EMCAL/AliEMCALPIDUtils.h | 6 +- EMCAL/AliEMCALRecoUtils.cxx | 4 +- HLT/BASE/util/AliHLTCaloClusterDataStruct.h | 4 +- HLT/CALO/AliHLTCaloClusterAnalyser.cxx | 2 +- HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx | 2 +- PHOS/AliPHOSFastRecParticle.cxx | 6 +- PHOS/AliPHOSFastRecParticle.h | 2 +- PHOS/AliPHOSPIDv1.cxx | 8 +- PHOS/AliPHOSPIDv1.h | 2 +- PHOS/AliPHOSRecParticle.cxx | 4 +- PWG/CaloTrackCorrBase/AliCaloPID.cxx | 6 +- PWGHF/hfe/AliHFEdca.cxx | 2 +- STEER/AOD/AliAODPid.cxx | 12 +- STEER/AOD/AliAODPid.h | 2 +- STEER/AOD/AliAODTrack.cxx | 23 +- STEER/AOD/AliAODTrack.h | 8 +- STEER/AOD/AliAODpidUtil.cxx | 48 ++- STEER/AOD/AliAODpidUtil.h | 30 -- STEER/CMakelibSTEERBase.pkg | 2 + STEER/ESD/AliESDCaloCluster.cxx | 8 +- STEER/ESD/AliESDCaloCluster.h | 6 +- STEER/ESD/AliESDpid.cxx | 19 +- STEER/ESD/AliESDpid.h | 13 +- STEER/ESD/AliESDtrack.cxx | 23 +- STEER/ESD/AliESDtrack.h | 11 +- STEER/STEERBase/AliEMCALPIDResponse.cxx | 23 +- STEER/STEERBase/AliEMCALPIDResponse.h | 4 +- STEER/STEERBase/AliITSPIDResponse.h | 3 +- STEER/STEERBase/AliPID.cxx | 93 +++--- STEER/STEERBase/AliPID.h | 55 ++-- STEER/STEERBase/AliPIDCombined.cxx | 35 +- STEER/STEERBase/AliPIDCombined.h | 6 +- STEER/STEERBase/AliPIDResponse.cxx | 309 ++++++++++++++---- STEER/STEERBase/AliPIDResponse.h | 79 +++-- STEER/STEERBase/AliTPCPIDResponse.cxx | 7 +- STEER/STEERBase/AliVTrack.h | 3 + STEER/STEERBaseLinkDef.h | 2 + 48 files changed, 664 insertions(+), 364 deletions(-) diff --git a/ANALYSIS/AliAnalysisTaskESDfilter.cxx b/ANALYSIS/AliAnalysisTaskESDfilter.cxx index 4dac86a92ae..90e4ae6c110 100644 --- a/ANALYSIS/AliAnalysisTaskESDfilter.cxx +++ b/ANALYSIS/AliAnalysisTaskESDfilter.cxx @@ -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()); diff --git a/ANALYSIS/AliAnalysisTaskPIDResponse.cxx b/ANALYSIS/AliAnalysisTaskPIDResponse.cxx index 4ff0d91cd97..f0a2f086322 100644 --- a/ANALYSIS/AliAnalysisTaskPIDResponse.cxx +++ b/ANALYSIS/AliAnalysisTaskPIDResponse.cxx @@ -27,7 +27,6 @@ #include #include - #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(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(); + } } //______________________________________________________________________________ diff --git a/ANALYSIS/AliAnalysisTaskPIDResponse.h b/ANALYSIS/AliAnalysisTaskPIDResponse.h index 2894ae7a551..71cdd660397 100644 --- a/ANALYSIS/AliAnalysisTaskPIDResponse.h +++ b/ANALYSIS/AliAnalysisTaskPIDResponse.h @@ -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 diff --git a/ANALYSIS/AliAnalysisTaskPIDqa.cxx b/ANALYSIS/AliAnalysisTaskPIDqa.cxx index e1eff32802d..54b8ff1d8b6 100644 --- a/ANALYSIS/AliAnalysisTaskPIDqa.cxx +++ b/ANALYSIS/AliAnalysisTaskPIDqa.cxx @@ -381,13 +381,13 @@ void AliAnalysisTaskPIDqa::FillITSqa() } - for (Int_t ispecie=0; ispecieAt(ispecie); if (!h) continue; Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie); h->Fill(mom,nSigma); } - TH2 *h=(TH2*)theList->At(AliPID::kSPECIES); + TH2 *h=(TH2*)theList->At(AliPID::kSPECIESC); if (h) { Double_t sig=track->GetITSsignal(); h->Fill(mom,sig); @@ -429,14 +429,14 @@ void AliAnalysisTaskPIDqa::FillTPCqa() Double_t mom=track->GetTPCmomentum(); - for (Int_t ispecie=0; ispecieAt(ispecie); if (!h) continue; Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie); h->Fill(mom,nSigma); } - TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES); + TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIESC); if (h) { Double_t sig=track->GetTPCsignal(); h->Fill(mom,sig); @@ -529,7 +529,7 @@ void AliAnalysisTaskPIDqa::FillTOFqa() Double_t mom=track->P(); - for (Int_t ispecie=0; ispecieAt(ispecie); if (!h) continue; Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie); @@ -875,7 +875,7 @@ void AliAnalysisTaskPIDqa::FillTPCTOFqa() Double_t mom=track->P(); Double_t momTPC=track->GetTPCmomentum(); - for (Int_t ispecie=0; ispecieNumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie); Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie); @@ -885,11 +885,11 @@ void AliAnalysisTaskPIDqa::FillTPCTOFqa() if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC); //TOF after TPC cut - h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES); + h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC); if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF); //EMCAL after TOF and TPC cut - h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES); + h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC); if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){ Int_t nMatchClus = track->GetEMCALcluster(); @@ -929,7 +929,7 @@ void AliAnalysisTaskPIDqa::SetupITSqa() TVectorD *vX=MakeLogBinning(200,.1,30); //ITS+TPC tracks - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), @@ -943,7 +943,7 @@ void AliAnalysisTaskPIDqa::SetupITSqa() fListQAits->Add(hSig); //ITS Standalone tracks - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), @@ -957,7 +957,7 @@ void AliAnalysisTaskPIDqa::SetupITSqa() fListQAitsSA->Add(hSigSA); //ITS Pure Standalone tracks - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), @@ -982,7 +982,7 @@ void AliAnalysisTaskPIDqa::SetupTPCqa() TVectorD *vX=MakeLogBinning(200,.1,30); - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), @@ -1028,7 +1028,7 @@ void AliAnalysisTaskPIDqa::SetupTOFqa() TVectorD *vX=MakeLogBinning(200,.1,30); - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), @@ -1158,7 +1158,7 @@ void AliAnalysisTaskPIDqa::SetupTPCTOFqa() TVectorD *vX=MakeLogBinning(200,.1,30); //TPC signals after TOF cut - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), @@ -1167,7 +1167,7 @@ void AliAnalysisTaskPIDqa::SetupTPCTOFqa() } //TOF signals after TPC cut - for (Int_t ispecie=0; ispecieGetNrows()-1,vX->GetMatrixArray(), diff --git a/ANALYSIS/TenderSupplies/AddTaskTender.C b/ANALYSIS/TenderSupplies/AddTaskTender.C index e1d15751fa0..2831f2bca40 100644 --- a/ANALYSIS/TenderSupplies/AddTaskTender.C +++ b/ANALYSIS/TenderSupplies/AddTaskTender.C @@ -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 //================================================ diff --git a/ANALYSIS/TenderSupplies/AliPIDTenderSupply.cxx b/ANALYSIS/TenderSupplies/AliPIDTenderSupply.cxx index 6b567c90ab8..cadf2133ffe 100644 --- a/ANALYSIS/TenderSupplies/AliPIDTenderSupply.cxx +++ b/ANALYSIS/TenderSupplies/AliPIDTenderSupply.cxx @@ -29,8 +29,11 @@ #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 // diff --git a/ANALYSIS/TenderSupplies/AliPIDTenderSupply.h b/ANALYSIS/TenderSupplies/AliPIDTenderSupply.h index a82dc492751..59d0330b887 100644 --- a/ANALYSIS/TenderSupplies/AliPIDTenderSupply.h +++ b/ANALYSIS/TenderSupplies/AliPIDTenderSupply.h @@ -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 }; diff --git a/CORRFW/AliCFTrackCutPid.cxx b/CORRFW/AliCFTrackCutPid.cxx index a57cdf8528f..7d94b0efc4c 100644 --- a/CORRFW/AliCFTrackCutPid.cxx +++ b/CORRFW/AliCFTrackCutPid.cxx @@ -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++){ diff --git a/CORRFW/AliCFTrackCutPid.h b/CORRFW/AliCFTrackCutPid.h index 3e82fcc9f22..2b4df8dc07a 100644 --- a/CORRFW/AliCFTrackCutPid.h +++ b/CORRFW/AliCFTrackCutPid.h @@ -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 diff --git a/EMCAL/AliEMCALPID.cxx b/EMCAL/AliEMCALPID.cxx index 71cd6631c4a..ec26872d7b9 100644 --- a/EMCAL/AliEMCALPID.cxx +++ b/EMCAL/AliEMCALPID.cxx @@ -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; iGetPIDFinal(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=0&&idx<3) return fPID[idx]; else return 0.;} - Double_t GetPIDFinal(Int_t idx) const {if (idx>=0&&idx=0&&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=0&&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 diff --git a/EMCAL/AliEMCALRecoUtils.cxx b/EMCAL/AliEMCALRecoUtils.cxx index dac3c432e0a..56cbc677c9e 100644 --- a/EMCAL/AliEMCALRecoUtils.cxx +++ b/EMCAL/AliEMCALRecoUtils.cxx @@ -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); } diff --git a/HLT/BASE/util/AliHLTCaloClusterDataStruct.h b/HLT/BASE/util/AliHLTCaloClusterDataStruct.h index ab25cc9fe99..fefa58ddc5f 100644 --- a/HLT/BASE/util/AliHLTCaloClusterDataStruct.h +++ b/HLT/BASE/util/AliHLTCaloClusterDataStruct.h @@ -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 diff --git a/HLT/CALO/AliHLTCaloClusterAnalyser.cxx b/HLT/CALO/AliHLTCaloClusterAnalyser.cxx index 04267b63cb7..66fee315f83 100644 --- a/HLT/CALO/AliHLTCaloClusterAnalyser.cxx +++ b/HLT/CALO/AliHLTCaloClusterAnalyser.cxx @@ -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; } diff --git a/HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx b/HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx index a0d18a57250..d1302d69665 100644 --- a/HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx +++ b/HLT/PHOS/AliHLTPHOSClusterAnalyser.cxx @@ -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; } diff --git a/PHOS/AliPHOSFastRecParticle.cxx b/PHOS/AliPHOSFastRecParticle.cxx index 2bea0db4762..6ef6bfbf9b0 100644 --- a/PHOS/AliPHOSFastRecParticle.cxx +++ b/PHOS/AliPHOSFastRecParticle.cxx @@ -52,7 +52,7 @@ AliPHOSFastRecParticle::AliPHOSFastRecParticle() : { // ctor - for(Int_t i=0; iSetParameters(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; iIsEMCAL() && 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); diff --git a/PWGHF/hfe/AliHFEdca.cxx b/PWGHF/hfe/AliHFEdca.cxx index 5223146d7cf..1f27dec33fc 100644 --- a/PWGHF/hfe/AliHFEdca.cxx +++ b/PWGHF/hfe/AliHFEdca.cxx @@ -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; diff --git a/STEER/AOD/AliAODPid.cxx b/STEER/AOD/AliAODPid.cxx index 81a385781b9..321fcc412e5 100644 --- a/STEER/AOD/AliAODPid.cxx +++ b/STEER/AOD/AliAODPid.cxx @@ -43,7 +43,7 @@ AliAODPid::AliAODPid(): fTPCdEdxInfo(0) { // default constructor - for(Int_t i=0; i(*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 { diff --git a/STEER/AOD/AliAODTrack.h b/STEER/AOD/AliAODTrack.h index 7d4e44c17f2..1a5487e4644 100644 --- a/STEER/AOD/AliAODTrack.h +++ b/STEER/AOD/AliAODTrack.h @@ -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 diff --git a/STEER/AOD/AliAODpidUtil.cxx b/STEER/AOD/AliAODpidUtil.cxx index 34c022fbeb3..d3c1d4035a1 100644 --- a/STEER/AOD/AliAODpidUtil.cxx +++ b/STEER/AOD/AliAODpidUtil.cxx @@ -34,6 +34,8 @@ #include "AliAODMCHeader.h" #include "AliAODMCParticle.h" +#include + 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.; +} diff --git a/STEER/AOD/AliAODpidUtil.h b/STEER/AOD/AliAODpidUtil.h index 7f3718a5080..c962d9bfb1a 100644 --- a/STEER/AOD/AliAODpidUtil.h +++ b/STEER/AOD/AliAODpidUtil.h @@ -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 diff --git a/STEER/CMakelibSTEERBase.pkg b/STEER/CMakelibSTEERBase.pkg index ef4d53ded76..c3187b35ac8 100644 --- a/STEER/CMakelibSTEERBase.pkg +++ b/STEER/CMakelibSTEERBase.pkg @@ -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 diff --git a/STEER/ESD/AliESDCaloCluster.cxx b/STEER/ESD/AliESDCaloCluster.cxx index 4aaa0bc03c2..e05f7d46e84 100644 --- a/STEER/ESD/AliESDCaloCluster.cxx +++ b/STEER/ESD/AliESDCaloCluster.cxx @@ -53,7 +53,7 @@ AliESDCaloCluster::AliESDCaloCluster() : // The default ESD constructor // fGlobalPos[0] = fGlobalPos[1] = fGlobalPos[2] = 0.; - for(Int_t i=0; i 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=0 && i=0 && i 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)); +} diff --git a/STEER/ESD/AliESDpid.h b/STEER/ESD/AliESDpid.h index 4d272fd528e..4d9a3a423ea 100644 --- a/STEER/ESD/AliESDpid.h +++ b/STEER/ESD/AliESDpid.h @@ -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 diff --git a/STEER/ESD/AliESDtrack.cxx b/STEER/ESD/AliESDtrack.cxx index 0a148c2b643..68a77bda57d 100644 --- a/STEER/ESD/AliESDtrack.cxx +++ b/STEER/ESD/AliESDtrack.cxx @@ -128,6 +128,7 @@ #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 Probability from Gaussian distribution // // (proper normalization to each other?) // // // -// NO Parametrization (outside pT range): --> return 999 // +// NO Parametrization (outside pT range): --> return kFALSE // ////////////////////////////////////////////////////////////////////////// #include @@ -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; } //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/STEER/STEERBase/AliEMCALPIDResponse.h b/STEER/STEERBase/AliEMCALPIDResponse.h index 52cf1b8c1a6..57c17fb8abb 100644 --- a/STEER/STEERBase/AliEMCALPIDResponse.h +++ b/STEER/STEERBase/AliEMCALPIDResponse.h @@ -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: diff --git a/STEER/STEERBase/AliITSPIDResponse.h b/STEER/STEERBase/AliITSPIDResponse.h index c4b63d60d52..3556220b597 100644 --- a/STEER/STEERBase/AliITSPIDResponse.h +++ b/STEER/STEERBase/AliITSPIDResponse.h @@ -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; diff --git a/STEER/STEERBase/AliPID.cxx b/STEER/STEERBase/AliPID.cxx index 7c5b2613a9e..96d203023de 100644 --- a/STEER/STEERBase/AliPID.cxx +++ b/STEER/STEERBase/AliPID.cxx @@ -50,79 +50,90 @@ 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; diff --git a/STEER/STEERBase/AliPID.h b/STEER/STEERBase/AliPID.h index e2544121573..c4b6875f02b 100644 --- a/STEER/STEERBase/AliPID.h +++ b/STEER/STEERBase/AliPID.h @@ -14,14 +14,14 @@ #include - +#include 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 }; diff --git a/STEER/STEERBase/AliPIDCombined.cxx b/STEER/STEERBase/AliPIDCombined.cxx index 15a6556d529..248a7f894b7 100644 --- a/STEER/STEERBase/AliPIDCombined.cxx +++ b/STEER/STEERBase/AliPIDCombined.cxx @@ -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::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; } diff --git a/STEER/STEERBase/AliPIDCombined.h b/STEER/STEERBase/AliPIDCombined.h index ce1a52b2ec7..d7704b26d4a 100644 --- a/STEER/STEERBase/AliPIDCombined.h +++ b/STEER/STEERBase/AliPIDCombined.h @@ -60,11 +60,11 @@ private: Int_t fRejectMismatchMask; // Detectors set return flat prob. if mismatch detected Bool_t fEnablePriors; // Enable bayesian PID (if kFALSE priors set flat) Int_t fSelectedSpecies; // Number of selected species to study - TH1F *fPriorsDistributions[AliPID::kSPECIES+AliPID::kSPECIESLN]; // priors + TH1F *fPriorsDistributions[AliPID::kSPECIESC]; // priors Bool_t fUseDefaultTPCPriors; // switch to use Defaul TPC Priors - static TH2F *fDefaultPriorsTPC[5]; // Default priors for TPC tracks + static TH2F *fDefaultPriorsTPC[AliPID::kSPECIES]; // Default priors for TPC tracks - ClassDef(AliPIDCombined,2); + ClassDef(AliPIDCombined, 3); // Combined PID using priors }; #endif diff --git a/STEER/STEERBase/AliPIDResponse.cxx b/STEER/STEERBase/AliPIDResponse.cxx index 491b6b0eb6a..08f33d27cdb 100644 --- a/STEER/STEERBase/AliPIDResponse.cxx +++ b/STEER/STEERBase/AliPIDResponse.cxx @@ -29,6 +29,7 @@ #include #include #include +#include #include #include @@ -40,6 +41,7 @@ #include #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<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<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; @@ -306,6 +389,16 @@ 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<GetDetectorPID()){ + return track->GetDetectorPID()->GetRawProbability(kITS, p, nSpecies); + } + // set flat distribution (no decision) for (Int_t j=0; jGetStatus()&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 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; jGetTPCsignalTunedOnData(track); - for (Int_t j=0; jGetTPCsignalN(),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; jP(),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; jGetStatus()&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; jGetEMCALcluster(); 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; jGetDetectorPID()){ +// return track->GetDetectorPID()->GetRawProbability(kPHOS, p, nSpecies); +// } + // set flat distribution (no decision) for (Int_t j=0; jGetDetectorPID()){ + return track->GetDetectorPID()->GetRawProbability(kHMPID, p, nSpecies); + } + // set flat distribution (no decision) for (Int_t j=0; jGetStatus()&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; ispecAt(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; ispecGetName())); + // 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; ispecGetName())); @@ -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(oadbc->GetObject(fRun,"TOFparams")); + if (oadbc) fTOFPIDParams = dynamic_cast(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; itrackGetNumberOfTracks(); ++itrack){ + AliVTrack *track=dynamic_cast(fCurrentEvent->GetTrack(itrack)); + if (!track) continue; + + AliDetectorPID *detPID=new AliDetectorPID; + for (Int_t idet=0; idetSetNumberOfSigmas((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){ // diff --git a/STEER/STEERBase/AliPIDResponse.h b/STEER/STEERBase/AliPIDResponse.h index df581906260..307d5b497eb 100644 --- a/STEER/STEERBase/AliPIDResponse.h +++ b/STEER/STEERBase/AliPIDResponse.h @@ -23,15 +23,29 @@ #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<P(); - Bool_t isSA=kTRUE; - if(track->GetTPCNcls()>0) isSA=kFALSE; - return fITSResponse.GetNumberOfSigmas(mom,dEdx,type,nPointsForPid,isSA); -} - #endif diff --git a/STEER/STEERBase/AliTPCPIDResponse.cxx b/STEER/STEERBase/AliTPCPIDResponse.cxx index 93811c2bc0e..b0bf074436f 100644 --- a/STEER/STEERBase/AliTPCPIDResponse.cxx +++ b/STEER/STEERBase/AliTPCPIDResponse.cxx @@ -25,6 +25,7 @@ #include #include #include +#include #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; } diff --git a/STEER/STEERBase/AliVTrack.h b/STEER/STEERBase/AliVTrack.h index ae1af45f366..bb40010cd6f 100644 --- a/STEER/STEERBase/AliVTrack.h +++ b/STEER/STEERBase/AliVTrack.h @@ -16,6 +16,7 @@ class AliVVertex; class AliExternalTrackParam; class AliTPCdEdxInfo; +class AliDetectorPID; class AliVTrack: public AliVParticle { @@ -60,6 +61,8 @@ public: virtual Double_t GetTRDslice(Int_t /*plane*/, Int_t /*slice*/) const { return -1.; } virtual Int_t GetNumberOfTRDslices() const { return 0; } virtual UChar_t GetTRDntrackletsPID() const { return 0;} + virtual void SetDetectorPID(const AliDetectorPID */*pid*/) {;} + virtual const AliDetectorPID* GetDetectorPID() const { return 0x0; } virtual Double_t GetTRDchi2() const { return -1;} virtual Int_t GetEMCALcluster() const {return -1;} diff --git a/STEER/STEERBaseLinkDef.h b/STEER/STEERBaseLinkDef.h index e9906bd8d9f..e4240fa7142 100644 --- a/STEER/STEERBaseLinkDef.h +++ b/STEER/STEERBaseLinkDef.h @@ -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+; -- 2.39.3