From 0792aa82787a0a04e9a9bea3564868dae5b71eba Mon Sep 17 00:00:00 2001 From: dainese Date: Fri, 23 Oct 2009 13:40:06 +0000 Subject: [PATCH] Updates in PID usage (Markus) --- PWG3/hfe/AddTaskHFE.C | 6 +- PWG3/hfe/AliAnalysisTaskHFE.cxx | 96 ++++++++--- PWG3/hfe/AliAnalysisTaskHFE.h | 3 + PWG3/hfe/AliHFEcuts.cxx | 67 +++++++- PWG3/hfe/AliHFEcuts.h | 10 +- PWG3/hfe/AliHFEmcQA.cxx | 275 +++++++++++++++++++++++++++++++- PWG3/hfe/AliHFEmcQA.h | 19 ++- PWG3/hfe/AliHFEpid.cxx | 155 ++++++++++++++++-- PWG3/hfe/AliHFEpid.h | 15 +- PWG3/hfe/AliHFEpidTPC.cxx | 52 ++++-- PWG3/hfe/AliHFEpidTPC.h | 57 ++++--- PWG3/hfe/AliHFEpidTRD.cxx | 37 ++++- PWG3/hfe/AliHFEpidTRD.h | 10 +- PWG3/hfe/AliHFEsecVtx.cxx | 1 + PWG3/hfe/runElectronTask.C | 4 +- 15 files changed, 721 insertions(+), 86 deletions(-) diff --git a/PWG3/hfe/AddTaskHFE.C b/PWG3/hfe/AddTaskHFE.C index 890f34695cb..174987e67b7 100644 --- a/PWG3/hfe/AddTaskHFE.C +++ b/PWG3/hfe/AddTaskHFE.C @@ -22,11 +22,11 @@ AliAnalysisTask *AddTaskHFE(){ AliHFEcuts *hfecuts = new AliHFEcuts; hfecuts->CreateStandardCuts(); hfecuts->SetCutITSpixel(AliHFEextraCuts::kFirst); + hfecuts->SetMinNTrackletsTRD(1); //hfecuts->SetCheckITSLayerStatus(kFALSE); - AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE; - task->AddPIDdetector("TPC"); - //task->SetPIDStrategy(3); + AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE("Heavy Flavour Electron Analysis"); + task->SetPIDStrategy(4); task->SetHFECuts(hfecuts); task->SetQAOn(AliAnalysisTaskHFE::kMCqa); task->SetSecVtxOn(); diff --git a/PWG3/hfe/AliAnalysisTaskHFE.cxx b/PWG3/hfe/AliAnalysisTaskHFE.cxx index 398aef540e2..ce077540a16 100644 --- a/PWG3/hfe/AliAnalysisTaskHFE.cxx +++ b/PWG3/hfe/AliAnalysisTaskHFE.cxx @@ -69,6 +69,43 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(): AliAnalysisTask("PID efficiency Analysis", "") , fQAlevel(0) , fPIDdetectors("") + , fPIDstrategy(0) + , fESD(0x0) + , fMC(0x0) + , fCFM(0x0) + , fCorrelation(0x0) + , fPIDperformance(0x0) + , fPID(0x0) + , fCuts(0x0) + , fSecVtx(0x0) + , fMCQA(0x0) + , fNEvents(0x0) + , fNElectronTracksEvent(0x0) + , fQA(0x0) + , fOutput(0x0) + , fHistMCQA(0x0) + , fHistSECVTX(0x0) +{ + // + // Default constructor + // + DefineInput(0, TChain::Class()); + DefineOutput(0, TH1I::Class()); + DefineOutput(1, TList::Class()); + DefineOutput(2, TList::Class()); + + // Initialize cuts + fPID = new AliHFEpid; + + SetHasMCData(); +} + +//____________________________________________________________ +AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name): + AliAnalysisTask(name, "") + , fQAlevel(0) + , fPIDdetectors("") + , fPIDstrategy(0) , fESD(0x0) , fMC(0x0) , fCFM(0x0) @@ -104,6 +141,7 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref): AliAnalysisTask(ref) , fQAlevel(ref.fQAlevel) , fPIDdetectors(ref.fPIDdetectors) + , fPIDstrategy(ref.fPIDstrategy) , fESD(ref.fESD) , fMC(ref.fMC) , fCFM(ref.fCFM) @@ -134,6 +172,7 @@ AliAnalysisTaskHFE &AliAnalysisTaskHFE::operator=(const AliAnalysisTaskHFE &ref) AliAnalysisTask::operator=(ref); fQAlevel = ref.fQAlevel; fPIDdetectors = ref.fPIDdetectors; + fPIDstrategy = ref.fPIDstrategy; fESD = ref.fESD; fMC = ref.fMC; fCFM = ref.fCFM; @@ -233,14 +272,14 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){ // Temporary fix: Initialize particle cuts with 0x0 for(Int_t istep = 0; istep < fCFM->GetParticleContainer()->GetNStep(); istep++) fCFM->SetParticleCutsList(istep, 0x0); - if(fCuts){ - fCuts->Initialize(fCFM); - if(fCuts->IsInDebugMode()){ - fQA->Add(fCuts->GetQAhistograms()); - } - } else { - AliError("Cuts not available. Output will be meaningless"); + if(!fCuts){ + AliWarning("Cuts not available. Default cuts will be used"); + fCuts = new AliHFEcuts; + fCuts->CreateStandardCuts(); } + fCuts->Initialize(fCFM); + if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms()); + // add output objects to the List fOutput->AddAt(fCFM->GetParticleContainer(), 0); fOutput->AddAt(fCorrelation, 1); @@ -255,8 +294,11 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){ fQA->Add(fPID->GetQAhistograms()); } fPID->SetHasMCData(HasMCData()); - if(!fPIDdetectors.Length()) AddPIDdetector("TPC"); - fPID->InitializePID(fPIDdetectors.Data()); // Only restrictions to TPC allowed + if(!fPIDdetectors.Length() && ! fPIDstrategy) AddPIDdetector("TPC"); + if(fPIDstrategy) + fPID->InitializePID(Form("Strategy%d", fPIDstrategy)); + else + fPID->InitializePID(fPIDdetectors.Data()); // Only restrictions to TPC allowed // mcQA---------------------------------- if (IsQAOn(kMCqa)) { @@ -268,8 +310,12 @@ void AliAnalysisTaskHFE::CreateOutputObjects(){ fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_"); // create histograms for beauty fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_"); // create histograms for charm fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_"); // create histograms for beauty - fMCQA->CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm - fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty + fMCQA->CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_"); // create histograms for charm + fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_"); // create histograms for beauty + fMCQA->CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_"); // create histograms for charm + fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_"); // create histograms for beauty + fMCQA->CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_"); // create histograms for charm + fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_"); // create histograms for beauty TIter next(gDirectory->GetList()); TObject *obj; int counter = 0; @@ -326,7 +372,6 @@ void AliAnalysisTaskHFE::Exec(Option_t *){ AliError("HFE cuts not available"); return; } - fCFM->SetEventInfo(fMC); //fCFM->CheckEventCuts(AliCFManager::kEvtGenCuts, fMC); Double_t container[6]; @@ -378,6 +423,7 @@ void AliAnalysisTaskHFE::Exec(Option_t *){ // // Loop MC // + fCuts->SetEventInfo(fMC); for(Int_t imc = fMC->GetNumberOfTracks(); imc--;){ mctrack = dynamic_cast(fMC->GetTrack(imc)); @@ -411,6 +457,7 @@ void AliAnalysisTaskHFE::Exec(Option_t *){ Bool_t signal = kTRUE; + fCuts->SetEventInfo(fESD); for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){ track = fESD->GetTrack(itrack); @@ -432,6 +479,7 @@ void AliAnalysisTaskHFE::Exec(Option_t *){ // Check if it is electrons near the vertex if(!(mctrack = dynamic_cast(fMC->GetTrack(TMath::Abs(track->GetLabel()))))) continue; + TParticle* mctrack4QA = fMC->Stack()->Particle(TMath::Abs(track->GetLabel())); container[3] = mctrack->Pt(); container[4] = mctrack->Eta(); @@ -492,6 +540,12 @@ void AliAnalysisTaskHFE::Exec(Option_t *){ ((THnSparseF *)fCorrelation->At(0))->Fill(container); } + if (IsQAOn(kMCqa)) { + // mc qa for after the reconstruction cuts + AliDebug(2, "Running MC QA"); + fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 3); // charm + fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 3); // beauty + } // track accepted, do PID AliHFEpidObject hfetrack; @@ -501,6 +555,12 @@ void AliAnalysisTaskHFE::Exec(Option_t *){ if(!fPID->IsSelected(&hfetrack)) continue; nElectronCandidates++; + if (IsQAOn(kMCqa)) { + // mc qa for after the reconstruction and pid cuts + AliDebug(2, "Running MC QA"); + fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kCharm, AliHFEmcQA::kElectronPDG, 4); // charm + fMCQA->GetDecayedKine(mctrack4QA, AliHFEmcQA::kBeauty, AliHFEmcQA::kElectronPDG, 4); // beauty + } // Fill Containers if(signal) { @@ -907,15 +967,7 @@ void AliAnalysisTaskHFE::PrintStatus(){ printf("\n"); printf("\tQA: \n"); printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" : "NO"); - if(fCuts) { - if(fCuts->IsInDebugMode()) { - printf("\t\tCUTS: YES\n"); - } else { - printf("\t\tCUTS: NO\n"); - } - } else { - printf("\t\tCUTS: NO\n"); - } + printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsInDebugMode()) ? "YES" : "NO"); printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO"); printf("\n"); } @@ -1022,7 +1074,7 @@ Float_t AliAnalysisTaskHFE::GetRapidity(TParticle *part) // return rapidity Float_t rapidity; - if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999; + if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999; else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); return rapidity; } diff --git a/PWG3/hfe/AliAnalysisTaskHFE.h b/PWG3/hfe/AliAnalysisTaskHFE.h index 5433ad44cdc..f6d8c725196 100644 --- a/PWG3/hfe/AliAnalysisTaskHFE.h +++ b/PWG3/hfe/AliAnalysisTaskHFE.h @@ -48,6 +48,7 @@ class AliAnalysisTaskHFE : public AliAnalysisTask{ kMCqa =1 }; AliAnalysisTaskHFE(); + AliAnalysisTaskHFE(const char * name); AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref); AliAnalysisTaskHFE& operator=(const AliAnalysisTaskHFE &ref); virtual ~AliAnalysisTaskHFE(); @@ -72,6 +73,7 @@ class AliAnalysisTaskHFE : public AliAnalysisTask{ void SetSecVtxOn(Bool_t option = kTRUE) { SetBit(kIsSecVtxOn, option); }; void SetRunPostProcess(Bool_t option = kTRUE) { SetBit(kIsRunningPostProcess, option); }; void SetPIDdetectors(Char_t *detectors){ fPIDdetectors = detectors; } + void SetPIDStrategy(UInt_t strategy) { fPIDstrategy = strategy; } void AddPIDdetector(Char_t *detector); void PrintStatus(); Float_t GetRapidity(TParticle *part); @@ -100,6 +102,7 @@ class AliAnalysisTaskHFE : public AliAnalysisTask{ ULong_t fQAlevel; // QA level TString fPIDdetectors; // Detectors for Particle Identification + UInt_t fPIDstrategy; // PID Strategy AliESDEvent *fESD; //! The ESD Event AliMCEvent *fMC; //! The MC Event AliCFManager *fCFM; //! Correction Framework Manager diff --git a/PWG3/hfe/AliHFEcuts.cxx b/PWG3/hfe/AliHFEcuts.cxx index 734f25651f1..1c7c8a04bb0 100644 --- a/PWG3/hfe/AliHFEcuts.cxx +++ b/PWG3/hfe/AliHFEcuts.cxx @@ -31,6 +31,8 @@ #include "AliCFAcceptanceCuts.h" #include "AliCFCutBase.h" +#include "AliCFEventGenCuts.h" +#include "AliCFEventRecCuts.h" #include "AliCFManager.h" #include "AliCFParticleGenCuts.h" #include "AliCFTrackIsPrimaryCuts.h" @@ -49,6 +51,7 @@ AliHFEcuts::AliHFEcuts(): fMinClustersTPC(0), fMinTrackletsTRD(0), fCutITSPixel(0), + fCheckITSLayerStatus(kTRUE), fMaxChi2clusterTPC(0.), fMinClusterRatioTPC(0.), fSigmaToVtx(0.), @@ -73,6 +76,7 @@ AliHFEcuts::AliHFEcuts(const AliHFEcuts &c): fMinClustersTPC(c.fMinClustersTPC), fMinTrackletsTRD(c.fMinTrackletsTRD), fCutITSPixel(c.fCutITSPixel), + fCheckITSLayerStatus(c.fCheckITSLayerStatus), fMaxChi2clusterTPC(c.fMaxChi2clusterTPC), fMinClusterRatioTPC(c.fMinClusterRatioTPC), fSigmaToVtx(c.fSigmaToVtx), @@ -120,8 +124,13 @@ void AliHFEcuts::Initialize(AliCFManager *cfm){ SetHFElectronITSCuts(); SetHFElectronTRDCuts(); + // Connect the event cuts + /*SetEventCutList(kEventStepGenerated); + SetEventCutList(kEventStepReconstructed); + cfm->SetEventCutsList(kEventStepGenerated, dynamic_cast(fCutList->FindObject("fEvGenCuts"))); + cfm->SetEventCutsList(kEventStepReconstructed, dynamic_cast(fCutList->FindObject("fEvRecCuts")));*/ - // Connect the cuts + // Connect the particle cuts cfm->SetParticleCutsList(kStepMCGenerated, dynamic_cast(fCutList->FindObject("fPartGenCuts"))); cfm->SetParticleCutsList(kStepMCInAcceptance, dynamic_cast(fCutList->FindObject("fPartAccCuts"))); cfm->SetParticleCutsList(kStepRecKineITSTPC, dynamic_cast(fCutList->FindObject("fPartRecKineITSTPCCuts"))); @@ -143,6 +152,61 @@ void AliHFEcuts::Initialize(){ } +//__________________________________________________________________ +void AliHFEcuts::SetEventInfo(TObject *ev){ + // + // Safe wrapper for AliCFManager::SetEventInfo that does a preselection of + // cut objects that need a special event info according to the object type + // in the argument. + // Gets rid of a few run time messages + // + if(!TString(ev->IsA()->GetName()).CompareTo("AliMCEvent")){ + TObjArray *genCuts = dynamic_cast(fCutList->FindObject("fPartGenCuts")); + if(genCuts){ + AliCFParticleGenCuts *myGenCut = dynamic_cast(genCuts->FindObject("fCutsGenMC")); + if(myGenCut) myGenCut->SetEvtInfo(ev); + } + } else if(!TString(ev->IsA()->GetName()).CompareTo("AliESDEvent")){ + TObjArray *primCuts = dynamic_cast(fCutList->FindObject("fPartPrimCuts")); + if(primCuts){ + AliCFTrackIsPrimaryCuts *myPrimCut = dynamic_cast(primCuts->FindObject("fCutsPrimaryCuts")); + if(myPrimCut) myPrimCut->SetEvtInfo(ev); + } + } +} + +//__________________________________________________________________ +void AliHFEcuts::SetEventCutList(Int_t istep){ + // + // Cuts for Event Selection + // + TObjArray *arr = new TObjArray; + if(istep == kEventStepGenerated){ + AliCFEventGenCuts *evGenCuts = new AliCFEventGenCuts("fCutsEvGen", "Event Generated cuts"); + evGenCuts->SetNTracksCut(1); + evGenCuts->SetRequireVtxCuts(kTRUE); + evGenCuts->SetVertexXCut(-1, 1); + evGenCuts->SetVertexYCut(-1, 1); + evGenCuts->SetVertexZCut(-30, 30); + + arr->SetName("fEvGenCuts"); + arr->AddLast(evGenCuts); + } else { + AliCFEventRecCuts *evRecCuts = new AliCFEventRecCuts("fCutsEvRec", "Event Reconstructed cuts"); + evRecCuts->SetUseSPDVertex(); + evRecCuts->SetNTracksCut(1); + evRecCuts->SetRequireVtxCuts(kTRUE); + evRecCuts->SetVertexXCut(-1, 1); + evRecCuts->SetVertexYCut(-1, 1); + evRecCuts->SetVertexZCut(-30, 30); + + + arr->SetName("fEvRecCuts"); + arr->AddLast(evRecCuts); + } + fCutList->AddLast(arr); +} + //__________________________________________________________________ void AliHFEcuts::SetParticleGenCutList(){ // @@ -284,6 +348,7 @@ void AliHFEcuts::SetHFElectronITSCuts(){ AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupPixels","Extra cuts from the HFE group"); if(IsRequireITSpixel()){ hfecuts->SetRequireITSpixel(AliHFEextraCuts::ITSPixel_t(fCutITSPixel)); + hfecuts->SetCheckITSstatus(fCheckITSLayerStatus); } if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA); diff --git a/PWG3/hfe/AliHFEcuts.h b/PWG3/hfe/AliHFEcuts.h index 3f3d66134cf..af188d66371 100644 --- a/PWG3/hfe/AliHFEcuts.h +++ b/PWG3/hfe/AliHFEcuts.h @@ -53,6 +53,10 @@ class AliHFEcuts : public TObject{ kStepHFEcutsITS = 5, kStepHFEcutsTRD = 6 } CutStep_t; + typedef enum{ + kEventStepGenerated = 0, + kEventStepReconstructed = 1 + } EventCutStep_t; enum{ kNcutSteps = 7, kNcutESDSteps = 4 @@ -71,6 +75,7 @@ class AliHFEcuts : public TObject{ TList *GetQAhistograms() const { return fHistQA; } void SetDebugMode(); + void SetEventInfo(TObject *ev); void UnsetDebugMode() { SetBit(kDebugMode, kFALSE); } Bool_t IsInDebugMode() const { return TestBit(kDebugMode); }; @@ -84,6 +89,7 @@ class AliHFEcuts : public TObject{ // Setters inline void SetCutITSpixel(UChar_t cut); + void SetCheckITSLayerStatus(Bool_t checkITSLayerStatus) { fCheckITSLayerStatus = checkITSLayerStatus; } void SetMinNClustersTPC(UChar_t minClustersTPC) { fMinClustersTPC = minClustersTPC; } void SetMinNTrackletsTRD(UChar_t minNtrackletsTRD) { fMinTrackletsTRD = minNtrackletsTRD; } void SetMaxChi2perClusterTPC(Double_t chi2) { fMaxChi2clusterTPC = chi2; }; @@ -112,6 +118,7 @@ class AliHFEcuts : public TObject{ void SetRecPrimaryCutList(); void SetHFElectronITSCuts(); void SetHFElectronTRDCuts(); + void SetEventCutList(Int_t istep); ULong64_t fRequirements; // Bitmap for requirements Double_t fDCAtoVtx[2]; // DCA to Vertex @@ -120,6 +127,7 @@ class AliHFEcuts : public TObject{ UChar_t fMinClustersTPC; // Min.Number of TPC clusters UChar_t fMinTrackletsTRD; // Min. Number of TRD tracklets UChar_t fCutITSPixel; // Cut on ITS pixel + Bool_t fCheckITSLayerStatus; // Check ITS layer status Double_t fMaxChi2clusterTPC; // Max Chi2 per TPC cluster Double_t fMinClusterRatioTPC; // Min. Ratio findable / found TPC clusters Double_t fSigmaToVtx; // Sigma To Vertex @@ -172,7 +180,7 @@ void AliHFEcuts::CreateStandardCuts(){ fProdVtx[2] = -3; fProdVtx[3] = 3; SetRequireDCAToVertex(); - fDCAtoVtx[0] = 4.; + fDCAtoVtx[0] = 2.; fDCAtoVtx[1] = 10.; fMinClustersTPC = 50; fMinTrackletsTRD = 0; diff --git a/PWG3/hfe/AliHFEmcQA.cxx b/PWG3/hfe/AliHFEmcQA.cxx index eac28fd64c6..4997af35b06 100644 --- a/PWG3/hfe/AliHFEmcQA.cxx +++ b/PWG3/hfe/AliHFEmcQA.cxx @@ -48,6 +48,7 @@ #include #include #include +#include #include "AliHFEmcQA.h" @@ -61,6 +62,7 @@ ClassImp(AliHFEmcQA) //_______________________________________________________________________________________________ AliHFEmcQA::AliHFEmcQA() : fStack(0x0) + ,fMCArray(0x0) ,fNparents(0) { // Default constructor @@ -71,6 +73,7 @@ AliHFEmcQA::AliHFEmcQA() : AliHFEmcQA::AliHFEmcQA(const AliHFEmcQA&p): TObject(p) ,fStack(0x0) + ,fMCArray(0x0) ,fNparents(p.fNparents) { // Copy constructor @@ -457,6 +460,7 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde return; } Int_t iq = kquark - kCharm; + Bool_t IsFinalOpenCharm = kFALSE; /* if (iTrack < 0) { @@ -488,6 +492,13 @@ void AliHFEmcQA::GetDecayedKine(TParticle* mcpart, const Int_t kquark, Int_t kde Bool_t isMotherDirectCharm = kFALSE; if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { + for (Int_t i=0; iGetPdgCode()) != kdecayed ) return; + + // mother + Int_t iLabel = mcpart->GetMother(); + if (iLabel<0){ + AliDebug(1, "Stack label is negative, return\n"); + return; + } + + AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel); + AliAODMCParticle *partMotherCopy = partMother; + Int_t maPdgcode = partMother->GetPdgCode(); + Int_t maPdgcodeCopy = maPdgcode; + + Bool_t isMotherDirectCharm = kFALSE; + if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { + + for (Int_t i=0; iGetMother(); + if (jLabel == -1){ + isMotherDirectCharm = kTRUE; + break; // if there is no ancester + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return; + } + + // if there is an ancester + AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t j=0; jFill(mcpart->GetPdgCode()); + fHist[iq][kElectron2nd][icut].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron2nd][icut].fY->Fill(GetRapidity(mcpart)); + fHist[iq][kElectron2nd][icut].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][kDeHadron][icut].fPdgCode->Fill(grandMaPDG); + fHist[iq][kDeHadron][icut].fPt->Fill(grandMa->Pt()); + fHist[iq][kDeHadron][icut].fY->Fill(GetRapidity(grandMa)); + fHist[iq][kDeHadron][icut].fEta->Fill(grandMa->Eta()); + + return; + } + } + + partMother = grandMa; + } // end of iteration + } // end of if + if ((isMotherDirectCharm == kTRUE && kquark == kCharm) || kquark == kBeauty) { + for (Int_t i=0; iFill(mcpart->GetPdgCode()); + fHist[iq][kElectron][icut].fPt->Fill(mcpart->Pt()); + fHist[iq][kElectron][icut].fY->Fill(GetRapidity(mcpart)); + fHist[iq][kElectron][icut].fEta->Fill(mcpart->Eta()); + + // fill mother hadron kinematics + fHist[iq][keHadron][icut].fPdgCode->Fill(maPdgcodeCopy); + fHist[iq][keHadron][icut].fPt->Fill(partMotherCopy->Pt()); + fHist[iq][keHadron][icut].fY->Fill(GetRapidity(partMotherCopy)); + fHist[iq][keHadron][icut].fEta->Fill(partMotherCopy->Eta()); + + } + } + } // end of if + +} //__________________________________________ void AliHFEmcQA::IdentifyMother(Int_t motherlabel, Int_t &motherpdg, Int_t &grandmotherlabel) @@ -666,7 +774,172 @@ Float_t AliHFEmcQA::GetRapidity(TParticle *part) // return rapidity Float_t rapidity; - if(part->Energy() - part->Pz() == 0 || part->Energy() + part->Pz() == 0) rapidity=-999; + if(!((part->Energy() - part->Pz())*(part->Energy() + part->Pz())>0)) rapidity=-999; else rapidity = 0.5*(TMath::Log((part->Energy()+part->Pz()) / (part->Energy()-part->Pz()))); return rapidity; } + +//__________________________________________ +Float_t AliHFEmcQA::GetRapidity(AliAODMCParticle *part) +{ + // return rapidity + + Float_t rapidity; + if(!((part->E() - part->Pz())*(part->E() + part->Pz())>0)) rapidity=-999; + else rapidity = 0.5*(TMath::Log((part->E()+part->Pz()) / (part->E()-part->Pz()))); + return rapidity; +} + +//__________________________________________ +Int_t AliHFEmcQA::GetElectronSource(AliAODMCParticle *mcpart) +{ + // decay electron's origin + + if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1; + + Int_t origin = -1; + Bool_t IsFinalOpenCharm = kFALSE; + + // mother + Int_t iLabel = mcpart->GetMother(); + if (iLabel<0){ + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + AliAODMCParticle *partMother = (AliAODMCParticle*)fMCArray->At(iLabel); + Int_t maPdgcode = partMother->GetPdgCode(); + + // if the mother is charmed hadron + if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { + + for (Int_t i=0; iGetMother(); + if (jLabel == -1){ + origin = kDirectCharm; + return origin; + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + // if there is an ancester + AliAODMCParticle* grandMa = (AliAODMCParticle*)fMCArray->At(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t i=0; iGetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1; + + Int_t origin = -1; + Bool_t IsFinalOpenCharm = kFALSE; + + Int_t iLabel = mcpart->GetFirstMother(); + if (iLabel<0){ + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + TParticle *partMother = fStack->Particle(iLabel); + Int_t maPdgcode = partMother->GetPdgCode(); + + // if the mother is charmed hadron + if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) { + + for (Int_t i=0; iGetFirstMother(); + if (jLabel == -1){ + origin = kDirectCharm; + return origin; + } + if (jLabel < 0){ // safety protection + AliDebug(1, "Stack label is negative, return\n"); + return -1; + } + + // if there is an ancester + TParticle* grandMa = fStack->Particle(jLabel); + Int_t grandMaPDG = grandMa->GetPdgCode(); + + for (Int_t i=0; iMakeIterator(); TObjString *det = 0x0; while((det = dynamic_cast(detIterator->Next()))){ if(det->String().CompareTo("TPC") == 0){ AliInfo("Doing TPC PID"); fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPC PID"); - (dynamic_cast(fDetectorPID[kTPCpid]))->SetAsymmetricTPCsigmaCut(2., 10., 0., 4.); SETBIT(fEnabledDetectors, kTPCpid); } else if(det->String().CompareTo("TRD") == 0){ fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID"); @@ -181,6 +200,18 @@ Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){ // MC Event return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11); } + if(fPIDstrategy > 0 && fPIDstrategy < 6){ + Int_t pid = 0; + switch(fPIDstrategy){ + case 1: pid = IdentifyStrategy1(track); break; + case 2: pid = IdentifyStrategy2(track); break; + case 3: pid = IdentifyStrategy3(track); break; + case 4: pid = IdentifyStrategy4(track); break; + case 5: pid = IdentifyStrategy5(track); break; + default: break; + } + return pid; + } if(TESTBIT(fEnabledDetectors, kTPCpid)){ if(IsQAOn() && fDebugLevel > 1){ AliInfo("Filling QA plots"); @@ -220,6 +251,7 @@ Bool_t AliHFEpid::MakePidTpcTrd(AliHFEpidObject *track){ // momentum slices // if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return kFALSE; //AOD based detector PID combination not yet implemented + AliDebug(1, "Analysis Type OK, do PID"); AliESDtrack *esdTrack = dynamic_cast(track->fRecTrack); AliHFEpidTRD *trdPid = dynamic_cast(fDetectorPID[kTRDpid]); Int_t pdg = TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)); @@ -232,18 +264,25 @@ Bool_t AliHFEpid::MakePidTpcTrd(AliHFEpidObject *track){ case 2212: pid = AliPID::kProton; break; default: pid = -1; }; + Double_t content[10]; + content[0] = pid; + content[1] = esdTrack->P(); + content[2] = esdTrack->GetTPCsignal(); + content[3] = trdPid->GetTRDSignalV1(esdTrack, pid); if(IsQAOn() && fDebugLevel > 1){ - Double_t content[10]; - content[0] = pid; - content[1] = esdTrack->P(); - content[2] = esdTrack->GetTPCsignal(); - content[3] = trdPid->GetTRDSignalV1(esdTrack, pid); content[4] = trdPid->GetTRDSignalV2(esdTrack, pid); AliDebug(1, Form("Momentum: %f, TRD Signal: Method 1[%f], Method 2[%f]", content[1], content[3], content[4])); (dynamic_cast(fQAlist->At(kTRDSignal)))->Fill(content); } - Int_t trdDecision = 0; // TRD decision 0 means outside the allowed momentum region - return ((trdDecision = trdPid->IsSelected(track)) == 11 || trdDecision == 0) && fDetectorPID[kTPCpid]->IsSelected(track) == 11; + if(content[1] > 2){ // perform combined + AliDebug(1, "Momentum bigger 2 GeV/c, doing combined PID"); + if(content[2] > 65 && content[3] > 500) return kTRUE; + else return kFALSE; + } + else { + AliDebug(1, "Momentum smaller 2GeV/c, doing TPC alone PID"); + return fDetectorPID[kTPCpid]->IsSelected(track) == 11; + } } //____________________________________________________________ @@ -327,3 +366,95 @@ void AliHFEpid::SetQAOn(){ } } +//____________________________________________________________ +void AliHFEpid::InitStrategy1(){ + // + // TPC alone, 3-sigma cut + // + AliHFEpidTPC *pid = new AliHFEpidTPC("strat1TPCpid"); + pid->SetTPCnSigma(3); + Bool_t status = pid->InitializePID(); + if(IsQAOn() && status) pid->SetQAOn(fQAlist); + if(HasMCData() && status) pid->SetHasMCData(); + fDetectorPID[kTPCpid] = pid; +} + +//____________________________________________________________ +void AliHFEpid::InitStrategy2(){ + // + // TPC alone, symmetric 3 sigma cut and asymmetric sigma cut in the momentum region between 2GeV/c and 10 GeV/c and sigma between -1 and 100 + // + AliHFEpidTPC *pid = new AliHFEpidTPC("strat2TPCpid"); + pid->SetTPCnSigma(3); + pid->SetAsymmetricTPCsigmaCut(2., 10., 0., 4.); + Bool_t status = pid->InitializePID(); + if(IsQAOn() && status) pid->SetQAOn(fQAlist); + if(HasMCData() && status) pid->SetHasMCData(); + fDetectorPID[kTPCpid] = pid; +} + +//____________________________________________________________ +void AliHFEpid::InitStrategy3(){ + // + // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection + // + AliHFEpidTPC *pid = new AliHFEpidTPC("strat3TPCpid"); + pid->SetTPCnSigma(3); + pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 3.); + Bool_t status = pid->InitializePID(); + if(IsQAOn() && status) pid->SetQAOn(fQAlist); + if(HasMCData() && status) pid->SetHasMCData(); + fDetectorPID[kTPCpid] = pid; +} + +//____________________________________________________________ +void AliHFEpid::InitStrategy4(){ + // + // TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c + // + InitStrategy1(); + AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat4TRDpid"); + Bool_t status = trdpid->InitializePID(); + if(IsQAOn() && status) trdpid->SetQAOn(fQAlist); + if(HasMCData() && status) trdpid->SetHasMCData(); + fDetectorPID[kTRDpid] = trdpid; +} + +//____________________________________________________________ +void AliHFEpid::InitStrategy5(){ + // + // TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c + // + InitStrategy1(); + AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat5TRDpid"); + Bool_t status = trdpid->InitializePID(); + if(IsQAOn() && status) trdpid->SetQAOn(fQAlist); + if(HasMCData() && status) trdpid->SetHasMCData(); + fDetectorPID[kTRDpid] = trdpid; +} + +//____________________________________________________________ +Bool_t AliHFEpid::IdentifyStrategy1(AliHFEpidObject *track){ + return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11; +} + +//____________________________________________________________ +Bool_t AliHFEpid::IdentifyStrategy2(AliHFEpidObject *track){ + return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11; +} + +//____________________________________________________________ +Bool_t AliHFEpid::IdentifyStrategy3(AliHFEpidObject *track){ + return TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11; +} + +//____________________________________________________________ +Bool_t AliHFEpid::IdentifyStrategy4(AliHFEpidObject *track){ + Int_t trdpid = TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track)); + return (trdpid == 0 || trdpid == 11) && (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11); +} + +//____________________________________________________________ +Bool_t AliHFEpid::IdentifyStrategy5(AliHFEpidObject *track){ + return MakePidTpcTrd(track); +} diff --git a/PWG3/hfe/AliHFEpid.h b/PWG3/hfe/AliHFEpid.h index 091de8f01f2..881abcb84fb 100644 --- a/PWG3/hfe/AliHFEpid.h +++ b/PWG3/hfe/AliHFEpid.h @@ -56,7 +56,7 @@ class AliHFEpid : public TObject{ AliHFEpid &operator=(const AliHFEpid &c); ~AliHFEpid(); - Bool_t InitializePID(TString detectors); + Bool_t InitializePID(TString argument); Bool_t IsSelected(AliHFEpidObject *track); Bool_t IsQAOn() const { return TestBit(kIsQAOn); }; @@ -70,9 +70,22 @@ class AliHFEpid : public TObject{ Bool_t MakePidTpcTof(AliHFEpidObject *track); Bool_t MakePidTpcTrd(AliHFEpidObject *track); void MakePlotsItsTpc(AliHFEpidObject *track); + + // Stratgies + void InitStrategy1(); + void InitStrategy2(); + void InitStrategy3(); + void InitStrategy4(); + void InitStrategy5(); + Bool_t IdentifyStrategy1(AliHFEpidObject *track); + Bool_t IdentifyStrategy2(AliHFEpidObject *track); + Bool_t IdentifyStrategy3(AliHFEpidObject *track); + Bool_t IdentifyStrategy4(AliHFEpidObject *track); + Bool_t IdentifyStrategy5(AliHFEpidObject *track); private: AliHFEpidBase *fDetectorPID[kNdetectorPID]; //! Detector PID classes UInt_t fEnabledDetectors; // Enabled Detectors + UInt_t fPIDstrategy; // PID Strategy TList *fQAlist; //! QA histograms Int_t fDebugLevel; // Debug Level diff --git a/PWG3/hfe/AliHFEpidTPC.cxx b/PWG3/hfe/AliHFEpidTPC.cxx index 13b51165f49..dbbbb3935c6 100644 --- a/PWG3/hfe/AliHFEpidTPC.cxx +++ b/PWG3/hfe/AliHFEpidTPC.cxx @@ -53,6 +53,7 @@ AliHFEpidTPC::AliHFEpidTPC(const char* name) : AliHFEpidBase(name) , fLineCrossingsEnabled(0) , fNsigmaTPC(3) + , fRejectionEnabled(0) , fPID(NULL) , fPIDtpcESD(NULL) , fQAList(NULL) @@ -72,6 +73,7 @@ AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) : AliHFEpidBase("") , fLineCrossingsEnabled(0) , fNsigmaTPC(2) + , fRejectionEnabled(0) , fPID(NULL) , fPIDtpcESD(NULL) , fQAList(NULL) @@ -103,6 +105,7 @@ void AliHFEpidTPC::Copy(TObject &o) const{ target.fLineCrossingsEnabled = fLineCrossingsEnabled; target.fNsigmaTPC = fNsigmaTPC; + target.fRejectionEnabled = fRejectionEnabled; target.fPID = new AliPID(*fPID); target.fPIDtpcESD = new AliTPCpidESD(*fPIDtpcESD); target.fQAList = dynamic_cast(fQAList->Clone()); @@ -176,17 +179,26 @@ Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){ } } if(isLineCrossing) return 0; + + // Check particle rejection + if(HasParticleRejection()){ + Int_t reject = Reject(esdTrack); + AliDebug(1, Form("PID code from Rejection: %d", reject)); + if(reject != 0) return reject; + } // Check whether distance from the electron line is smaller than n-sigma // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut Float_t nsigma = fPIDtpcESD->GetNumberOfSigmas(esdTrack, AliPID::kElectron); Float_t p = 0.; + Int_t pdg = 0; if(HasAsymmetricSigmaCut() && (p = esdTrack->P()) >= fPAsigCut[0] && p <= fPAsigCut[1]){ - if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) return 11; + if(nsigma >= fNAsigmaTPC[0] && nsigma <= fNAsigmaTPC[1]) pdg = 11; } else { - if(TMath::Abs(nsigma) < fNsigmaTPC ) return 11; + if(TMath::Abs(nsigma) < fNsigmaTPC ) pdg = 11; } - return 0; + if(IsQAon() && pdg != 0) (dynamic_cast(fQAList->At(kHistTPCselected)))->Fill(esdTrack->GetInnerParam() ? esdTrack->GetInnerParam()->P() : esdTrack->P(), esdTrack->GetTPCsignal()); + return pdg; } //___________________________________________________________________ @@ -195,6 +207,25 @@ Int_t AliHFEpidTPC::MakePIDaod(AliAODTrack * /*aodTrack*/, AliAODMCParticle * /* return 0; } +//___________________________________________________________________ +Int_t AliHFEpidTPC::Reject(AliESDtrack *track){ + // + // reject particles based on asymmetric sigma cut + // + Int_t pdc[AliPID::kSPECIES] = {11,13,211,321,2212}; + Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P(); + for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){ + if(!TESTBIT(fRejectionEnabled, ispec)) continue; + AliDebug(1, Form("Particle Rejection enabled for species %d", ispec)); + // Particle rejection enabled + if(p < fRejection[4*ispec] || p > fRejection[4*ispec+2]) continue; + Double_t sigma = fPIDtpcESD->GetNumberOfSigmas(track, static_cast(ispec)); + AliDebug(1, Form("Sigma %f, min %f, max %f", sigma, fRejection[4*ispec + 1], fRejection[4*ispec+3])); + if(sigma >= fRejection[4*ispec+1] && sigma <= fRejection[4*ispec+3]) return pdc[ispec] * track->Charge(); + } + return 0; +} + //___________________________________________________________________ void AliHFEpidTPC::AddTPCdEdxLineCrossing(Int_t species, Double_t sigma){ // @@ -317,13 +348,14 @@ void AliHFEpidTPC::AddQAhistograms(TList *qaList){ fQAList = new TList; fQAList->SetName("fTPCqaHistos"); - fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 60, 0, 600), kHistTPCelectron); - fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 60, 0, 600), kHistTPCmuon); - fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 60, 0, 600), kHistTPCpion); - fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 60, 0, 600), kHistTPCkaon); - fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 60, 0, 600), kHistTPCproton); - fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 60, 0, 600), kHistTPCothers); - fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 60, 0, 600), kHistTPCall); + fQAList->AddAt(new TH2I("fHistTPCelectron","TPC signal for Electrons", 200, 0, 20, 200, 0, 200), kHistTPCelectron); + fQAList->AddAt(new TH2I("fHistTPCmuon","TPC signal for Muons", 200, 0, 20, 200, 0, 200), kHistTPCmuon); + fQAList->AddAt(new TH2I("fHistTPCpion","TPC signal for Pions", 200, 0, 20, 200, 0, 200), kHistTPCpion); + fQAList->AddAt(new TH2I("fHistTPCkaon","TPC signal for Kaons", 200, 0, 20, 200, 0, 200), kHistTPCkaon); + fQAList->AddAt(new TH2I("fHistTPCproton","TPC signal for Protons", 200, 0, 20, 200, 0, 200), kHistTPCproton); + fQAList->AddAt(new TH2I("fHistTPCothers","TPC signal for other species", 200, 0, 20, 200, 0, 200), kHistTPCothers); + fQAList->AddAt(new TH2I("fHistTPCall","TPC signal for all species", 200, 0, 20, 200, 0, 200), kHistTPCall); + fQAList->AddAt(new TH2I("fHistTPCselected","TPC signal for all selected particles", 200, 0, 20, 200, 0, 200), kHistTPCselected); fQAList->AddAt(new TH2F("fHistTPCprobEl","TPC likelihood for electrons to be an electron vs. p", 200, 0.,20.,200,0.,1.), kHistTPCprobEl); fQAList->AddAt(new TH2F("fHistTPCprobPi","TPC likelihood for pions to be an electron vs. p", 200, 0.,20.,200, 0.,1.), kHistTPCprobPi); diff --git a/PWG3/hfe/AliHFEpidTPC.h b/PWG3/hfe/AliHFEpidTPC.h index 7f457b2d07d..2e82dd5b4ed 100644 --- a/PWG3/hfe/AliHFEpidTPC.h +++ b/PWG3/hfe/AliHFEpidTPC.h @@ -40,28 +40,30 @@ class AliHFEpidTPC : public AliHFEpidBase{ kHistTPCproton = 4, kHistTPCothers = 5, kHistTPCall = 6, - kHistTPCprobEl = 7, - kHistTPCprobPi = 8, - kHistTPCprobMu = 9, - kHistTPCprobKa = 10, - kHistTPCprobPro = 11, - kHistTPCprobOth = 12, - kHistTPCprobAll = 13, - kHistTPCsuppressPi = 14, - kHistTPCsuppressMu = 15, - kHistTPCsuppressKa = 16, - kHistTPCsuppressPro = 17, - kHistTPCenhanceElPi = 18, - kHistTPCenhanceElMu = 19, - kHistTPCenhanceElKa = 20, - kHistTPCenhanceElPro = 21, - kHistTPCElprobPi = 22, - kHistTPCElprobMu = 23, - kHistTPCElprobKa = 24, - kHistTPCElprobPro = 25 + kHistTPCselected = 7, + kHistTPCprobEl = 8, + kHistTPCprobPi = 9, + kHistTPCprobMu = 10, + kHistTPCprobKa = 11, + kHistTPCprobPro = 12, + kHistTPCprobOth = 13, + kHistTPCprobAll = 14, + kHistTPCsuppressPi = 15, + kHistTPCsuppressMu = 16, + kHistTPCsuppressKa = 17, + kHistTPCsuppressPro = 18, + kHistTPCenhanceElPi = 19, + kHistTPCenhanceElMu = 20, + kHistTPCenhanceElKa = 21, + kHistTPCenhanceElPro = 22, + kHistTPCElprobPi = 23, + kHistTPCElprobMu = 24, + kHistTPCElprobKa = 25, + kHistTPCElprobPro = 26 } QAHist_t; enum{ - kAsymmetricSigmaCut = BIT(20) + kAsymmetricSigmaCut = BIT(20), + kRejection = BIT(21) }; public: AliHFEpidTPC(const Char_t *name); @@ -75,8 +77,10 @@ class AliHFEpidTPC : public AliHFEpidBase{ void AddTPCdEdxLineCrossing(Int_t species, Double_t sigma); Bool_t HasAsymmetricSigmaCut() const { return TestBit(kAsymmetricSigmaCut);} + Bool_t HasParticleRejection() const { return TestBit(kRejection); } void SetTPCnSigma(Short_t nSigma) { fNsigmaTPC = nSigma; }; inline void SetAsymmetricTPCsigmaCut(Float_t pmin, Float_t pmax, Float_t sigmaMin, Float_t sigmaMax); + inline void SetRejectParticle(Int_t species, Float_t pmin, Float_t sigmaMin, Float_t pmax, Float_t sigmaMax); protected: void Copy(TObject &o) const; @@ -84,6 +88,7 @@ class AliHFEpidTPC : public AliHFEpidBase{ void FillTPChistograms(const AliESDtrack *track, const AliMCParticle *mctrack); Int_t MakePIDaod(AliAODTrack *aodTrack, AliAODMCParticle *mcTrack); Int_t MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mcTrack); + Int_t Reject(AliESDtrack *track); Double_t Likelihood(const AliESDtrack *track, Int_t species, Float_t rsig = 2.); Double_t Suppression(const AliESDtrack *track, Int_t species); @@ -93,6 +98,8 @@ class AliHFEpidTPC : public AliHFEpidBase{ Float_t fPAsigCut[2]; // Momentum region where to perform asymmetric sigma cut Float_t fNAsigmaTPC[2]; // Asymmetric TPC Sigma band Short_t fNsigmaTPC; // TPC sigma band + Float_t fRejection[4*AliPID::kSPECIES]; // All informations for Particle Rejection, order pmin, sigmin, pmax, sigmax + UChar_t fRejectionEnabled; // Bitmap for enabled particle rejection AliPID *fPID; //! PID Object AliTPCpidESD *fPIDtpcESD; //! TPC PID object TList *fQAList; //! QA histograms @@ -107,5 +114,15 @@ inline void AliHFEpidTPC::SetAsymmetricTPCsigmaCut(Float_t pmin, Float_t pmax, F fNAsigmaTPC[1] = sigmaMax; SetBit(kAsymmetricSigmaCut, kTRUE); } + +inline void AliHFEpidTPC::SetRejectParticle(Int_t species, Float_t pmin, Float_t sigmaMin, Float_t pmax, Float_t sigmaMax){ + if(species < 0 || species >= AliPID::kSPECIES) return; + fRejection[4*species] = pmin; + fRejection[4*species+1] = sigmaMin; + fRejection[4*species+2] = pmax; + fRejection[4*species+3] = sigmaMax; + SETBIT(fRejectionEnabled, species); + SetBit(kRejection, kTRUE); +} #endif diff --git a/PWG3/hfe/AliHFEpidTRD.cxx b/PWG3/hfe/AliHFEpidTRD.cxx index b3af3a7dcc1..005eabfaf3b 100644 --- a/PWG3/hfe/AliHFEpidTRD.cxx +++ b/PWG3/hfe/AliHFEpidTRD.cxx @@ -151,8 +151,15 @@ Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle * /*mcTrack* Double_t pidProbs[AliPID::kSPECIES]; esdTrack->GetTRDpid(pidProbs); - if(pidProbs[AliPID::kElectron] > GetTRDthresholds(0.91, p)) return 11; - return 0; + if(IsQAon())FillHistogramsLikelihood(0, p, pidProbs[AliPID::kElectron]); + Double_t threshold = GetTRDthresholds(0.91, p); + AliDebug(1, Form("Threshold: %f\n", threshold)); + if(IsQAon()) (dynamic_cast(fContainer->At(kHistTRDthresholds)))->Fill(p, threshold); + if(pidProbs[AliPID::kElectron] > threshold){ + if(IsQAon()) FillHistogramsLikelihood(1, p, pidProbs[AliPID::kElectron]); + return 11; + } + return 211; } //___________________________________________________________________ @@ -162,7 +169,8 @@ Double_t AliHFEpidTRD::GetTRDthresholds(Double_t electronEff, Double_t p){ // Double_t params[4]; GetParameters(electronEff, params); - return 1. - params[0] * p - params[1] * TMath::Exp(-params[2] * p) - params[3]; + Double_t threshold = 1. - params[0] - params[1] * p - params[2] * TMath::Exp(-params[3] * p); + return TMath::Max(TMath::Min(threshold, 0.99), 0.2); // truncate the threshold upperwards to 0.999 and lowerwards to 0.2 and exclude unphysical values } //___________________________________________________________________ @@ -203,7 +211,7 @@ void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){ // // return parameter set for the given efficiency bin // - Int_t effbin = static_cast((electronEff - 0.7)/0.3 * 6.); + Int_t effbin = static_cast((electronEff - 0.7)/0.05); memcpy(parameters, fThreshParams + effbin * 4, sizeof(Double_t) * 4); } @@ -322,6 +330,11 @@ void AliHFEpidTRD::AddQAhistograms(TList *l){ Double_t momentumBins[kMomentumBins]; for(Int_t ibin = 0; ibin < kMomentumBins; ibin++) momentumBins[ibin] = static_cast(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast(ibin))); + // Likelihood Histograms + fContainer->AddAt(new TH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeBefore); + fContainer->AddAt(new TH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeAfter); + fContainer->AddAt(new TH2F("fTRDthesholds", "TRD Electron thresholds", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDthresholds); + // Signal Histograms fContainer->AddAt(new TH2F("fTRDSigV1all", "TRD Signal (all particles, Method 1)", kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistTRDSigV1); fContainer->AddAt(new TH2F("fTRDSigV2all", "TRD Signal (all particles, Method 2)", kMomentumBins - 1, momentumBins, kSigBinsMeth2, kMinSig, kMaxSigMeth2), kHistTRDSigV2); for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){ @@ -331,3 +344,19 @@ void AliHFEpidTRD::AddQAhistograms(TList *l){ l->AddLast(fContainer); } +//___________________________________________________________________ +void AliHFEpidTRD::FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb){ + // + // Fill Likelihood Histogram before respectively after decision + // + TH2F *histo = NULL; + if(whenFilled) + histo = dynamic_cast(fContainer->At(kHistTRDlikeAfter)); + else + histo = dynamic_cast(fContainer->At(kHistTRDlikeBefore)); + if(!histo){ + AliError("QA histograms not found"); + return; + } + histo->Fill(p, elProb); +} diff --git a/PWG3/hfe/AliHFEpidTRD.h b/PWG3/hfe/AliHFEpidTRD.h index 54768f70f95..fe512802982 100644 --- a/PWG3/hfe/AliHFEpidTRD.h +++ b/PWG3/hfe/AliHFEpidTRD.h @@ -37,9 +37,12 @@ class AliHFEpidTRD : public AliHFEpidBase{ kThreshParams = 24 }; enum{ - kHistTRDSigV1 = 0, - kHistTRDSigV2 = 1, - kHistOverallSpecies = 2 + kHistTRDlikeBefore = 0, + kHistTRDlikeAfter = 1, + kHistTRDthresholds = 2, + kHistTRDSigV1 = 3, + kHistTRDSigV2 = 4, + kHistOverallSpecies = 5 }; AliHFEpidTRD(const Char_t *name); AliHFEpidTRD(const AliHFEpidTRD &ref); @@ -64,6 +67,7 @@ class AliHFEpidTRD : public AliHFEpidBase{ virtual void AddQAhistograms(TList *l); void GetParameters(Double_t electronEff, Double_t *parameters); + void FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb); void FillHistogramsTRDSignalV1(Double_t signal, Double_t p, Int_t species); void FillHistogramsTRDSignalV2(Double_t signal, Double_t p, Int_t species); private: diff --git a/PWG3/hfe/AliHFEsecVtx.cxx b/PWG3/hfe/AliHFEsecVtx.cxx index 208611ea917..cd632f72a2b 100644 --- a/PWG3/hfe/AliHFEsecVtx.cxx +++ b/PWG3/hfe/AliHFEsecVtx.cxx @@ -819,6 +819,7 @@ Int_t AliHFEsecVtx::PairOrigin(AliESDtrack* track1, AliESDtrack* track2) // caution : it can also return parton pdg code if it originated from same string or gluon spliting. // + if (track1->GetLabel()<0 || track2->GetLabel()<0) return 0; TParticle* part1 = fStack->Particle(TMath::Abs(track1->GetLabel())); TParticle* part2 = fStack->Particle(TMath::Abs(track2->GetLabel())); if (!(part1)) return 0; diff --git a/PWG3/hfe/runElectronTask.C b/PWG3/hfe/runElectronTask.C index 2982d3952aa..950fcd3df74 100644 --- a/PWG3/hfe/runElectronTask.C +++ b/PWG3/hfe/runElectronTask.C @@ -46,12 +46,12 @@ void runElectronTask(const char *treelist = 0x0){ hfecuts->CreateStandardCuts(); AliAnalysisTaskHFE *task = new AliAnalysisTaskHFE; task->SetHFECuts(hfecuts); - task->AddPIDdetector("TPC"); + task->SetPIDStrategy("Strategy4"); task->SetQAOn(AliAnalysisTaskHFE::kPIDqa); task->SetQAOn(AliAnalysisTaskHFE::kMCqa); task->SetSecVtxOn(); pidEffManager->AddTask(task); - task->ConnectInput(0, pidEffManager->CreateContainer("esdTree", TChain::Class(), AliAnalysisManager::kInputContainer)); + task->ConnectInput(0, pidEffManager->GetCommonInputContainer()); task->ConnectOutput(0, pidEffManager->CreateContainer("nEvents", TH1I::Class(), AliAnalysisManager::kOutputContainer, "HFEtask.root")); task->ConnectOutput(1, pidEffManager->CreateContainer("Results", TList::Class(), AliAnalysisManager::kOutputContainer, "HFEtask.root")); task->ConnectOutput(2, pidEffManager->CreateContainer("QA", TList::Class(), AliAnalysisManager::kOutputContainer, "HFEtask.root")); -- 2.43.0