From 7e57ef8fafe1e6fa3037628e946ea0de6ba33552 Mon Sep 17 00:00:00 2001 From: Raphaelle Bailhache Date: Thu, 14 Jul 2022 16:45:03 +0200 Subject: [PATCH 1/4] Dielectron efficiency task --- PWGEM/CMakeLists.txt | 2 +- PWGEM/Dilepton/CMakeLists.txt | 2 +- PWGEM/Dilepton/Tasks/CMakeLists.txt | 16 + PWGEM/Dilepton/Tasks/CMakeLists.txt~ | 45 + PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx | 1154 +++++++++++++++++++++ PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ | 1189 ++++++++++++++++++++++ 6 files changed, 2406 insertions(+), 2 deletions(-) create mode 100644 PWGEM/Dilepton/Tasks/CMakeLists.txt create mode 100644 PWGEM/Dilepton/Tasks/CMakeLists.txt~ create mode 100644 PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx create mode 100644 PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ diff --git a/PWGEM/CMakeLists.txt b/PWGEM/CMakeLists.txt index c6eb3a2030a..2b881e0ca4a 100644 --- a/PWGEM/CMakeLists.txt +++ b/PWGEM/CMakeLists.txt @@ -9,7 +9,7 @@ # granted to it by virtue of its status as an Intergovernmental Organization # or submit itself to any jurisdiction. -# add_subdirectory(Dilepton) +add_subdirectory(Dilepton) add_subdirectory(PhotonMeson) add_subdirectory(Tasks) diff --git a/PWGEM/Dilepton/CMakeLists.txt b/PWGEM/Dilepton/CMakeLists.txt index 626558b5be9..7565dace413 100644 --- a/PWGEM/Dilepton/CMakeLists.txt +++ b/PWGEM/Dilepton/CMakeLists.txt @@ -11,6 +11,6 @@ # add_subdirectory(Core) # add_subdirectory(DataModel) -# add_subdirectory(Tasks) +add_subdirectory(Tasks) # add_subdirectory(TableProducer) diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt new file mode 100644 index 00000000000..5c2c5e400fe --- /dev/null +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -0,0 +1,16 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + + +o2physics_add_dpl_workflow(efficiency-ee + SOURCES emEfficiencyEE.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt~ b/PWGEM/Dilepton/Tasks/CMakeLists.txt~ new file mode 100644 index 00000000000..56e990ba20b --- /dev/null +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt~ @@ -0,0 +1,45 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2physics_add_dpl_workflow(table-reader + SOURCES tableReader.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(efficiency + SOURCES dqEfficiency.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dilepton-ee + SOURCES dileptonEE.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dilepton-mumu + SOURCES dileptonMuMu.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(filter-pp + SOURCES filterPP.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(v0-selector + SOURCES v0selector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore O2::DetectorsVertexing + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(flow + SOURCES dqFlow.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx new file mode 100644 index 00000000000..5064efa6869 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx @@ -0,0 +1,1154 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// +// Analysis task for calculating single electron and dielectron efficiency +// +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "PWGDQ/Core/MCSignal.h" +#include "PWGDQ/Core/MCSignalLibrary.h" +#include +#include +#include +#include +#include +#include +#include + +using std::cout; +using std::endl; +using std::string; + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod; + +// Some definitions +namespace o2::aod +{ + +namespace emanalysisflags +{ + DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); + DECLARE_SOA_COLUMN(IsBarrelSelected, isBarrelSelected, int); +} // namespace reducedevent + + DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", emanalysisflags::IsEventSelected); + DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "BARRELTRACKCUTS", emanalysisflags::IsBarrelSelected); +} // namespace o2::aod + +// +using MyEvents = soa::Join; +using MyEventsSelected = soa::Join; +using MyBarrelTracks = soa::Join; +//using MyBarrelTracksWithCov = soa::Join; +using MyBarrelTracksSelected = soa::Join; +//using MyBarrelTracksSelectedWithCov = soa::Join; + + +// +constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended; +constexpr static uint32_t gkMCEventFillMap = VarManager::ObjTypes::ReducedEventMC; +constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; +//constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; +constexpr static uint32_t gkParticleMCFillMap = VarManager::ObjTypes::ParticleMC; + + + +void DefineHistograms(HistogramManager* histMan, TString histClasses); +void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps); + +struct AnalysisEventSelection { + + Produces eventSel; + OutputObj fOutputList{"output"}; + Configurable fConfigEventCuts{"cfgEventCuts", "eventStandard", "Event selection"}; + Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; + + HistogramManager* fHistMan; + AnalysisCompositeCut* fEventCut; + + void init(o2::framework::InitContext&) + { + fEventCut = new AnalysisCompositeCut(true); + TString eventCutStr = fConfigEventCuts.value; + fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + + VarManager::SetDefaultVarNames(); + if (fConfigQA) { + fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + fHistMan->SetUseDefaultVariableNames(kTRUE); + fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + DefineHistograms(fHistMan, "Event_BeforeCuts;Event_AfterCuts;"); // define all histograms + VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + fOutputList.setObject(fHistMan->GetMainHistogramList()); + } + } + + template + void runSelection(TEvent const& event, TEventsMC const& mcEvents) + { + // Reset the values array + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + + VarManager::FillEvent(event); + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + VarManager::FillEvent(event.reducedMCevent()); + } + if (fConfigQA) { + fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); // automatically fill all the histograms in the class Event + } + if (fEventCut->IsSelected(VarManager::fgValues)) { + if (fConfigQA) { + fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + } + eventSel(1); + } else { + eventSel(0); + } + } + + void processSkimmed(MyEvents::iterator const& event, aod::ReducedMCEvents const& mcEvents) + { + runSelection(event, mcEvents); + } + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisEventSelection, processSkimmed, "Run event selection on DQ skimmed events", false); + PROCESS_SWITCH(AnalysisEventSelection, processDummy, "Dummy process function", false); +}; +struct AnalysisEventQa { + + Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; + OutputObj fOutputList{"output"}; + THashList *fMainList; + TH1D *fNbRecCollisionPerMCCollision; + TProfile *fNbmctrack; + TH1D *fNbMcEvent; + TH1D *fNbRecEvent; + + void init(o2::framework::InitContext&) + { + + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + VarManager::SetDefaultVarNames(); + + fMainList = new THashList; + fMainList->SetOwner(kTRUE); + fMainList->SetName("OAnalysisEventQA"); + fNbRecCollisionPerMCCollision = new TH1D("MCEvRecEv","",20,0.,20.); + fNbmctrack = new TProfile("mctrack","",20,0.,20.); + fNbMcEvent = new TH1D("MCEvent","",1,0.,1.); + fNbRecEvent = new TH1D("RecEvent","",1,0.,1.); + fMainList->Add(fNbRecCollisionPerMCCollision); + fMainList->Add(fNbmctrack); + fMainList->Add(fNbMcEvent); + fMainList->Add(fNbRecEvent); + fOutputList.setObject(fMainList); + + } + + Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + + template + void runSelection(TEvents const& events, TEventsMC const& eventsMC, TTracksMC const& tracksMC) + { + + uint8_t eventFilter = 0; + std::map fMCEventNbmctrack; + std::map fMCEventNbReco; + std::map fMCEventLabels; + + + int fMCCounters = 0; //! - mc event counter + int fEvCounters = 0; // ! - rec event counter + + + // First loop + + for (auto& event : events) { + + + // Reset the values array + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event); + eventFilter = uint32_t(event.isEventSelected()); + if(!eventFilter) continue; + fEvCounters++; + fNbRecEvent->Fill(0.5); + + + Int_t midrap = 0; + Int_t globalindexmc = -1; + + // skimmed data + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + midrap = dNdetach(groupedMCTracks); + + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + + } + + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fMCCounters; + fMCEventNbReco[globalindexmc] = 1; + fMCEventNbmctrack[globalindexmc] = midrap; + fNbMcEvent->Fill(0.5); + fMCCounters++; + + } + else { + fMCEventNbReco[globalindexmc] = fMCEventNbReco.find(globalindexmc)->second + 1; + } + + + } // end loop over events + + for (const auto& [mcEv, NbRecEv] : fMCEventNbReco) { + + fNbRecCollisionPerMCCollision->Fill(fMCEventNbReco.find(mcEv)->second); + + } + + for (const auto& [mcEv, NbRecEv] : fMCEventNbmctrack) { + + fNbmctrack->Fill(fMCEventNbReco.find(mcEv)->second,fMCEventNbmctrack.find(mcEv)->second); + + } + + } + + template + Int_t dNdetach(TTracksMC const& groupedMCTracks){ + + Int_t midrap = 0; + for(auto mctrack : groupedMCTracks){ + if(TMath::Abs(mctrack.eta())<0.5 && mctrack.isPhysicalPrimary() && (TMath::Abs(mctrack.pdgCode())==211 || mctrack.pdgCode()==111)) { + midrap++; + } + } + return midrap; + + } + + + void processSkimmed(soa::Filtered const&events, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + runSelection(events, eventsMC, tracksMC); + } + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisEventQa, processSkimmed, "Run event selection on DQ skimmed events", false); + PROCESS_SWITCH(AnalysisEventQa, processDummy, "Dummy process function", false); + +}; + +struct AnalysisTrackSelection { + + Produces trackSel; + Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; + + // configurables + //Configurable fConfigEventCuts{"cfgEventCuts", "eventStandardNoINT7", "Event selection"}; + Configurable fConfigCuts{"cfgTrackCuts", "jpsiPID1", "Comma separated list of barrel track cuts"}; + Configurable fConfigMCSignals{"cfgTrackMCSignals", "", "Comma separated list of MC signals"}; + + // 3D histos + Configurable fConfigUsePtVec{"cfgUsePtVec", true, "If true, non-linear pt bins predefined"}; + Configurable fConfigMinPt{"cfgMinPt", 0., "min Pt in 3D histos"}; + Configurable fConfigMaxPt{"cfgMaxPt", 10., "max Pt in 3D histos"}; + Configurable fConfigStepPt{"cfgStepPt", 1, "Nb of steps in pt in 3D histos"}; + Configurable fConfigMinEta{"cfgMinEta", -0.8, "min Eta in 3D histos"}; + Configurable fConfigMaxEta{"cfgMaxEta", 0.8, "max Eta in 3D histos"}; + Configurable fConfigStepEta{"cfgStepEta", 16, "Nb of steps in Eta in 3D histos"}; + Configurable fConfigMinPhi{"cfgMinPhi", 0., "min Phi in 3D histos"}; + Configurable fConfigMaxPhi{"cfgMaxPhi", 6.3, "max Phi in 3D histos"}; + Configurable fConfigStepPhi{"cfgStepPhi", 63, "Nb of steps in Phi in 3D histos"}; + + // Resolution histos + Configurable fConfigResolutionOn{"cfgResolution", false, "If true, fill resolution histograms"}; + Configurable fConfigUsePtVecRes{"cfgUsePtVecRes", true, "If true, non-linear pt bins predefined in res histos"}; + Configurable fConfigMinPtRes{"cfgMinPtRes", 0., "min Pt in res histos"}; + Configurable fConfigMaxPtRes{"cfgMaxPtRes", 20., "max Pt in res histos"}; + Configurable fConfigStepPtRes{"cfgStepPtRes", 1, "Nb of steps in pt in res histos"}; + Configurable fConfigStepDeltaPt{"cfgStepDeltaPt", 1, "Nb of steps in delta pt in res histos"}; + Configurable fConfigMinDeltaEta{"cfgMinDeltaEta", -0.5, "min delta Eta in res histos"}; + Configurable fConfigMaxDeltaEta{"cfgMaxDeltaEta", 0.5, "max delta Eta in res histos"}; + Configurable fConfigStepDeltaEta{"cfgStepDeltaEta", 500, "Nb of steps in detla Eta in res histos"}; + Configurable fConfigMinDeltaPhi{"cfgMinDeltaPhi", -0.5, "min delta Phi in res histos"}; + Configurable fConfigMaxDeltaPhi{"cfgMaxDeltaPhi", 0.5, "max delta Phi in res histos"}; + Configurable fConfigStepDeltaPhi{"cfgStepDeltaPhi", 500, "Nb of steps in delta Phi in res histos"}; + + + Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; + + // output lists + OutputObj fOutputList{"output"}; + THashList *fMainList; // Main list + THashList *fSingleElectronList; // 3D histos for MC and reconstructed signals + THashList *fResolutionList; // Resolution histograms + THashList *fQASingleElectronList; // QA in case on with histo manager outputs + + // Cuts and signals + //AnalysisCompositeCut* fEventCut; // Taken from event selection part + std::vector fTrackCuts; // list of track cuts + AnalysisCompositeCut* fTrackCutsRes; // track cut for resolution map + std::vector fMCSignals; // list of signals to be checked + MCSignal* fMCSignalRes; // signal for res + + //3D histos + std::vector fHistGenPosPart; + std::vector fHistGenNegPart; + std::vector fHistGenSmearedPosPart; + std::vector fHistGenSmearedNegPart; + std::vector fHistRecPosPart; + std::vector fHistRecNegPart; + // Binning + std::vector fPtBins; + std::vector fEtaBins; + std::vector fPhiBins; + + // Res histos + std::vector fHistRes; + // Binning + std::vector fPtResBins; + std::vector fDeltaEtaBins; + std::vector fDeltaPhiBins; + + // QA + HistogramManager* fHistManQA; // histo manager + std::vector fHistNamesRecoQA; // list of histo names for all reconstructed tracks in histo manager + std::vector> fHistNamesMCMatchedQA; // list of histo names for reconstructed signals in histo manager + std::vector> fHistNamesMCQA; // list of histo names for generated signals in histo manager + + void init(o2::framework::InitContext&) + { + + // Create list output + fMainList = new THashList; + fMainList->SetOwner(kTRUE); + fMainList->SetName("trackselection"); + + // Create list output for 3D eta,phi,pt + fSingleElectronList = new THashList; + fSingleElectronList->SetOwner(kTRUE); + fSingleElectronList->SetName("SingleElectron"); + + // Binning 3D histos + if(fConfigUsePtVec){ + const Int_t Npt = 68; + Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.}; + std::vector v_pte(pte,std::end(pte)); + fPtBins = v_pte; + } + else { + SetBinsLinear(fPtBins, fConfigMinPt, fConfigMaxPt, fConfigStepPt); + } + SetBinsLinear(fEtaBins, fConfigMinEta, fConfigMaxEta, fConfigStepEta); + SetBinsLinear(fPhiBins, fConfigMinPhi, fConfigMaxPhi, fConfigStepPhi); + const int fNptBins = fPtBins.size()-1; + const int fNetaBins = fEtaBins.size()-1; + const int fNphiBins = fPhiBins.size()-1; + + // Event cut: taken from event selection part + //fEventCut = new AnalysisCompositeCut(true); + //TString eventCutStr = fConfigEventCuts.value; + //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + + // List of track cuts + TString cutNamesStr = fConfigCuts.value; + if (!cutNamesStr.IsNull()) { + std::unique_ptr objArray(cutNamesStr.Tokenize(",")); + for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + fTrackCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); + } + } + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + VarManager::SetDefaultVarNames(); + + // List of MC signals + TString configSigNamesStr = fConfigMCSignals.value; + std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); + for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() != 1) { // NOTE: only 1 prong signals + continue; + } + // List of signal to be checked + fMCSignals.push_back(*sig); + } + } + + //Configure 3D histograms + // Create List with generated particles + TList* Generated = new TList(); + Generated->SetOwner(); + Generated->SetName("Generated"); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_pos->Sumw2(); + fHistGenPosPart.push_back(th3_tmp_pos); + Generated->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_neg->Sumw2(); + fHistGenNegPart.push_back(th3_tmp_neg); + Generated->Add(th3_tmp_neg); + } + // Create List with generated+smeared particles + TList* GeneratedSmeared = new TList(); + GeneratedSmeared->SetName("GeneratedSmeared"); + GeneratedSmeared->SetOwner(); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_pos->Sumw2(); + fHistGenSmearedPosPart.push_back(th3_tmp_pos); + GeneratedSmeared->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_neg->Sumw2(); + fHistGenSmearedNegPart.push_back(th3_tmp_neg); + GeneratedSmeared->Add(th3_tmp_neg); + } + + fSingleElectronList->Add(Generated); + fSingleElectronList->Add(GeneratedSmeared); + + // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge + for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ + TList* list = new TList(); + list->SetName(fTrackCuts.at(list_i).GetName()); + list->SetOwner(); + + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH3D* th3_tmp_pos = new TH3D(Form("Nrec_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_pos->Sumw2(); + th3_tmp_pos->SetDirectory(0x0); + fHistRecPosPart.push_back(th3_tmp_pos); + list->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Nrec_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_neg->Sumw2(); + th3_tmp_neg->SetDirectory(0x0); + fHistRecNegPart.push_back(th3_tmp_neg); + list->Add(th3_tmp_neg); + + } + fSingleElectronList->Add(list); + } + fMainList->Add(fSingleElectronList); + + // Resolution histogramms + if(fConfigResolutionOn) { + + // Binning 3D histos + if(fConfigUsePtVecRes){ + const Int_t Npt = 73; + Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.,12.0,14.,16.,18.,20.}; + std::vector v_pte(pte,std::end(pte)); + fPtResBins = v_pte; + } + else { + SetBinsLinear(fPtResBins, fConfigMinPtRes, fConfigMaxPtRes, fConfigStepPtRes); + } + SetBinsLinear(fDeltaEtaBins, fConfigMinDeltaEta, fConfigMaxDeltaEta, fConfigStepDeltaEta); + SetBinsLinear(fDeltaPhiBins, fConfigMinDeltaPhi, fConfigMaxDeltaPhi, fConfigStepDeltaPhi); + const int fNptresBins = fPtResBins.size()-1; + const int fNDeltaetaBins = fDeltaEtaBins.size()-1; + const int fNDeltaphiBins = fDeltaPhiBins.size()-1; + + + fResolutionList = new THashList; + fResolutionList->SetOwner(kTRUE); + fResolutionList->SetName("Resolution"); + + printf("Histos\n"); + TH2D *thPtGen_DeltaPtOverPtGen = new TH2D("PtGen_DeltaPtOverPtGen", "", fNptresBins, fPtResBins.data(), fConfigStepDeltaPt, -1., +1.); + TH2D *thPtGen_DeltaEta = new TH2D("PtGen_DeltaEta", "", fNptresBins, fPtResBins.data(), fNDeltaetaBins, fDeltaEtaBins.data()); + TH2D *thPtGen_DeltaPhi_Ele = new TH2D("PtGen_DeltaPhi_Ele", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); + TH2D *thPtGen_DeltaPhi_Pos = new TH2D("PtGen_DeltaPhi_Pos", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); + + thPtGen_DeltaPtOverPtGen ->Sumw2(); + thPtGen_DeltaEta ->Sumw2(); + thPtGen_DeltaPhi_Ele ->Sumw2(); + thPtGen_DeltaPhi_Pos ->Sumw2(); + + thPtGen_DeltaPtOverPtGen ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPtOverPtGen ->GetYaxis()->SetTitle("(p^{gen}_{T} - p^{rec}_{T}) / p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaEta ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaEta ->GetYaxis()->SetTitle("#eta^{gen} - #eta^{rec}"); + thPtGen_DeltaPhi_Ele ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPhi_Ele ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); + thPtGen_DeltaPhi_Pos ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPhi_Pos ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); + + printf("Add\n"); + fHistRes.push_back(thPtGen_DeltaPtOverPtGen); + fHistRes.push_back(thPtGen_DeltaEta); + fHistRes.push_back(thPtGen_DeltaPhi_Ele); + fHistRes.push_back(thPtGen_DeltaPhi_Pos); + + fResolutionList->Add(thPtGen_DeltaPtOverPtGen); + fResolutionList->Add(thPtGen_DeltaEta); + fResolutionList->Add(thPtGen_DeltaPhi_Ele); + fResolutionList->Add(thPtGen_DeltaPhi_Pos); + + fMainList->Add(fResolutionList); + + } + + + // Configure QA histogram classes + if (fConfigQA){ + // Create list output for QA + fQASingleElectronList = new THashList; + fQASingleElectronList->SetOwner(kTRUE); + fQASingleElectronList->SetName("SingleElectronQA"); + TString histClassesQA = ""; + for (auto& cut : fTrackCuts) { + + // All reconstructed leptons + TString nameStr = Form("TrackBarrel_%s", cut.GetName()); + fHistNamesRecoQA.push_back(nameStr); + histClassesQA += Form("%s;", nameStr.Data()); + + // All reconstructed leptons matched to a 1 prong signal or MC 1 prong signal directly + std::vector mcnamesreco; + for (int isig = 0; isig < fMCSignals.size(); ++isig) { + TString nameStr2 = Form("TrackBarrel_%s_%s", cut.GetName(), fMCSignals.at(isig).GetName()); + mcnamesreco.push_back(nameStr2); + histClassesQA += Form("%s;", nameStr2.Data()); + } + fHistNamesMCMatchedQA.push_back(mcnamesreco); + } + + // Add histogram classes for each MC signal at generated level + std::vector mcnamesgen; + for (int isig = 0; isig < fMCSignals.size(); ++isig) { + TString nameStr2 = Form("MCTruthGen_%s", fMCSignals.at(isig).GetName()); + mcnamesgen.push_back(nameStr2); + histClassesQA += Form("%s;", nameStr2.Data()); + + } + fHistNamesMCQA.push_back(mcnamesgen); + + + fHistManQA = new HistogramManager("SingleElectronQA", "aa", VarManager::kNVars); + fHistManQA->SetUseDefaultVariableNames(kTRUE); + fHistManQA->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + DefineHistograms(fHistManQA, histClassesQA.Data()); // define all histograms + VarManager::SetUseVars(fHistManQA->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + fQASingleElectronList = fHistManQA->GetMainHistogramList(); + fMainList->Add(fQASingleElectronList); + + } + + fOutputList.setObject(fMainList); + +} + + Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; + + template + void runSelection(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) + { + + uint8_t eventFilter = 0; + std::map fMCEventLabels; + int fCounters = 0; //! [0] - particle counter, [1] - event counter + + + for (auto& event : events) { + + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::ResetValues(0, VarManager::kNMCParticleVariables); + // fill event information which might be needed in histograms that combine track and event properties + VarManager::FillEvent(event); + //if(!fEventCut->IsSelected(VarManager::fgValues)) continue; + eventFilter = uint32_t(event.isEventSelected()); + if(!eventFilter) continue; + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + VarManager::FillEvent(event.reducedMCevent()); + } + + // Look if we did not already saw the collision and fill the denominator of the single electron efficiency + Int_t globalindexmc = -1; + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + } + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fCounters; + fCounters++; + // skimmed data + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGenTrack(groupedMCTracks); + } + } + + + // Loop over reconstructed tracks belonging to the event and fill the numerator of the efficiency as well as the resolution map + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); + runRecTrack(groupedTracks,tracksMC); + } + + }// end loop over events + + } + + template + void runMCGenTrack(TTracksMC const& groupedMCTracks) + { + + + // loop over mc stack and fill histograms for pure MC truth signals + // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event + //auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); + for (auto& mctrack : groupedMCTracks) { + VarManager::FillTrack(mctrack); + // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. + // NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member. + // TODO: Use the mcReducedFlags to select signals + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + if ((*sig).CheckSignal(true, groupedMCTracks, mctrack)) { + if(mctrack.pdgCode() > 0) { + dynamic_cast(fHistGenNegPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); + } + else { + dynamic_cast(fHistGenPosPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); + } + if (fConfigQA) fHistManQA->FillHistClass(Form("MCTruthGen_%s", (*sig).GetName()), VarManager::fgValues); + } + } + } + } + + template + void runRecTrack(TTracks const& groupedTracks, TTracksMC const& tracksMC) + { + + uint32_t filterMap = 0; + trackSel.reserve(groupedTracks.size()); + for (auto& track : groupedTracks) { + filterMap = 0; + + VarManager::FillTrack(track); // compute track quantities + + // compute MC matched quantities + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + VarManager::FillTrack(track.reducedMCTrack()); + } + + + // compute track selection and publish the bit map + int i = 0; + for (auto cut = fTrackCuts.begin(); cut != fTrackCuts.end(); cut++, i++) { + if ((*cut).IsSelected(VarManager::fgValues)) { + filterMap |= (uint32_t(1) << i); + if (fConfigQA) { + fHistManQA->FillHistClass(fHistNamesRecoQA[i].Data(), VarManager::fgValues); + } + } + } + trackSel(static_cast(filterMap)); + if (!filterMap) { + continue; + } + + // compute MC matching decisions + uint32_t mcDecision = 0; + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + if ((*sig).CheckSignal(true, tracksMC, track.reducedMCTrack())) { + mcDecision |= (uint32_t(1) << isig); + } + } + + } + + + // fill histograms + for (unsigned int i = 0; i < fMCSignals.size(); i++) { + if (!(mcDecision & (uint32_t(1) << i))) { + continue; + } + for (unsigned int j = 0; j < fTrackCuts.size(); j++) { + if (filterMap & (uint8_t(1) << j)) { + if(track.sign() < 0) { + dynamic_cast(fHistRecNegPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); + } + else { + dynamic_cast(fHistRecPosPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); + } + if(fConfigResolutionOn && (i==0) && (j==0)){ + + Double_t mcpt = -10000.; + Double_t mceta = -10000.; + Double_t mcphi = -1000.; + Int_t mcpdg = -10000.; + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto mctrack = track.reducedMCTrack(); + mcpt = mctrack.pt(); + mceta = mctrack.eta(); + mcphi = mctrack.phi(); + mcpdg = mctrack.pdgCode(); + } + + Double_t deltaptoverpt = -1000.; + if(mcpt>0.) deltaptoverpt = (mcpt-track.pt())/mcpt; + Double_t deltaeta = mceta-track.eta(); + Double_t deltaphi = mcphi-track.phi(); + dynamic_cast(fHistRes.at(0))->Fill(mcpt,deltaptoverpt); + dynamic_cast(fHistRes.at(1))->Fill(mcpt,deltaeta); + if(mcpdg<0){ + dynamic_cast(fHistRes.at(2))->Fill(mcpt,deltaphi); + }else { + dynamic_cast(fHistRes.at(3))->Fill(mcpt,deltaphi); + } + } + if(fConfigQA) fHistManQA->FillHistClass(fHistNamesMCMatchedQA[j][i].Data(), VarManager::fgValues); + } + } // end loop over cuts + } // end loop over MC signals + + + } // end loop over reconstructed track belonging to the events + + } + + + + void processSkimmed(soa::Filtered const& events, MyBarrelTracks const& tracks, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + runSelection(events, tracks, eventsMC, tracksMC); + + + } + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisTrackSelection, processSkimmed, "Run barrel track selection on DQ skimmed tracks", false); + PROCESS_SWITCH(AnalysisTrackSelection, processDummy, "Dummy process function", false); +}; + + +struct AnalysisSameEventPairing { + + + //Produces dileptonList; + + // Filter based on previous components in the task + Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; + Filter filterBarrelTrackSelected = aod::emanalysisflags::isBarrelSelected > 0; + + Configurable fConfigTrackCuts{"cfgTrackCuts", "", "Comma separated list of barrel track cuts"}; + Configurable fConfigMCSignals{"cfgBarrelMCSignals", "", "Comma separated list of MC signals"}; + Configurable fConfigMinPt{"cfgMinPt", 0., "Fiducial min Pt for MC signal"}; + Configurable fConfigMaxPt{"cfgMaxPt", 10., "Fiducial max Pt for MC signal"}; + Configurable fConfigMinEta{"cfgMinEta", -0.8, "Fiducial min eta for MC signal"}; + Configurable fConfigMaxEta{"cfgMaxEta", 0.8, "Fiducial max eta for MC signal"}; + Configurable fConfigFlatTables{"cfgFlatTables", false, "Produce a single flat tables with all relevant information of the pairs and single tracks"}; + + // TODO: here we specify signals, however signal decisions are precomputed and stored in mcReducedFlags + // TODO: The tasks based on skimmed MC could/should rely ideally just on these flags + // TODO: special AnalysisCuts to be prepared in this direction + // TODO: cuts on the MC truth information to be added if needed + + + // 2D histos: mee and ptee + Configurable fConfigMinPtee{"cfgMinPtee", 0., "min Ptee in 2D histos"}; + Configurable fConfigMaxPtee{"cfgMaxPtee", 10., "max Ptee in 2D histos"}; + Configurable fConfigStepPtee{"cfgStepPtee", 100, "Nb of steps in ptee in 2D histos"}; + Configurable fConfigMinMee{"cfgMinMee", 0., "min Mee in 2D histos"}; + Configurable fConfigMaxMee{"cfgMaxMee", 3.5, "max Mee in 2D histos"}; + Configurable fConfigStepMee{"cfgStepMee", 600, "Nb of steps in Mee in 2D histos"}; + + // output lists + OutputObj fOutputList{"output"}; + THashList *fMainList; // Main list + THashList *fPairList; // 2D histos for MC and reconstructed signals + THashList *fQAPairList; // QA in case on with histo manager outputs + + // Cuts and signals + //AnalysisCompositeCut* fEventCut; // Taken from event selection part + std::vector fTrackCuts; // list of track cuts + std::vector fMCSignals; // list of signals with one prong to be checked: ULS 2D histos + + //2D histo vectors + std::vector fHistGenPair; + std::vector fHistGenSmearedPair; + std::vector fHistRecPair; + // Binning + std::vector fPteeBins; + std::vector fMeeBins; + + // QA: to be defined + + + void init(o2::framework::InitContext& context) + { + + printf("Start init\n"); + + // Create list output + fMainList = new THashList; + fMainList->SetOwner(kTRUE); + fMainList->SetName("pairselection"); + + // Create list output for 3D eta,phi,pt + fPairList = new THashList; + fPairList->SetOwner(kTRUE); + fPairList->SetName("Dielectron"); + + // Binning 2D histos + SetBinsLinear(fPteeBins, fConfigMinPtee, fConfigMaxPtee, fConfigStepPtee); + SetBinsLinear(fMeeBins, fConfigMinMee, fConfigMaxMee, fConfigStepMee); + const int fNpteeBins = fPteeBins.size()-1; + const int fNmeeBins = fMeeBins.size()-1; + + // Event cut: taken from event selection part + //fEventCut = new AnalysisCompositeCut(true); + //TString eventCutStr = fConfigEventCuts.value; + //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + + // List of track cuts + TString cutNamesStr = fConfigTrackCuts.value; + if (!cutNamesStr.IsNull()) { + std::unique_ptr objArray(cutNamesStr.Tokenize(",")); + for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + fTrackCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); + } + } + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + VarManager::SetDefaultVarNames(); + + // List of MC signals + TString configSigNamesStr = fConfigMCSignals.value; + std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); + for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() == 2) { // only 2 prong signals + fMCSignals.push_back(*sig); + } + // List of signal to be checked + } + } + + //Configure 2D histograms + // Create List with generated particles + TList* Generated = new TList(); + Generated->SetOwner(); + Generated->SetName("Generated"); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH2D* th2_tmp = new TH2D(Form("Ngen_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + th2_tmp->Sumw2(); + fHistGenPair.push_back(th2_tmp); + Generated->Add(th2_tmp); + } + + // Create List with generated+smeared particles + TList* GeneratedSmeared = new TList(); + GeneratedSmeared->SetName("GeneratedSmeared"); + GeneratedSmeared->SetOwner(); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH2D* th2_tmp = new TH2D(Form("NgenSmeared_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + th2_tmp->Sumw2(); + fHistGenSmearedPair.push_back(th2_tmp); + GeneratedSmeared->Add(th2_tmp); + } + + + fPairList->Add(Generated); + fPairList->Add(GeneratedSmeared); + + // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge + for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ + TList* list = new TList(); + list->SetName(fTrackCuts.at(list_i).GetName()); + list->SetOwner(); + + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH2D* th2_tmp = new TH2D(Form("Nrec_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + th2_tmp->Sumw2(); + th2_tmp->SetDirectory(0x0); + fHistRecPair.push_back(th2_tmp); + list->Add(th2_tmp); + } + fPairList->Add(list); + } + fMainList->Add(fPairList); + + fOutputList.setObject(fMainList); + + } + + Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; + + template < uint32_t TEventFillMap, uint32_t TEventMCFillMap, uint32_t TTrackFillMap, typename TEvents, typename TTracks, typename TEventsMC, typename TTracksMC> + void runPairing(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) + { + + + std::map fMCEventLabels; + int fCounters = 0; //! [0] - particle counter, [1] - event counter + + + for (auto& event : events) { + + if (!event.isEventSelected()) { + return; + } + + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::ResetValues(0, VarManager::kNMCParticleVariables); + // fill event information which might be needed in histograms that combine track and event properties + VarManager::FillEvent(event); + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + VarManager::FillEvent(event.reducedMCevent()); + } + + // Look if we did not already saw the collision and fill the denominator of the single electron efficiency + Int_t globalindexmc = -1; + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + } + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fCounters; + fCounters++; + // skimmed data + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGenPair(groupedMCTracks); + } + } + + + //auto groupedTrackspos = posTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); + //groupedTrackspos.bindInternalIndicesTo(&tracks); + //auto groupedTracksneg = negTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); + //groupedTracksneg.bindInternalIndicesTo(&tracks); + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); + groupedTracks.bindInternalIndicesTo(&tracks); + runRecPair(groupedTracks,tracksMC); + } + + } // end loop over reconstructed event + } // end loop pairing function + + + template + void runMCGenPair(TTracksMC const& groupedMCTracks) + { + // + Double_t masse = 0.00051099895; // 0.5 MeV/c2 -> 0.0005 GeV/c2 + + for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { + + if((abs(t1.pdgCode()) != 11) || (abs(t2.pdgCode()) != 11)) continue; + if(t1.pdgCode()*t2.pdgCode() > 0) continue; // ULS only + + TLorentzVector Lvec1; + TLorentzVector Lvec2; + Lvec1.SetPtEtaPhiM(t1.pt(),t1.eta(),t1.phi(),masse); + Lvec2.SetPtEtaPhiM(t2.pt(),t2.eta(),t2.phi(),masse); + TLorentzVector LvecM = Lvec1 + Lvec2; + double mass = LvecM.M(); + double pairpt = LvecM.Pt(); + //double opangle = Lvec1.Angle(Lvec2.Vect()); + + //printf("Check before\n"); + // Fiducial cut + Bool_t genfidcut = kTRUE; + if((t1.eta()> fConfigMaxEta) || (t2.eta() > fConfigMaxEta) || (t1.eta() < fConfigMinEta) || (t2.eta() < fConfigMinEta) || (t1.pt()> fConfigMaxPt) || (t2.pt() > fConfigMaxPt) || (t1.pt() < fConfigMinPt) || (t2.pt() < fConfigMinPt)) genfidcut = kFALSE; + + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + if ((*sig).CheckSignal(true, groupedMCTracks, t1, t2)) { + + // not smeared + if(genfidcut) dynamic_cast(fHistGenPair.at(isig))->Fill(mass,pairpt); + // need to implement smeared + + } + } + } //end of true pairing loop + } // end runMCGen + + + + template < uint32_t TTrackFillMap, typename TTracks, typename TTracksMC> + void runRecPair(TTracks const& tracks, TTracksMC const& tracksMC) + { + + // Loop over two track combinations + uint8_t twoTrackFilter = 0; + //uint32_t dileptonFilterMap = 0; + // uint32_t dileptonMcDecision = 0; + //dileptonList.reserve(1); + + + + for (auto& [t1, t2] : combinations(tracks, tracks)) { + + twoTrackFilter = uint32_t(t1.isBarrelSelected()) & uint32_t(t2.isBarrelSelected()); + + if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + continue; + } + VarManager::FillPair(t1, t2); + + // run MC matching for this pair + uint32_t mcDecision = 0; + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { // for skimmed DQ model + if ((*sig).CheckSignal(true, tracksMC, t1.reducedMCTrack(), t2.reducedMCTrack())) { + mcDecision |= (uint32_t(1) << isig); + } + } + + }// end of loop MC signals + + //dileptonFilterMap = twoTrackFilter; + //dileptonMcDecision = mcDecision; + //dileptonList(event, VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), dileptonFilterMap, dileptonMcDecision); + + + for (unsigned int i = 0; i < fMCSignals.size(); i++) { + if (!(mcDecision & (uint32_t(1) << i))) { + continue; + } + for (unsigned int j = 0; j < fTrackCuts.size(); j++) { + if (twoTrackFilter & (uint8_t(1) << j)) { + dynamic_cast(fHistRecPair.at(j * fMCSignals.size() + i))->Fill(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt]); + } + } + } + } + + } + + void processToEESkimmed(soa::Filtered const& events, + soa::Filtered const& tracks, + ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + + runPairing< gkEventFillMap, gkMCEventFillMap, gkTrackFillMap>(events, tracks, eventsMC, tracksMC); + + } + + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisSameEventPairing, processToEESkimmed, "Run barrel barrel pairing on DQ skimmed tracks", false); + PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy process function", false); +}; + + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} + +void DefineHistograms(HistogramManager* histMan, TString histClasses) +{ + // + // Define here the histograms for all the classes required in analysis. + // The histogram classes are provided in the histClasses string, separated by semicolon ";" + // The histogram classes and their components histograms are defined below depending on the name of the histogram class + // + std::unique_ptr objArray(histClasses.Tokenize(";")); + for (Int_t iclass = 0; iclass < objArray->GetEntries(); ++iclass) { + TString classStr = objArray->At(iclass)->GetName(); + histMan->AddHistClass(classStr.Data()); + + // NOTE: The level of detail for histogramming can be controlled via configurables + if (classStr.Contains("Event")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", "trigger,cent,mc"); + } + + if (classStr.Contains("Track")) { + if (classStr.Contains("Barrel")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", "its,tpcpid,dca,tofpid,mc"); + } + } + + if (classStr.Contains("Pairs")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel", "vertexing-barrel"); + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_dimuon", "vertexing-forward"); + } + + if (classStr.Contains("MCTruthGenPair")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair"); + histMan->AddHistogram(objArray->At(iclass)->GetName(), "Pt", "MC generator p_{T} distribution", false, 200, 0.0, 20.0, VarManager::kMCPt); + histMan->AddHistogram(objArray->At(iclass)->GetName(), "Eta", "MC generator #eta distribution", false, 500, -5.0, 5.0, VarManager::kMCEta); + histMan->AddHistogram(objArray->At(iclass)->GetName(), "Phi", "MC generator #varphi distribution", false, 500, -6.3, 6.3, VarManager::kMCPhi); + } + if (classStr.Contains("MCTruthGen")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth"); + } + if (classStr.Contains("DileptonsSelected")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel"); + } + + + } +}// end loop over histogram classes + +void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps){ + fBins.clear(); + const double stepSize = (max - min) / steps; + for (unsigned int i = 0; i < steps+1; ++i){ + fBins.push_back(i * stepSize + min); + } +} diff --git a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ new file mode 100644 index 00000000000..50c71d4eeab --- /dev/null +++ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ @@ -0,0 +1,1189 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// +// Analysis task for calculating single electron and dielectron efficiency +// +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "PWGDQ/Core/MCSignal.h" +#include "PWGDQ/Core/MCSignalLibrary.h" +#include +#include +#include +#include +#include +#include +#include + +using std::cout; +using std::endl; +using std::string; + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod; + +// Some definitions +namespace o2::aod +{ + +namespace emanalysisflags +{ + DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); + DECLARE_SOA_COLUMN(IsBarrelSelected, isBarrelSelected, int); +} // namespace reducedevent + + DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", emanalysisflags::IsEventSelected); + DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "BARRELTRACKCUTS", emanalysisflags::IsBarrelSelected); +} // namespace o2::aod + +// +using MyEvents = soa::Join; +using MyEventsSelected = soa::Join; +using MyBarrelTracks = soa::Join; +//using MyBarrelTracksWithCov = soa::Join; +using MyBarrelTracksSelected = soa::Join; +//using MyBarrelTracksSelectedWithCov = soa::Join; + + +// +constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended; +constexpr static uint32_t gkMCEventFillMap = VarManager::ObjTypes::ReducedEventMC; +constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; +//constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; +constexpr static uint32_t gkParticleMCFillMap = VarManager::ObjTypes::ParticleMC; + + + +void DefineHistograms(HistogramManager* histMan, TString histClasses); +void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps); + +struct AnalysisEventSelection { + + Produces eventSel; + OutputObj fOutputList{"output"}; + Configurable fConfigEventCuts{"cfgEventCuts", "eventStandard", "Event selection"}; + Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; + + HistogramManager* fHistMan; + AnalysisCompositeCut* fEventCut; + + void init(o2::framework::InitContext&) + { + fEventCut = new AnalysisCompositeCut(true); + TString eventCutStr = fConfigEventCuts.value; + fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + + VarManager::SetDefaultVarNames(); + if (fConfigQA) { + fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + fHistMan->SetUseDefaultVariableNames(kTRUE); + fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + DefineHistograms(fHistMan, "Event_BeforeCuts;Event_AfterCuts;"); // define all histograms + VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + fOutputList.setObject(fHistMan->GetMainHistogramList()); + } + } + + template + void runSelection(TEvent const& event, TEventsMC const& mcEvents) + { + // Reset the values array + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + + VarManager::FillEvent(event); + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + VarManager::FillEvent(event.reducedMCevent()); + } + if (fConfigQA) { + fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); // automatically fill all the histograms in the class Event + } + if (fEventCut->IsSelected(VarManager::fgValues)) { + if (fConfigQA) { + fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + } + eventSel(1); + } else { + eventSel(0); + } + } + + void processSkimmed(MyEvents::iterator const& event, aod::ReducedMCEvents const& mcEvents) + { + runSelection(event, mcEvents); + } + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisEventSelection, processSkimmed, "Run event selection on DQ skimmed events", false); + PROCESS_SWITCH(AnalysisEventSelection, processDummy, "Dummy process function", false); +}; +struct AnalysisEventQa { + + Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; + OutputObj fOutputList{"output"}; + THashList *fMainList; + TH1D *fNbRecCollisionPerMCCollision; + TProfile *fNbmctrack; + TH1D *fNbMcEvent; + TH1D *fNbRecEvent; + + void init(o2::framework::InitContext&) + { + + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + VarManager::SetDefaultVarNames(); + + fMainList = new THashList; + fMainList->SetOwner(kTRUE); + fMainList->SetName("OAnalysisEventQA"); + fNbRecCollisionPerMCCollision = new TH1D("MCEvRecEv","",20,0.,20.); + fNbmctrack = new TProfile("mctrack","",20,0.,20.); + fNbMcEvent = new TH1D("MCEvent","",1,0.,1.); + fNbRecEvent = new TH1D("RecEvent","",1,0.,1.); + fMainList->Add(fNbRecCollisionPerMCCollision); + fMainList->Add(fNbmctrack); + fMainList->Add(fNbMcEvent); + fMainList->Add(fNbRecEvent); + fOutputList.setObject(fMainList); + + } + + Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + + template + void runSelection(TEvents const& events, TEventsMC const& eventsMC, TTracksMC const& tracksMC) + { + + uint8_t eventFilter = 0; + std::map fMCEventNbmctrack; + std::map fMCEventNbReco; + std::map fMCEventLabels; + + + int fMCCounters = 0; //! - mc event counter + int fEvCounters = 0; // ! - rec event counter + + + // First loop + + for (auto& event : events) { + + + // Reset the values array + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event); + eventFilter = uint32_t(event.isEventSelected()); + if(!eventFilter) continue; + fEvCounters++; + fNbRecEvent->Fill(0.5); + + + Int_t midrap = 0; + Int_t globalindexmc = -1; + + // skimmed data + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + midrap = dNdetach(groupedMCTracks); + + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + + } + + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fMCCounters; + fMCEventNbReco[globalindexmc] = 1; + fMCEventNbmctrack[globalindexmc] = midrap; + fNbMcEvent->Fill(0.5); + fMCCounters++; + + } + else { + fMCEventNbReco[globalindexmc] = fMCEventNbReco.find(globalindexmc)->second + 1; + } + + + } // end loop over events + + for (const auto& [mcEv, NbRecEv] : fMCEventNbReco) { + + fNbRecCollisionPerMCCollision->Fill(fMCEventNbReco.find(mcEv)->second); + + } + + for (const auto& [mcEv, NbRecEv] : fMCEventNbmctrack) { + + fNbmctrack->Fill(fMCEventNbReco.find(mcEv)->second,fMCEventNbmctrack.find(mcEv)->second); + + } + + } + + template + Int_t dNdetach(TTracksMC const& groupedMCTracks){ + + Int_t midrap = 0; + for(auto mctrack : groupedMCTracks){ + if(TMath::Abs(mctrack.eta())<0.5 && mctrack.isPhysicalPrimary() && (TMath::Abs(mctrack.pdgCode())==211 || mctrack.pdgCode()==111)) { + midrap++; + } + } + return midrap; + + } + + + void processSkimmed(soa::Filtered const&events, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + runSelection(events, eventsMC, tracksMC); + } + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisEventQa, processSkimmed, "Run event selection on DQ skimmed events", false); + PROCESS_SWITCH(AnalysisEventQa, processDummy, "Dummy process function", false); + +}; + +struct AnalysisTrackSelection { + + Produces trackSel; + Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; + + // configurables + //Configurable fConfigEventCuts{"cfgEventCuts", "eventStandardNoINT7", "Event selection"}; + Configurable fConfigCuts{"cfgTrackCuts", "jpsiPID1", "Comma separated list of barrel track cuts"}; + Configurable fConfigMCSignals{"cfgTrackMCSignals", "", "Comma separated list of MC signals"}; + + // 3D histos + Configurable fConfigUsePtVec{"cfgUsePtVec", true, "If true, non-linear pt bins predefined"}; + Configurable fConfigMinPt{"cfgMinPt", 0., "min Pt in 3D histos"}; + Configurable fConfigMaxPt{"cfgMaxPt", 10., "max Pt in 3D histos"}; + Configurable fConfigStepPt{"cfgStepPt", 1, "Nb of steps in pt in 3D histos"}; + Configurable fConfigMinEta{"cfgMinEta", -0.8, "min Eta in 3D histos"}; + Configurable fConfigMaxEta{"cfgMaxEta", 0.8, "max Eta in 3D histos"}; + Configurable fConfigStepEta{"cfgStepEta", 16, "Nb of steps in Eta in 3D histos"}; + Configurable fConfigMinPhi{"cfgMinPhi", 0., "min Phi in 3D histos"}; + Configurable fConfigMaxPhi{"cfgMaxPhi", 6.3, "max Phi in 3D histos"}; + Configurable fConfigStepPhi{"cfgStepPhi", 63, "Nb of steps in Phi in 3D histos"}; + + // Resolution histos + Configurable fConfigResolutionOn{"cfgResolution", false, "If true, fill resolution histograms"}; + Configurable fConfigUsePtVecRes{"cfgUsePtVecRes", true, "If true, non-linear pt bins predefined in res histos"}; + Configurable fConfigMinPtRes{"cfgMinPtRes", 0., "min Pt in res histos"}; + Configurable fConfigMaxPtRes{"cfgMaxPtRes", 20., "max Pt in res histos"}; + Configurable fConfigStepPtRes{"cfgStepPtRes", 1, "Nb of steps in pt in res histos"}; + Configurable fConfigStepDeltaPt{"cfgStepDeltaPt", 1, "Nb of steps in delta pt in res histos"}; + Configurable fConfigMinDeltaEta{"cfgMinDeltaEta", -0.5, "min delta Eta in res histos"}; + Configurable fConfigMaxDeltaEta{"cfgMaxDeltaEta", 0.5, "max delta Eta in res histos"}; + Configurable fConfigStepDeltaEta{"cfgStepDeltaEta", 500, "Nb of steps in detla Eta in res histos"}; + Configurable fConfigMinDeltaPhi{"cfgMinDeltaPhi", -0.5, "min delta Phi in res histos"}; + Configurable fConfigMaxDeltaPhi{"cfgMaxDeltaPhi", 0.5, "max delta Phi in res histos"}; + Configurable fConfigStepDeltaPhi{"cfgStepDeltaPhi", 500, "Nb of steps in delta Phi in res histos"}; + + + Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; + + // output lists + OutputObj fOutputList{"output"}; + THashList *fMainList; // Main list + THashList *fSingleElectronList; // 3D histos for MC and reconstructed signals + THashList *fResolutionList; // Resolution histograms + THashList *fQASingleElectronList; // QA in case on with histo manager outputs + + // Cuts and signals + //AnalysisCompositeCut* fEventCut; // Taken from event selection part + std::vector fTrackCuts; // list of track cuts + AnalysisCompositeCut* fTrackCutsRes; // track cut for resolution map + std::vector fMCSignals; // list of signals to be checked + MCSignal* fMCSignalRes; // signal for res + + //3D histos + std::vector fHistGenPosPart; + std::vector fHistGenNegPart; + std::vector fHistGenSmearedPosPart; + std::vector fHistGenSmearedNegPart; + std::vector fHistRecPosPart; + std::vector fHistRecNegPart; + // Binning + std::vector fPtBins; + std::vector fEtaBins; + std::vector fPhiBins; + + // Res histos + std::vector fHistRes; + // Binning + std::vector fPtResBins; + std::vector fDeltaEtaBins; + std::vector fDeltaPhiBins; + + // QA + HistogramManager* fHistManQA; // histo manager + std::vector fHistNamesRecoQA; // list of histo names for all reconstructed tracks in histo manager + std::vector> fHistNamesMCMatchedQA; // list of histo names for reconstructed signals in histo manager + std::vector> fHistNamesMCQA; // list of histo names for generated signals in histo manager + + void init(o2::framework::InitContext&) + { + + // Create list output + fMainList = new THashList; + fMainList->SetOwner(kTRUE); + fMainList->SetName("trackselection"); + + // Create list output for 3D eta,phi,pt + fSingleElectronList = new THashList; + fSingleElectronList->SetOwner(kTRUE); + fSingleElectronList->SetName("SingleElectron"); + + // Binning 3D histos + if(fConfigUsePtVec){ + const Int_t Npt = 68; + Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.}; + std::vector v_pte(pte,std::end(pte)); + fPtBins = v_pte; + } + else { + SetBinsLinear(fPtBins, fConfigMinPt, fConfigMaxPt, fConfigStepPt); + } + SetBinsLinear(fEtaBins, fConfigMinEta, fConfigMaxEta, fConfigStepEta); + SetBinsLinear(fPhiBins, fConfigMinPhi, fConfigMaxPhi, fConfigStepPhi); + const int fNptBins = fPtBins.size()-1; + const int fNetaBins = fEtaBins.size()-1; + const int fNphiBins = fPhiBins.size()-1; + + // Event cut: taken from event selection part + //fEventCut = new AnalysisCompositeCut(true); + //TString eventCutStr = fConfigEventCuts.value; + //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + + // List of track cuts + TString cutNamesStr = fConfigCuts.value; + if (!cutNamesStr.IsNull()) { + std::unique_ptr objArray(cutNamesStr.Tokenize(",")); + for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + fTrackCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); + } + } + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + VarManager::SetDefaultVarNames(); + + // List of MC signals + TString configSigNamesStr = fConfigMCSignals.value; + std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); + for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() != 1) { // NOTE: only 1 prong signals + continue; + } + // List of signal to be checked + fMCSignals.push_back(*sig); + } + } + + //Configure 3D histograms + // Create List with generated particles + TList* Generated = new TList(); + Generated->SetOwner(); + Generated->SetName("Generated"); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_pos->Sumw2(); + fHistGenPosPart.push_back(th3_tmp_pos); + Generated->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_neg->Sumw2(); + fHistGenNegPart.push_back(th3_tmp_neg); + Generated->Add(th3_tmp_neg); + } + // Create List with generated+smeared particles + TList* GeneratedSmeared = new TList(); + GeneratedSmeared->SetName("GeneratedSmeared"); + GeneratedSmeared->SetOwner(); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_pos->Sumw2(); + fHistGenSmearedPosPart.push_back(th3_tmp_pos); + GeneratedSmeared->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_neg->Sumw2(); + fHistGenSmearedNegPart.push_back(th3_tmp_neg); + GeneratedSmeared->Add(th3_tmp_neg); + } + + fSingleElectronList->Add(Generated); + fSingleElectronList->Add(GeneratedSmeared); + + // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge + for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ + TList* list = new TList(); + list->SetName(fTrackCuts.at(list_i).GetName()); + list->SetOwner(); + + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH3D* th3_tmp_pos = new TH3D(Form("Nrec_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_pos->Sumw2(); + th3_tmp_pos->SetDirectory(0x0); + fHistRecPosPart.push_back(th3_tmp_pos); + list->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Nrec_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + th3_tmp_neg->Sumw2(); + th3_tmp_neg->SetDirectory(0x0); + fHistRecNegPart.push_back(th3_tmp_neg); + list->Add(th3_tmp_neg); + + } + fSingleElectronList->Add(list); + } + fMainList->Add(fSingleElectronList); + + // Resolution histogramms + if(fConfigResolutionOn) { + + // Binning 3D histos + if(fConfigUsePtVecRes){ + const Int_t Npt = 73; + Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.,12.0,14.,16.,18.,20.}; + std::vector v_pte(pte,std::end(pte)); + fPtResBins = v_pte; + } + else { + SetBinsLinear(fPtResBins, fConfigMinPtRes, fConfigMaxPtRes, fConfigStepPtRes); + } + SetBinsLinear(fDeltaEtaBins, fConfigMinDeltaEta, fConfigMaxDeltaEta, fConfigStepDeltaEta); + SetBinsLinear(fDeltaPhiBins, fConfigMinDeltaPhi, fConfigMaxDeltaPhi, fConfigStepDeltaPhi); + const int fNptresBins = fPtResBins.size()-1; + const int fNDeltaetaBins = fDeltaEtaBins.size()-1; + const int fNDeltaphiBins = fDeltaPhiBins.size()-1; + + + fResolutionList = new THashList; + fResolutionList->SetOwner(kTRUE); + fResolutionList->SetName("Resolution"); + + printf("Histos\n"); + TH2D *thPtGen_DeltaPtOverPtGen = new TH2D("PtGen_DeltaPtOverPtGen", "", fNptresBins, fPtResBins.data(), fConfigStepDeltaPt, -1., +1.); + TH2D *thPtGen_DeltaEta = new TH2D("PtGen_DeltaEta", "", fNptresBins, fPtResBins.data(), fNDeltaetaBins, fDeltaEtaBins.data()); + TH2D *thPtGen_DeltaPhi_Ele = new TH2D("PtGen_DeltaPhi_Ele", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); + TH2D *thPtGen_DeltaPhi_Pos = new TH2D("PtGen_DeltaPhi_Pos", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); + + thPtGen_DeltaPtOverPtGen ->Sumw2(); + thPtGen_DeltaEta ->Sumw2(); + thPtGen_DeltaPhi_Ele ->Sumw2(); + thPtGen_DeltaPhi_Pos ->Sumw2(); + + thPtGen_DeltaPtOverPtGen ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPtOverPtGen ->GetYaxis()->SetTitle("(p^{gen}_{T} - p^{rec}_{T}) / p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaEta ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaEta ->GetYaxis()->SetTitle("#eta^{gen} - #eta^{rec}"); + thPtGen_DeltaPhi_Ele ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPhi_Ele ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); + thPtGen_DeltaPhi_Pos ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPhi_Pos ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); + + printf("Add\n"); + fHistRes.push_back(thPtGen_DeltaPtOverPtGen); + fHistRes.push_back(thPtGen_DeltaEta); + fHistRes.push_back(thPtGen_DeltaPhi_Ele); + fHistRes.push_back(thPtGen_DeltaPhi_Pos); + + fResolutionList->Add(thPtGen_DeltaPtOverPtGen); + fResolutionList->Add(thPtGen_DeltaEta); + fResolutionList->Add(thPtGen_DeltaPhi_Ele); + fResolutionList->Add(thPtGen_DeltaPhi_Pos); + + fMainList->Add(fResolutionList); + + } + + + // Configure QA histogram classes + if (fConfigQA){ + // Create list output for QA + fQASingleElectronList = new THashList; + fQASingleElectronList->SetOwner(kTRUE); + fQASingleElectronList->SetName("SingleElectronQA"); + TString histClassesQA = ""; + for (auto& cut : fTrackCuts) { + + // All reconstructed leptons + TString nameStr = Form("TrackBarrel_%s", cut.GetName()); + fHistNamesRecoQA.push_back(nameStr); + histClassesQA += Form("%s;", nameStr.Data()); + + // All reconstructed leptons matched to a 1 prong signal or MC 1 prong signal directly + std::vector mcnamesreco; + for (int isig = 0; isig < fMCSignals.size(); ++isig) { + TString nameStr2 = Form("TrackBarrel_%s_%s", cut.GetName(), fMCSignals.at(isig).GetName()); + mcnamesreco.push_back(nameStr2); + histClassesQA += Form("%s;", nameStr2.Data()); + } + fHistNamesMCMatchedQA.push_back(mcnamesreco); + } + + // Add histogram classes for each MC signal at generated level + std::vector mcnamesgen; + for (int isig = 0; isig < fMCSignals.size(); ++isig) { + TString nameStr2 = Form("MCTruthGen_%s", fMCSignals.at(isig).GetName()); + mcnamesgen.push_back(nameStr2); + histClassesQA += Form("%s;", nameStr2.Data()); + + } + fHistNamesMCQA.push_back(mcnamesgen); + + + fHistManQA = new HistogramManager("SingleElectronQA", "aa", VarManager::kNVars); + fHistManQA->SetUseDefaultVariableNames(kTRUE); + fHistManQA->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + DefineHistograms(fHistManQA, histClassesQA.Data()); // define all histograms + VarManager::SetUseVars(fHistManQA->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + fQASingleElectronList = fHistManQA->GetMainHistogramList(); + fMainList->Add(fQASingleElectronList); + + } + + fOutputList.setObject(fMainList); + +} + + Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; + + template + void runSelection(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) + { + + uint8_t eventFilter = 0; + std::map fMCEventLabels; + int fCounters = 0; //! [0] - particle counter, [1] - event counter + + + for (auto& event : events) { + + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::ResetValues(0, VarManager::kNMCParticleVariables); + // fill event information which might be needed in histograms that combine track and event properties + VarManager::FillEvent(event); + //if(!fEventCut->IsSelected(VarManager::fgValues)) continue; + eventFilter = uint32_t(event.isEventSelected()); + if(!eventFilter) continue; + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + VarManager::FillEvent(event.reducedMCevent()); + } + + // Look if we did not already saw the collision and fill the denominator of the single electron efficiency + Int_t globalindexmc = -1; + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + } + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fCounters; + fCounters++; + // skimmed data + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGenTrack(groupedMCTracks); + } + } + + + // Loop over reconstructed tracks belonging to the event and fill the numerator of the efficiency as well as the resolution map + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); + runRecTrack(groupedTracks,tracksMC); + } + + }// end loop over events + + } + + template + void runMCGenTrack(TTracksMC const& groupedMCTracks) + { + + + // loop over mc stack and fill histograms for pure MC truth signals + // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event + //auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); + for (auto& mctrack : groupedMCTracks) { + VarManager::FillTrack(mctrack); + // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. + // NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member. + // TODO: Use the mcReducedFlags to select signals + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + if ((*sig).CheckSignal(true, groupedMCTracks, mctrack)) { + if(mctrack.pdgCode() > 0) { + dynamic_cast(fHistGenNegPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); + } + else { + dynamic_cast(fHistGenPosPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); + } + if (fConfigQA) fHistManQA->FillHistClass(Form("MCTruthGen_%s", (*sig).GetName()), VarManager::fgValues); + } + } + } + } + + template + void runRecTrack(TTracks const& groupedTracks, TTracksMC const& tracksMC) + { + + uint32_t filterMap = 0; + trackSel.reserve(groupedTracks.size()); + for (auto& track : groupedTracks) { + filterMap = 0; + + VarManager::FillTrack(track); // compute track quantities + + // compute MC matched quantities + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + VarManager::FillTrack(track.reducedMCTrack()); + } + if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { + if (!track.has_mcParticle()) { + continue; + } + auto mctrack = track.template mcParticle_as(); + VarManager::FillTrack(mctrack); + } + + // compute track selection and publish the bit map + int i = 0; + for (auto cut = fTrackCuts.begin(); cut != fTrackCuts.end(); cut++, i++) { + if ((*cut).IsSelected(VarManager::fgValues)) { + filterMap |= (uint32_t(1) << i); + if (fConfigQA) { + fHistManQA->FillHistClass(fHistNamesRecoQA[i].Data(), VarManager::fgValues); + } + } + } + trackSel(static_cast(filterMap)); + if (!filterMap) { + continue; + } + + // compute MC matching decisions + uint32_t mcDecision = 0; + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + if ((*sig).CheckSignal(true, tracksMC, track.reducedMCTrack())) { + mcDecision |= (uint32_t(1) << isig); + } + } + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { + if (track.has_mcParticle()) { + auto mctrack = track.template mcParticle_as(); + if ((*sig).CheckSignal(true, tracksMC, mctrack)) { + mcDecision |= (uint32_t(1) << isig); + } + } + } + + } + + + // fill histograms + for (unsigned int i = 0; i < fMCSignals.size(); i++) { + if (!(mcDecision & (uint32_t(1) << i))) { + continue; + } + for (unsigned int j = 0; j < fTrackCuts.size(); j++) { + if (filterMap & (uint8_t(1) << j)) { + if(track.sign() < 0) { + dynamic_cast(fHistRecNegPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); + } + else { + dynamic_cast(fHistRecPosPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); + } + if(fConfigResolutionOn && (i==0) && (j==0)){ + + Double_t mcpt = -10000.; + Double_t mceta = -10000.; + Double_t mcphi = -1000.; + Int_t mcpdg = -10000.; + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto mctrack = track.reducedMCTrack(); + mcpt = mctrack.pt(); + mceta = mctrack.pt(); + mcphi = mctrack.pt(); + mcpdg = mctrack.pdgCode(); + } + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { + if (track.has_mcParticle()) { + auto mctrack = track.template mcParticle_as(); + mcpt = mctrack.pt(); + mceta = mctrack.pt(); + mcphi = mctrack.pt(); + mcpdg = mctrack.pdgCode(); + } + } + + + Double_t deltaptoverpt = -1000.; + if(mcpt>0.) deltaptoverpt = (mcpt-track.pt())/mcpt; + Double_t deltaeta = mceta-track.eta(); + Double_t deltaphi = mcphi-track.phi(); + dynamic_cast(fHistRes.at(0))->Fill(mcpt,deltaptoverpt); + dynamic_cast(fHistRes.at(1))->Fill(mcpt,deltaeta); + if(mcpdg<0){ + dynamic_cast(fHistRes.at(2))->Fill(mcpt,deltaphi); + }else { + dynamic_cast(fHistRes.at(3))->Fill(mcpt,deltaphi); + } + } + if(fConfigQA) fHistManQA->FillHistClass(fHistNamesMCMatchedQA[j][i].Data(), VarManager::fgValues); + } + } // end loop over cuts + } // end loop over MC signals + + + } // end loop over reconstructed track belonging to the events + + } + + + + void processSkimmed(soa::Filtered const& events, MyBarrelTracks const& tracks, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + runSelection(events, tracks, eventsMC, tracksMC); + + + } + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisTrackSelection, processSkimmed, "Run barrel track selection on DQ skimmed tracks", false); + PROCESS_SWITCH(AnalysisTrackSelection, processDummy, "Dummy process function", false); +}; + + +struct AnalysisSameEventPairing { + + + //Produces dileptonList; + + // Filter based on previous components in the task + Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; + Filter filterBarrelTrackSelected = aod::emanalysisflags::isBarrelSelected > 0; + + Configurable fConfigTrackCuts{"cfgTrackCuts", "", "Comma separated list of barrel track cuts"}; + Configurable fConfigMCSignals{"cfgBarrelMCSignals", "", "Comma separated list of MC signals"}; + Configurable fConfigMinPt{"cfgMinPt", 0., "Fiducial min Pt for MC signal"}; + Configurable fConfigMaxPt{"cfgMaxPt", 10., "Fiducial max Pt for MC signal"}; + Configurable fConfigMinEta{"cfgMinEta", -0.8, "Fiducial min eta for MC signal"}; + Configurable fConfigMaxEta{"cfgMaxEta", 0.8, "Fiducial max eta for MC signal"}; + Configurable fConfigFlatTables{"cfgFlatTables", false, "Produce a single flat tables with all relevant information of the pairs and single tracks"}; + + // TODO: here we specify signals, however signal decisions are precomputed and stored in mcReducedFlags + // TODO: The tasks based on skimmed MC could/should rely ideally just on these flags + // TODO: special AnalysisCuts to be prepared in this direction + // TODO: cuts on the MC truth information to be added if needed + + + // 2D histos: mee and ptee + Configurable fConfigMinPtee{"cfgMinPtee", 0., "min Ptee in 2D histos"}; + Configurable fConfigMaxPtee{"cfgMaxPtee", 10., "max Ptee in 2D histos"}; + Configurable fConfigStepPtee{"cfgStepPtee", 100, "Nb of steps in ptee in 2D histos"}; + Configurable fConfigMinMee{"cfgMinMee", 0., "min Mee in 2D histos"}; + Configurable fConfigMaxMee{"cfgMaxMee", 3.5, "max Mee in 2D histos"}; + Configurable fConfigStepMee{"cfgStepMee", 600, "Nb of steps in Mee in 2D histos"}; + + // output lists + OutputObj fOutputList{"output"}; + THashList *fMainList; // Main list + THashList *fPairList; // 2D histos for MC and reconstructed signals + THashList *fQAPairList; // QA in case on with histo manager outputs + + // Cuts and signals + //AnalysisCompositeCut* fEventCut; // Taken from event selection part + std::vector fTrackCuts; // list of track cuts + std::vector fMCSignals; // list of signals with one prong to be checked: ULS 2D histos + + //2D histo vectors + std::vector fHistGenPair; + std::vector fHistGenSmearedPair; + std::vector fHistRecPair; + // Binning + std::vector fPteeBins; + std::vector fMeeBins; + + // QA: to be defined + + + void init(o2::framework::InitContext& context) + { + + printf("Start init\n"); + + // Create list output + fMainList = new THashList; + fMainList->SetOwner(kTRUE); + fMainList->SetName("pairselection"); + + // Create list output for 3D eta,phi,pt + fPairList = new THashList; + fPairList->SetOwner(kTRUE); + fPairList->SetName("Dielectron"); + + // Binning 2D histos + SetBinsLinear(fPteeBins, fConfigMinPtee, fConfigMaxPtee, fConfigStepPtee); + SetBinsLinear(fMeeBins, fConfigMinMee, fConfigMaxMee, fConfigStepMee); + const int fNpteeBins = fPteeBins.size()-1; + const int fNmeeBins = fMeeBins.size()-1; + + // Event cut: taken from event selection part + //fEventCut = new AnalysisCompositeCut(true); + //TString eventCutStr = fConfigEventCuts.value; + //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + + // List of track cuts + TString cutNamesStr = fConfigTrackCuts.value; + if (!cutNamesStr.IsNull()) { + std::unique_ptr objArray(cutNamesStr.Tokenize(",")); + for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + fTrackCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); + } + } + VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + VarManager::SetDefaultVarNames(); + + // List of MC signals + TString configSigNamesStr = fConfigMCSignals.value; + std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); + for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() == 2) { // only 2 prong signals + fMCSignals.push_back(*sig); + } + // List of signal to be checked + } + } + + //Configure 2D histograms + // Create List with generated particles + TList* Generated = new TList(); + Generated->SetOwner(); + Generated->SetName("Generated"); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH2D* th2_tmp = new TH2D(Form("Ngen_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + th2_tmp->Sumw2(); + fHistGenPair.push_back(th2_tmp); + Generated->Add(th2_tmp); + } + + // Create List with generated+smeared particles + TList* GeneratedSmeared = new TList(); + GeneratedSmeared->SetName("GeneratedSmeared"); + GeneratedSmeared->SetOwner(); + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH2D* th2_tmp = new TH2D(Form("NgenSmeared_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + th2_tmp->Sumw2(); + fHistGenSmearedPair.push_back(th2_tmp); + GeneratedSmeared->Add(th2_tmp); + } + + + fPairList->Add(Generated); + fPairList->Add(GeneratedSmeared); + + // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge + for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ + TList* list = new TList(); + list->SetName(fTrackCuts.at(list_i).GetName()); + list->SetOwner(); + + for (unsigned int i = 0; i < fMCSignals.size(); ++i){ + TH2D* th2_tmp = new TH2D(Form("Nrec_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + th2_tmp->Sumw2(); + th2_tmp->SetDirectory(0x0); + fHistRecPair.push_back(th2_tmp); + list->Add(th2_tmp); + } + fPairList->Add(list); + } + fMainList->Add(fPairList); + + fOutputList.setObject(fMainList); + + } + + Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; + + template < uint32_t TEventFillMap, uint32_t TEventMCFillMap, uint32_t TTrackFillMap, typename TEvents, typename TTracks, typename TEventsMC, typename TTracksMC> + void runPairing(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) + { + + + std::map fMCEventLabels; + int fCounters = 0; //! [0] - particle counter, [1] - event counter + + + for (auto& event : events) { + + if (!event.isEventSelected()) { + return; + } + + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::ResetValues(0, VarManager::kNMCParticleVariables); + // fill event information which might be needed in histograms that combine track and event properties + VarManager::FillEvent(event); + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + VarManager::FillEvent(event.reducedMCevent()); + } + + // Look if we did not already saw the collision and fill the denominator of the single electron efficiency + Int_t globalindexmc = -1; + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + } + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fCounters; + fCounters++; + // skimmed data + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGenPair(groupedMCTracks); + } + } + + + //auto groupedTrackspos = posTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); + //groupedTrackspos.bindInternalIndicesTo(&tracks); + //auto groupedTracksneg = negTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); + //groupedTracksneg.bindInternalIndicesTo(&tracks); + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); + groupedTracks.bindInternalIndicesTo(&tracks); + runRecPair(groupedTracks,tracksMC); + } + + } // end loop over reconstructed event + } // end loop pairing function + + + template + void runMCGenPair(TTracksMC const& groupedMCTracks) + { + // + Double_t masse = 0.00051099895; // 0.5 MeV/c2 -> 0.0005 GeV/c2 + + for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { + + if((abs(t1.pdgCode()) != 11) || (abs(t2.pdgCode()) != 11)) continue; + if(t1.pdgCode()*t2.pdgCode() > 0) continue; // ULS only + + TLorentzVector Lvec1; + TLorentzVector Lvec2; + Lvec1.SetPtEtaPhiM(t1.pt(),t1.eta(),t1.phi(),masse); + Lvec2.SetPtEtaPhiM(t2.pt(),t2.eta(),t2.phi(),masse); + TLorentzVector LvecM = Lvec1 + Lvec2; + double mass = LvecM.M(); + double pairpt = LvecM.Pt(); + //double opangle = Lvec1.Angle(Lvec2.Vect()); + + //printf("Check before\n"); + // Fiducial cut + Bool_t genfidcut = kTRUE; + if((t1.eta()> fConfigMaxEta) || (t2.eta() > fConfigMaxEta) || (t1.eta() < fConfigMinEta) || (t2.eta() < fConfigMinEta) || (t1.pt()> fConfigMaxPt) || (t2.pt() > fConfigMaxPt) || (t1.pt() < fConfigMinPt) || (t2.pt() < fConfigMinPt)) genfidcut = kFALSE; + + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + if ((*sig).CheckSignal(true, groupedMCTracks, t1, t2)) { + + // not smeared + if(genfidcut) dynamic_cast(fHistGenPair.at(isig))->Fill(mass,pairpt); + // need to implement smeared + + } + } + } //end of true pairing loop + } // end runMCGen + + + + template < uint32_t TTrackFillMap, typename TTracks, typename TTracksMC> + void runRecPair(TTracks const& tracks, TTracksMC const& tracksMC) + { + + // Loop over two track combinations + uint8_t twoTrackFilter = 0; + //uint32_t dileptonFilterMap = 0; + // uint32_t dileptonMcDecision = 0; + //dileptonList.reserve(1); + + + + for (auto& [t1, t2] : combinations(tracks, tracks)) { + + twoTrackFilter = uint32_t(t1.isBarrelSelected()) & uint32_t(t2.isBarrelSelected()); + + if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + continue; + } + VarManager::FillPair(t1, t2); + + // run MC matching for this pair + uint32_t mcDecision = 0; + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { // for skimmed DQ model + if ((*sig).CheckSignal(true, tracksMC, t1.reducedMCTrack(), t2.reducedMCTrack())) { + mcDecision |= (uint32_t(1) << isig); + } + } + if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { // for Ao2d + if (t1.has_mcParticle() && t2.has_mcParticle()) { + auto mctrack1 = t1.template mcParticle_as(); + auto mctrack2 = t2.template mcParticle_as(); + if ((*sig).CheckSignal(true, tracksMC, mctrack1, mctrack2)) { + mcDecision |= (uint32_t(1) << isig); + } + } + } + + }// end of loop MC signals + + //dileptonFilterMap = twoTrackFilter; + //dileptonMcDecision = mcDecision; + //dileptonList(event, VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), dileptonFilterMap, dileptonMcDecision); + + + for (unsigned int i = 0; i < fMCSignals.size(); i++) { + if (!(mcDecision & (uint32_t(1) << i))) { + continue; + } + for (unsigned int j = 0; j < fTrackCuts.size(); j++) { + if (twoTrackFilter & (uint8_t(1) << j)) { + dynamic_cast(fHistRecPair.at(j * fMCSignals.size() + i))->Fill(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt]); + } + } + } + } + + } + + void processToEESkimmed(soa::Filtered const& events, + soa::Filtered const& tracks, + ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + { + + runPairing< gkEventFillMap, gkMCEventFillMap, gkTrackFillMap>(events, tracks, eventsMC, tracksMC); + + } + + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisSameEventPairing, processToEESkimmed, "Run barrel barrel pairing on DQ skimmed tracks", false); + PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy process function", false); +}; + + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +} + +void DefineHistograms(HistogramManager* histMan, TString histClasses) +{ + // + // Define here the histograms for all the classes required in analysis. + // The histogram classes are provided in the histClasses string, separated by semicolon ";" + // The histogram classes and their components histograms are defined below depending on the name of the histogram class + // + std::unique_ptr objArray(histClasses.Tokenize(";")); + for (Int_t iclass = 0; iclass < objArray->GetEntries(); ++iclass) { + TString classStr = objArray->At(iclass)->GetName(); + histMan->AddHistClass(classStr.Data()); + + // NOTE: The level of detail for histogramming can be controlled via configurables + if (classStr.Contains("Event")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", "trigger,cent,mc"); + } + + if (classStr.Contains("Track")) { + if (classStr.Contains("Barrel")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", "its,tpcpid,dca,tofpid,mc"); + } + } + + if (classStr.Contains("Pairs")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel", "vertexing-barrel"); + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_dimuon", "vertexing-forward"); + } + + if (classStr.Contains("MCTruthGenPair")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair"); + histMan->AddHistogram(objArray->At(iclass)->GetName(), "Pt", "MC generator p_{T} distribution", false, 200, 0.0, 20.0, VarManager::kMCPt); + histMan->AddHistogram(objArray->At(iclass)->GetName(), "Eta", "MC generator #eta distribution", false, 500, -5.0, 5.0, VarManager::kMCEta); + histMan->AddHistogram(objArray->At(iclass)->GetName(), "Phi", "MC generator #varphi distribution", false, 500, -6.3, 6.3, VarManager::kMCPhi); + } + if (classStr.Contains("MCTruthGen")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth"); + } + if (classStr.Contains("DileptonsSelected")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel"); + } + + + } +}// end loop over histogram classes + +void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps){ + fBins.clear(); + const double stepSize = (max - min) / steps; + for (unsigned int i = 0; i < steps+1; ++i){ + fBins.push_back(i * stepSize + min); + } +} From 836b4239cb391fb96a27a299daff09f128c24618 Mon Sep 17 00:00:00 2001 From: Raphaelle Bailhache Date: Thu, 14 Jul 2022 22:13:49 +0200 Subject: [PATCH 2/4] Clear up --- PWGEM/Dilepton/Tasks/CMakeLists.txt~ | 45 - PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ | 1189 ---------------------- 2 files changed, 1234 deletions(-) delete mode 100644 PWGEM/Dilepton/Tasks/CMakeLists.txt~ delete mode 100644 PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt~ b/PWGEM/Dilepton/Tasks/CMakeLists.txt~ deleted file mode 100644 index 56e990ba20b..00000000000 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt~ +++ /dev/null @@ -1,45 +0,0 @@ -# Copyright 2019-2020 CERN and copyright holders of ALICE O2. -# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -# All rights not expressly granted are reserved. -# -# This software is distributed under the terms of the GNU General Public -# License v3 (GPL Version 3), copied verbatim in the file "COPYING". -# -# In applying this license CERN does not waive the privileges and immunities -# granted to it by virtue of its status as an Intergovernmental Organization -# or submit itself to any jurisdiction. - -o2physics_add_dpl_workflow(table-reader - SOURCES tableReader.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(efficiency - SOURCES dqEfficiency.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(dilepton-ee - SOURCES dileptonEE.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(dilepton-mumu - SOURCES dileptonMuMu.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(filter-pp - SOURCES filterPP.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(v0-selector - SOURCES v0selector.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore O2::DetectorsVertexing - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(flow - SOURCES dqFlow.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore - COMPONENT_NAME Analysis) diff --git a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ deleted file mode 100644 index 50c71d4eeab..00000000000 --- a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx~ +++ /dev/null @@ -1,1189 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. -// -// -// Analysis task for calculating single electron and dielectron efficiency -// -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "PWGDQ/DataModel/ReducedInfoTables.h" -#include "PWGDQ/Core/VarManager.h" -#include "PWGDQ/Core/HistogramManager.h" -#include "PWGDQ/Core/AnalysisCut.h" -#include "PWGDQ/Core/AnalysisCompositeCut.h" -#include "PWGDQ/Core/HistogramsLibrary.h" -#include "PWGDQ/Core/CutsLibrary.h" -#include "PWGDQ/Core/MCSignal.h" -#include "PWGDQ/Core/MCSignalLibrary.h" -#include -#include -#include -#include -#include -#include -#include - -using std::cout; -using std::endl; -using std::string; - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::aod; - -// Some definitions -namespace o2::aod -{ - -namespace emanalysisflags -{ - DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); - DECLARE_SOA_COLUMN(IsBarrelSelected, isBarrelSelected, int); -} // namespace reducedevent - - DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", emanalysisflags::IsEventSelected); - DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "BARRELTRACKCUTS", emanalysisflags::IsBarrelSelected); -} // namespace o2::aod - -// -using MyEvents = soa::Join; -using MyEventsSelected = soa::Join; -using MyBarrelTracks = soa::Join; -//using MyBarrelTracksWithCov = soa::Join; -using MyBarrelTracksSelected = soa::Join; -//using MyBarrelTracksSelectedWithCov = soa::Join; - - -// -constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended; -constexpr static uint32_t gkMCEventFillMap = VarManager::ObjTypes::ReducedEventMC; -constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; -//constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; -constexpr static uint32_t gkParticleMCFillMap = VarManager::ObjTypes::ParticleMC; - - - -void DefineHistograms(HistogramManager* histMan, TString histClasses); -void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps); - -struct AnalysisEventSelection { - - Produces eventSel; - OutputObj fOutputList{"output"}; - Configurable fConfigEventCuts{"cfgEventCuts", "eventStandard", "Event selection"}; - Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; - - HistogramManager* fHistMan; - AnalysisCompositeCut* fEventCut; - - void init(o2::framework::InitContext&) - { - fEventCut = new AnalysisCompositeCut(true); - TString eventCutStr = fConfigEventCuts.value; - fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); - VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill - - VarManager::SetDefaultVarNames(); - if (fConfigQA) { - fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); - fHistMan->SetUseDefaultVariableNames(kTRUE); - fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); - DefineHistograms(fHistMan, "Event_BeforeCuts;Event_AfterCuts;"); // define all histograms - VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill - fOutputList.setObject(fHistMan->GetMainHistogramList()); - } - } - - template - void runSelection(TEvent const& event, TEventsMC const& mcEvents) - { - // Reset the values array - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); - - VarManager::FillEvent(event); - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - VarManager::FillEvent(event.reducedMCevent()); - } - if (fConfigQA) { - fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); // automatically fill all the histograms in the class Event - } - if (fEventCut->IsSelected(VarManager::fgValues)) { - if (fConfigQA) { - fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); - } - eventSel(1); - } else { - eventSel(0); - } - } - - void processSkimmed(MyEvents::iterator const& event, aod::ReducedMCEvents const& mcEvents) - { - runSelection(event, mcEvents); - } - void processDummy(MyEvents&) - { - // do nothing - } - - PROCESS_SWITCH(AnalysisEventSelection, processSkimmed, "Run event selection on DQ skimmed events", false); - PROCESS_SWITCH(AnalysisEventSelection, processDummy, "Dummy process function", false); -}; -struct AnalysisEventQa { - - Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; - OutputObj fOutputList{"output"}; - THashList *fMainList; - TH1D *fNbRecCollisionPerMCCollision; - TProfile *fNbmctrack; - TH1D *fNbMcEvent; - TH1D *fNbRecEvent; - - void init(o2::framework::InitContext&) - { - - VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill - VarManager::SetDefaultVarNames(); - - fMainList = new THashList; - fMainList->SetOwner(kTRUE); - fMainList->SetName("OAnalysisEventQA"); - fNbRecCollisionPerMCCollision = new TH1D("MCEvRecEv","",20,0.,20.); - fNbmctrack = new TProfile("mctrack","",20,0.,20.); - fNbMcEvent = new TH1D("MCEvent","",1,0.,1.); - fNbRecEvent = new TH1D("RecEvent","",1,0.,1.); - fMainList->Add(fNbRecCollisionPerMCCollision); - fMainList->Add(fNbmctrack); - fMainList->Add(fNbMcEvent); - fMainList->Add(fNbRecEvent); - fOutputList.setObject(fMainList); - - } - - Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; - - template - void runSelection(TEvents const& events, TEventsMC const& eventsMC, TTracksMC const& tracksMC) - { - - uint8_t eventFilter = 0; - std::map fMCEventNbmctrack; - std::map fMCEventNbReco; - std::map fMCEventLabels; - - - int fMCCounters = 0; //! - mc event counter - int fEvCounters = 0; // ! - rec event counter - - - // First loop - - for (auto& event : events) { - - - // Reset the values array - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); - VarManager::FillEvent(event); - eventFilter = uint32_t(event.isEventSelected()); - if(!eventFilter) continue; - fEvCounters++; - fNbRecEvent->Fill(0.5); - - - Int_t midrap = 0; - Int_t globalindexmc = -1; - - // skimmed data - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - - auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); - midrap = dNdetach(groupedMCTracks); - - auto mcEvent = event.reducedMCevent(); - globalindexmc = mcEvent.globalIndex(); - - } - - if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { - fMCEventLabels[globalindexmc] = fMCCounters; - fMCEventNbReco[globalindexmc] = 1; - fMCEventNbmctrack[globalindexmc] = midrap; - fNbMcEvent->Fill(0.5); - fMCCounters++; - - } - else { - fMCEventNbReco[globalindexmc] = fMCEventNbReco.find(globalindexmc)->second + 1; - } - - - } // end loop over events - - for (const auto& [mcEv, NbRecEv] : fMCEventNbReco) { - - fNbRecCollisionPerMCCollision->Fill(fMCEventNbReco.find(mcEv)->second); - - } - - for (const auto& [mcEv, NbRecEv] : fMCEventNbmctrack) { - - fNbmctrack->Fill(fMCEventNbReco.find(mcEv)->second,fMCEventNbmctrack.find(mcEv)->second); - - } - - } - - template - Int_t dNdetach(TTracksMC const& groupedMCTracks){ - - Int_t midrap = 0; - for(auto mctrack : groupedMCTracks){ - if(TMath::Abs(mctrack.eta())<0.5 && mctrack.isPhysicalPrimary() && (TMath::Abs(mctrack.pdgCode())==211 || mctrack.pdgCode()==111)) { - midrap++; - } - } - return midrap; - - } - - - void processSkimmed(soa::Filtered const&events, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) - { - runSelection(events, eventsMC, tracksMC); - } - void processDummy(MyEvents&) - { - // do nothing - } - - PROCESS_SWITCH(AnalysisEventQa, processSkimmed, "Run event selection on DQ skimmed events", false); - PROCESS_SWITCH(AnalysisEventQa, processDummy, "Dummy process function", false); - -}; - -struct AnalysisTrackSelection { - - Produces trackSel; - Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; - - // configurables - //Configurable fConfigEventCuts{"cfgEventCuts", "eventStandardNoINT7", "Event selection"}; - Configurable fConfigCuts{"cfgTrackCuts", "jpsiPID1", "Comma separated list of barrel track cuts"}; - Configurable fConfigMCSignals{"cfgTrackMCSignals", "", "Comma separated list of MC signals"}; - - // 3D histos - Configurable fConfigUsePtVec{"cfgUsePtVec", true, "If true, non-linear pt bins predefined"}; - Configurable fConfigMinPt{"cfgMinPt", 0., "min Pt in 3D histos"}; - Configurable fConfigMaxPt{"cfgMaxPt", 10., "max Pt in 3D histos"}; - Configurable fConfigStepPt{"cfgStepPt", 1, "Nb of steps in pt in 3D histos"}; - Configurable fConfigMinEta{"cfgMinEta", -0.8, "min Eta in 3D histos"}; - Configurable fConfigMaxEta{"cfgMaxEta", 0.8, "max Eta in 3D histos"}; - Configurable fConfigStepEta{"cfgStepEta", 16, "Nb of steps in Eta in 3D histos"}; - Configurable fConfigMinPhi{"cfgMinPhi", 0., "min Phi in 3D histos"}; - Configurable fConfigMaxPhi{"cfgMaxPhi", 6.3, "max Phi in 3D histos"}; - Configurable fConfigStepPhi{"cfgStepPhi", 63, "Nb of steps in Phi in 3D histos"}; - - // Resolution histos - Configurable fConfigResolutionOn{"cfgResolution", false, "If true, fill resolution histograms"}; - Configurable fConfigUsePtVecRes{"cfgUsePtVecRes", true, "If true, non-linear pt bins predefined in res histos"}; - Configurable fConfigMinPtRes{"cfgMinPtRes", 0., "min Pt in res histos"}; - Configurable fConfigMaxPtRes{"cfgMaxPtRes", 20., "max Pt in res histos"}; - Configurable fConfigStepPtRes{"cfgStepPtRes", 1, "Nb of steps in pt in res histos"}; - Configurable fConfigStepDeltaPt{"cfgStepDeltaPt", 1, "Nb of steps in delta pt in res histos"}; - Configurable fConfigMinDeltaEta{"cfgMinDeltaEta", -0.5, "min delta Eta in res histos"}; - Configurable fConfigMaxDeltaEta{"cfgMaxDeltaEta", 0.5, "max delta Eta in res histos"}; - Configurable fConfigStepDeltaEta{"cfgStepDeltaEta", 500, "Nb of steps in detla Eta in res histos"}; - Configurable fConfigMinDeltaPhi{"cfgMinDeltaPhi", -0.5, "min delta Phi in res histos"}; - Configurable fConfigMaxDeltaPhi{"cfgMaxDeltaPhi", 0.5, "max delta Phi in res histos"}; - Configurable fConfigStepDeltaPhi{"cfgStepDeltaPhi", 500, "Nb of steps in delta Phi in res histos"}; - - - Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; - - // output lists - OutputObj fOutputList{"output"}; - THashList *fMainList; // Main list - THashList *fSingleElectronList; // 3D histos for MC and reconstructed signals - THashList *fResolutionList; // Resolution histograms - THashList *fQASingleElectronList; // QA in case on with histo manager outputs - - // Cuts and signals - //AnalysisCompositeCut* fEventCut; // Taken from event selection part - std::vector fTrackCuts; // list of track cuts - AnalysisCompositeCut* fTrackCutsRes; // track cut for resolution map - std::vector fMCSignals; // list of signals to be checked - MCSignal* fMCSignalRes; // signal for res - - //3D histos - std::vector fHistGenPosPart; - std::vector fHistGenNegPart; - std::vector fHistGenSmearedPosPart; - std::vector fHistGenSmearedNegPart; - std::vector fHistRecPosPart; - std::vector fHistRecNegPart; - // Binning - std::vector fPtBins; - std::vector fEtaBins; - std::vector fPhiBins; - - // Res histos - std::vector fHistRes; - // Binning - std::vector fPtResBins; - std::vector fDeltaEtaBins; - std::vector fDeltaPhiBins; - - // QA - HistogramManager* fHistManQA; // histo manager - std::vector fHistNamesRecoQA; // list of histo names for all reconstructed tracks in histo manager - std::vector> fHistNamesMCMatchedQA; // list of histo names for reconstructed signals in histo manager - std::vector> fHistNamesMCQA; // list of histo names for generated signals in histo manager - - void init(o2::framework::InitContext&) - { - - // Create list output - fMainList = new THashList; - fMainList->SetOwner(kTRUE); - fMainList->SetName("trackselection"); - - // Create list output for 3D eta,phi,pt - fSingleElectronList = new THashList; - fSingleElectronList->SetOwner(kTRUE); - fSingleElectronList->SetName("SingleElectron"); - - // Binning 3D histos - if(fConfigUsePtVec){ - const Int_t Npt = 68; - Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.}; - std::vector v_pte(pte,std::end(pte)); - fPtBins = v_pte; - } - else { - SetBinsLinear(fPtBins, fConfigMinPt, fConfigMaxPt, fConfigStepPt); - } - SetBinsLinear(fEtaBins, fConfigMinEta, fConfigMaxEta, fConfigStepEta); - SetBinsLinear(fPhiBins, fConfigMinPhi, fConfigMaxPhi, fConfigStepPhi); - const int fNptBins = fPtBins.size()-1; - const int fNetaBins = fEtaBins.size()-1; - const int fNphiBins = fPhiBins.size()-1; - - // Event cut: taken from event selection part - //fEventCut = new AnalysisCompositeCut(true); - //TString eventCutStr = fConfigEventCuts.value; - //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); - - // List of track cuts - TString cutNamesStr = fConfigCuts.value; - if (!cutNamesStr.IsNull()) { - std::unique_ptr objArray(cutNamesStr.Tokenize(",")); - for (int icut = 0; icut < objArray->GetEntries(); ++icut) { - fTrackCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); - } - } - VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill - VarManager::SetDefaultVarNames(); - - // List of MC signals - TString configSigNamesStr = fConfigMCSignals.value; - std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); - for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { - MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); - if (sig) { - if (sig->GetNProngs() != 1) { // NOTE: only 1 prong signals - continue; - } - // List of signal to be checked - fMCSignals.push_back(*sig); - } - } - - //Configure 3D histograms - // Create List with generated particles - TList* Generated = new TList(); - Generated->SetOwner(); - Generated->SetName("Generated"); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_pos->Sumw2(); - fHistGenPosPart.push_back(th3_tmp_pos); - Generated->Add(th3_tmp_pos); - TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_neg->Sumw2(); - fHistGenNegPart.push_back(th3_tmp_neg); - Generated->Add(th3_tmp_neg); - } - // Create List with generated+smeared particles - TList* GeneratedSmeared = new TList(); - GeneratedSmeared->SetName("GeneratedSmeared"); - GeneratedSmeared->SetOwner(); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_pos->Sumw2(); - fHistGenSmearedPosPart.push_back(th3_tmp_pos); - GeneratedSmeared->Add(th3_tmp_pos); - TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_neg->Sumw2(); - fHistGenSmearedNegPart.push_back(th3_tmp_neg); - GeneratedSmeared->Add(th3_tmp_neg); - } - - fSingleElectronList->Add(Generated); - fSingleElectronList->Add(GeneratedSmeared); - - // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge - for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ - TList* list = new TList(); - list->SetName(fTrackCuts.at(list_i).GetName()); - list->SetOwner(); - - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH3D* th3_tmp_pos = new TH3D(Form("Nrec_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_pos->Sumw2(); - th3_tmp_pos->SetDirectory(0x0); - fHistRecPosPart.push_back(th3_tmp_pos); - list->Add(th3_tmp_pos); - TH3D* th3_tmp_neg = new TH3D(Form("Nrec_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_neg->Sumw2(); - th3_tmp_neg->SetDirectory(0x0); - fHistRecNegPart.push_back(th3_tmp_neg); - list->Add(th3_tmp_neg); - - } - fSingleElectronList->Add(list); - } - fMainList->Add(fSingleElectronList); - - // Resolution histogramms - if(fConfigResolutionOn) { - - // Binning 3D histos - if(fConfigUsePtVecRes){ - const Int_t Npt = 73; - Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.,12.0,14.,16.,18.,20.}; - std::vector v_pte(pte,std::end(pte)); - fPtResBins = v_pte; - } - else { - SetBinsLinear(fPtResBins, fConfigMinPtRes, fConfigMaxPtRes, fConfigStepPtRes); - } - SetBinsLinear(fDeltaEtaBins, fConfigMinDeltaEta, fConfigMaxDeltaEta, fConfigStepDeltaEta); - SetBinsLinear(fDeltaPhiBins, fConfigMinDeltaPhi, fConfigMaxDeltaPhi, fConfigStepDeltaPhi); - const int fNptresBins = fPtResBins.size()-1; - const int fNDeltaetaBins = fDeltaEtaBins.size()-1; - const int fNDeltaphiBins = fDeltaPhiBins.size()-1; - - - fResolutionList = new THashList; - fResolutionList->SetOwner(kTRUE); - fResolutionList->SetName("Resolution"); - - printf("Histos\n"); - TH2D *thPtGen_DeltaPtOverPtGen = new TH2D("PtGen_DeltaPtOverPtGen", "", fNptresBins, fPtResBins.data(), fConfigStepDeltaPt, -1., +1.); - TH2D *thPtGen_DeltaEta = new TH2D("PtGen_DeltaEta", "", fNptresBins, fPtResBins.data(), fNDeltaetaBins, fDeltaEtaBins.data()); - TH2D *thPtGen_DeltaPhi_Ele = new TH2D("PtGen_DeltaPhi_Ele", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); - TH2D *thPtGen_DeltaPhi_Pos = new TH2D("PtGen_DeltaPhi_Pos", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); - - thPtGen_DeltaPtOverPtGen ->Sumw2(); - thPtGen_DeltaEta ->Sumw2(); - thPtGen_DeltaPhi_Ele ->Sumw2(); - thPtGen_DeltaPhi_Pos ->Sumw2(); - - thPtGen_DeltaPtOverPtGen ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaPtOverPtGen ->GetYaxis()->SetTitle("(p^{gen}_{T} - p^{rec}_{T}) / p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaEta ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaEta ->GetYaxis()->SetTitle("#eta^{gen} - #eta^{rec}"); - thPtGen_DeltaPhi_Ele ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaPhi_Ele ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); - thPtGen_DeltaPhi_Pos ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaPhi_Pos ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); - - printf("Add\n"); - fHistRes.push_back(thPtGen_DeltaPtOverPtGen); - fHistRes.push_back(thPtGen_DeltaEta); - fHistRes.push_back(thPtGen_DeltaPhi_Ele); - fHistRes.push_back(thPtGen_DeltaPhi_Pos); - - fResolutionList->Add(thPtGen_DeltaPtOverPtGen); - fResolutionList->Add(thPtGen_DeltaEta); - fResolutionList->Add(thPtGen_DeltaPhi_Ele); - fResolutionList->Add(thPtGen_DeltaPhi_Pos); - - fMainList->Add(fResolutionList); - - } - - - // Configure QA histogram classes - if (fConfigQA){ - // Create list output for QA - fQASingleElectronList = new THashList; - fQASingleElectronList->SetOwner(kTRUE); - fQASingleElectronList->SetName("SingleElectronQA"); - TString histClassesQA = ""; - for (auto& cut : fTrackCuts) { - - // All reconstructed leptons - TString nameStr = Form("TrackBarrel_%s", cut.GetName()); - fHistNamesRecoQA.push_back(nameStr); - histClassesQA += Form("%s;", nameStr.Data()); - - // All reconstructed leptons matched to a 1 prong signal or MC 1 prong signal directly - std::vector mcnamesreco; - for (int isig = 0; isig < fMCSignals.size(); ++isig) { - TString nameStr2 = Form("TrackBarrel_%s_%s", cut.GetName(), fMCSignals.at(isig).GetName()); - mcnamesreco.push_back(nameStr2); - histClassesQA += Form("%s;", nameStr2.Data()); - } - fHistNamesMCMatchedQA.push_back(mcnamesreco); - } - - // Add histogram classes for each MC signal at generated level - std::vector mcnamesgen; - for (int isig = 0; isig < fMCSignals.size(); ++isig) { - TString nameStr2 = Form("MCTruthGen_%s", fMCSignals.at(isig).GetName()); - mcnamesgen.push_back(nameStr2); - histClassesQA += Form("%s;", nameStr2.Data()); - - } - fHistNamesMCQA.push_back(mcnamesgen); - - - fHistManQA = new HistogramManager("SingleElectronQA", "aa", VarManager::kNVars); - fHistManQA->SetUseDefaultVariableNames(kTRUE); - fHistManQA->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); - DefineHistograms(fHistManQA, histClassesQA.Data()); // define all histograms - VarManager::SetUseVars(fHistManQA->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill - fQASingleElectronList = fHistManQA->GetMainHistogramList(); - fMainList->Add(fQASingleElectronList); - - } - - fOutputList.setObject(fMainList); - -} - - Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; - Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; - - template - void runSelection(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) - { - - uint8_t eventFilter = 0; - std::map fMCEventLabels; - int fCounters = 0; //! [0] - particle counter, [1] - event counter - - - for (auto& event : events) { - - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); - VarManager::ResetValues(0, VarManager::kNMCParticleVariables); - // fill event information which might be needed in histograms that combine track and event properties - VarManager::FillEvent(event); - //if(!fEventCut->IsSelected(VarManager::fgValues)) continue; - eventFilter = uint32_t(event.isEventSelected()); - if(!eventFilter) continue; - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - VarManager::FillEvent(event.reducedMCevent()); - } - - // Look if we did not already saw the collision and fill the denominator of the single electron efficiency - Int_t globalindexmc = -1; - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - auto mcEvent = event.reducedMCevent(); - globalindexmc = mcEvent.globalIndex(); - } - if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { - fMCEventLabels[globalindexmc] = fCounters; - fCounters++; - // skimmed data - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); - groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenTrack(groupedMCTracks); - } - } - - - // Loop over reconstructed tracks belonging to the event and fill the numerator of the efficiency as well as the resolution map - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); - runRecTrack(groupedTracks,tracksMC); - } - - }// end loop over events - - } - - template - void runMCGenTrack(TTracksMC const& groupedMCTracks) - { - - - // loop over mc stack and fill histograms for pure MC truth signals - // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event - //auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); - for (auto& mctrack : groupedMCTracks) { - VarManager::FillTrack(mctrack); - // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. - // NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member. - // TODO: Use the mcReducedFlags to select signals - int isig = 0; - for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - if ((*sig).CheckSignal(true, groupedMCTracks, mctrack)) { - if(mctrack.pdgCode() > 0) { - dynamic_cast(fHistGenNegPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); - } - else { - dynamic_cast(fHistGenPosPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); - } - if (fConfigQA) fHistManQA->FillHistClass(Form("MCTruthGen_%s", (*sig).GetName()), VarManager::fgValues); - } - } - } - } - - template - void runRecTrack(TTracks const& groupedTracks, TTracksMC const& tracksMC) - { - - uint32_t filterMap = 0; - trackSel.reserve(groupedTracks.size()); - for (auto& track : groupedTracks) { - filterMap = 0; - - VarManager::FillTrack(track); // compute track quantities - - // compute MC matched quantities - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - VarManager::FillTrack(track.reducedMCTrack()); - } - if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { - if (!track.has_mcParticle()) { - continue; - } - auto mctrack = track.template mcParticle_as(); - VarManager::FillTrack(mctrack); - } - - // compute track selection and publish the bit map - int i = 0; - for (auto cut = fTrackCuts.begin(); cut != fTrackCuts.end(); cut++, i++) { - if ((*cut).IsSelected(VarManager::fgValues)) { - filterMap |= (uint32_t(1) << i); - if (fConfigQA) { - fHistManQA->FillHistClass(fHistNamesRecoQA[i].Data(), VarManager::fgValues); - } - } - } - trackSel(static_cast(filterMap)); - if (!filterMap) { - continue; - } - - // compute MC matching decisions - uint32_t mcDecision = 0; - int isig = 0; - for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - if ((*sig).CheckSignal(true, tracksMC, track.reducedMCTrack())) { - mcDecision |= (uint32_t(1) << isig); - } - } - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { - if (track.has_mcParticle()) { - auto mctrack = track.template mcParticle_as(); - if ((*sig).CheckSignal(true, tracksMC, mctrack)) { - mcDecision |= (uint32_t(1) << isig); - } - } - } - - } - - - // fill histograms - for (unsigned int i = 0; i < fMCSignals.size(); i++) { - if (!(mcDecision & (uint32_t(1) << i))) { - continue; - } - for (unsigned int j = 0; j < fTrackCuts.size(); j++) { - if (filterMap & (uint8_t(1) << j)) { - if(track.sign() < 0) { - dynamic_cast(fHistRecNegPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); - } - else { - dynamic_cast(fHistRecPosPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); - } - if(fConfigResolutionOn && (i==0) && (j==0)){ - - Double_t mcpt = -10000.; - Double_t mceta = -10000.; - Double_t mcphi = -1000.; - Int_t mcpdg = -10000.; - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto mctrack = track.reducedMCTrack(); - mcpt = mctrack.pt(); - mceta = mctrack.pt(); - mcphi = mctrack.pt(); - mcpdg = mctrack.pdgCode(); - } - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { - if (track.has_mcParticle()) { - auto mctrack = track.template mcParticle_as(); - mcpt = mctrack.pt(); - mceta = mctrack.pt(); - mcphi = mctrack.pt(); - mcpdg = mctrack.pdgCode(); - } - } - - - Double_t deltaptoverpt = -1000.; - if(mcpt>0.) deltaptoverpt = (mcpt-track.pt())/mcpt; - Double_t deltaeta = mceta-track.eta(); - Double_t deltaphi = mcphi-track.phi(); - dynamic_cast(fHistRes.at(0))->Fill(mcpt,deltaptoverpt); - dynamic_cast(fHistRes.at(1))->Fill(mcpt,deltaeta); - if(mcpdg<0){ - dynamic_cast(fHistRes.at(2))->Fill(mcpt,deltaphi); - }else { - dynamic_cast(fHistRes.at(3))->Fill(mcpt,deltaphi); - } - } - if(fConfigQA) fHistManQA->FillHistClass(fHistNamesMCMatchedQA[j][i].Data(), VarManager::fgValues); - } - } // end loop over cuts - } // end loop over MC signals - - - } // end loop over reconstructed track belonging to the events - - } - - - - void processSkimmed(soa::Filtered const& events, MyBarrelTracks const& tracks, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) - { - runSelection(events, tracks, eventsMC, tracksMC); - - - } - void processDummy(MyEvents&) - { - // do nothing - } - - PROCESS_SWITCH(AnalysisTrackSelection, processSkimmed, "Run barrel track selection on DQ skimmed tracks", false); - PROCESS_SWITCH(AnalysisTrackSelection, processDummy, "Dummy process function", false); -}; - - -struct AnalysisSameEventPairing { - - - //Produces dileptonList; - - // Filter based on previous components in the task - Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; - Filter filterBarrelTrackSelected = aod::emanalysisflags::isBarrelSelected > 0; - - Configurable fConfigTrackCuts{"cfgTrackCuts", "", "Comma separated list of barrel track cuts"}; - Configurable fConfigMCSignals{"cfgBarrelMCSignals", "", "Comma separated list of MC signals"}; - Configurable fConfigMinPt{"cfgMinPt", 0., "Fiducial min Pt for MC signal"}; - Configurable fConfigMaxPt{"cfgMaxPt", 10., "Fiducial max Pt for MC signal"}; - Configurable fConfigMinEta{"cfgMinEta", -0.8, "Fiducial min eta for MC signal"}; - Configurable fConfigMaxEta{"cfgMaxEta", 0.8, "Fiducial max eta for MC signal"}; - Configurable fConfigFlatTables{"cfgFlatTables", false, "Produce a single flat tables with all relevant information of the pairs and single tracks"}; - - // TODO: here we specify signals, however signal decisions are precomputed and stored in mcReducedFlags - // TODO: The tasks based on skimmed MC could/should rely ideally just on these flags - // TODO: special AnalysisCuts to be prepared in this direction - // TODO: cuts on the MC truth information to be added if needed - - - // 2D histos: mee and ptee - Configurable fConfigMinPtee{"cfgMinPtee", 0., "min Ptee in 2D histos"}; - Configurable fConfigMaxPtee{"cfgMaxPtee", 10., "max Ptee in 2D histos"}; - Configurable fConfigStepPtee{"cfgStepPtee", 100, "Nb of steps in ptee in 2D histos"}; - Configurable fConfigMinMee{"cfgMinMee", 0., "min Mee in 2D histos"}; - Configurable fConfigMaxMee{"cfgMaxMee", 3.5, "max Mee in 2D histos"}; - Configurable fConfigStepMee{"cfgStepMee", 600, "Nb of steps in Mee in 2D histos"}; - - // output lists - OutputObj fOutputList{"output"}; - THashList *fMainList; // Main list - THashList *fPairList; // 2D histos for MC and reconstructed signals - THashList *fQAPairList; // QA in case on with histo manager outputs - - // Cuts and signals - //AnalysisCompositeCut* fEventCut; // Taken from event selection part - std::vector fTrackCuts; // list of track cuts - std::vector fMCSignals; // list of signals with one prong to be checked: ULS 2D histos - - //2D histo vectors - std::vector fHistGenPair; - std::vector fHistGenSmearedPair; - std::vector fHistRecPair; - // Binning - std::vector fPteeBins; - std::vector fMeeBins; - - // QA: to be defined - - - void init(o2::framework::InitContext& context) - { - - printf("Start init\n"); - - // Create list output - fMainList = new THashList; - fMainList->SetOwner(kTRUE); - fMainList->SetName("pairselection"); - - // Create list output for 3D eta,phi,pt - fPairList = new THashList; - fPairList->SetOwner(kTRUE); - fPairList->SetName("Dielectron"); - - // Binning 2D histos - SetBinsLinear(fPteeBins, fConfigMinPtee, fConfigMaxPtee, fConfigStepPtee); - SetBinsLinear(fMeeBins, fConfigMinMee, fConfigMaxMee, fConfigStepMee); - const int fNpteeBins = fPteeBins.size()-1; - const int fNmeeBins = fMeeBins.size()-1; - - // Event cut: taken from event selection part - //fEventCut = new AnalysisCompositeCut(true); - //TString eventCutStr = fConfigEventCuts.value; - //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); - - // List of track cuts - TString cutNamesStr = fConfigTrackCuts.value; - if (!cutNamesStr.IsNull()) { - std::unique_ptr objArray(cutNamesStr.Tokenize(",")); - for (int icut = 0; icut < objArray->GetEntries(); ++icut) { - fTrackCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); - } - } - VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill - VarManager::SetDefaultVarNames(); - - // List of MC signals - TString configSigNamesStr = fConfigMCSignals.value; - std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); - for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { - MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); - if (sig) { - if (sig->GetNProngs() == 2) { // only 2 prong signals - fMCSignals.push_back(*sig); - } - // List of signal to be checked - } - } - - //Configure 2D histograms - // Create List with generated particles - TList* Generated = new TList(); - Generated->SetOwner(); - Generated->SetName("Generated"); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH2D* th2_tmp = new TH2D(Form("Ngen_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); - th2_tmp->Sumw2(); - fHistGenPair.push_back(th2_tmp); - Generated->Add(th2_tmp); - } - - // Create List with generated+smeared particles - TList* GeneratedSmeared = new TList(); - GeneratedSmeared->SetName("GeneratedSmeared"); - GeneratedSmeared->SetOwner(); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH2D* th2_tmp = new TH2D(Form("NgenSmeared_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); - th2_tmp->Sumw2(); - fHistGenSmearedPair.push_back(th2_tmp); - GeneratedSmeared->Add(th2_tmp); - } - - - fPairList->Add(Generated); - fPairList->Add(GeneratedSmeared); - - // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge - for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ - TList* list = new TList(); - list->SetName(fTrackCuts.at(list_i).GetName()); - list->SetOwner(); - - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH2D* th2_tmp = new TH2D(Form("Nrec_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); - th2_tmp->Sumw2(); - th2_tmp->SetDirectory(0x0); - fHistRecPair.push_back(th2_tmp); - list->Add(th2_tmp); - } - fPairList->Add(list); - } - fMainList->Add(fPairList); - - fOutputList.setObject(fMainList); - - } - - Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; - Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; - - template < uint32_t TEventFillMap, uint32_t TEventMCFillMap, uint32_t TTrackFillMap, typename TEvents, typename TTracks, typename TEventsMC, typename TTracksMC> - void runPairing(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) - { - - - std::map fMCEventLabels; - int fCounters = 0; //! [0] - particle counter, [1] - event counter - - - for (auto& event : events) { - - if (!event.isEventSelected()) { - return; - } - - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); - VarManager::ResetValues(0, VarManager::kNMCParticleVariables); - // fill event information which might be needed in histograms that combine track and event properties - VarManager::FillEvent(event); - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - VarManager::FillEvent(event.reducedMCevent()); - } - - // Look if we did not already saw the collision and fill the denominator of the single electron efficiency - Int_t globalindexmc = -1; - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - auto mcEvent = event.reducedMCevent(); - globalindexmc = mcEvent.globalIndex(); - } - if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { - fMCEventLabels[globalindexmc] = fCounters; - fCounters++; - // skimmed data - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); - groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenPair(groupedMCTracks); - } - } - - - //auto groupedTrackspos = posTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); - //groupedTrackspos.bindInternalIndicesTo(&tracks); - //auto groupedTracksneg = negTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); - //groupedTracksneg.bindInternalIndicesTo(&tracks); - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); - groupedTracks.bindInternalIndicesTo(&tracks); - runRecPair(groupedTracks,tracksMC); - } - - } // end loop over reconstructed event - } // end loop pairing function - - - template - void runMCGenPair(TTracksMC const& groupedMCTracks) - { - // - Double_t masse = 0.00051099895; // 0.5 MeV/c2 -> 0.0005 GeV/c2 - - for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { - - if((abs(t1.pdgCode()) != 11) || (abs(t2.pdgCode()) != 11)) continue; - if(t1.pdgCode()*t2.pdgCode() > 0) continue; // ULS only - - TLorentzVector Lvec1; - TLorentzVector Lvec2; - Lvec1.SetPtEtaPhiM(t1.pt(),t1.eta(),t1.phi(),masse); - Lvec2.SetPtEtaPhiM(t2.pt(),t2.eta(),t2.phi(),masse); - TLorentzVector LvecM = Lvec1 + Lvec2; - double mass = LvecM.M(); - double pairpt = LvecM.Pt(); - //double opangle = Lvec1.Angle(Lvec2.Vect()); - - //printf("Check before\n"); - // Fiducial cut - Bool_t genfidcut = kTRUE; - if((t1.eta()> fConfigMaxEta) || (t2.eta() > fConfigMaxEta) || (t1.eta() < fConfigMinEta) || (t2.eta() < fConfigMinEta) || (t1.pt()> fConfigMaxPt) || (t2.pt() > fConfigMaxPt) || (t1.pt() < fConfigMinPt) || (t2.pt() < fConfigMinPt)) genfidcut = kFALSE; - - int isig = 0; - for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - if ((*sig).CheckSignal(true, groupedMCTracks, t1, t2)) { - - // not smeared - if(genfidcut) dynamic_cast(fHistGenPair.at(isig))->Fill(mass,pairpt); - // need to implement smeared - - } - } - } //end of true pairing loop - } // end runMCGen - - - - template < uint32_t TTrackFillMap, typename TTracks, typename TTracksMC> - void runRecPair(TTracks const& tracks, TTracksMC const& tracksMC) - { - - // Loop over two track combinations - uint8_t twoTrackFilter = 0; - //uint32_t dileptonFilterMap = 0; - // uint32_t dileptonMcDecision = 0; - //dileptonList.reserve(1); - - - - for (auto& [t1, t2] : combinations(tracks, tracks)) { - - twoTrackFilter = uint32_t(t1.isBarrelSelected()) & uint32_t(t2.isBarrelSelected()); - - if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue - continue; - } - VarManager::FillPair(t1, t2); - - // run MC matching for this pair - uint32_t mcDecision = 0; - int isig = 0; - for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { // for skimmed DQ model - if ((*sig).CheckSignal(true, tracksMC, t1.reducedMCTrack(), t2.reducedMCTrack())) { - mcDecision |= (uint32_t(1) << isig); - } - } - if constexpr ((TTrackFillMap & VarManager::ObjTypes::Track) > 0) { // for Ao2d - if (t1.has_mcParticle() && t2.has_mcParticle()) { - auto mctrack1 = t1.template mcParticle_as(); - auto mctrack2 = t2.template mcParticle_as(); - if ((*sig).CheckSignal(true, tracksMC, mctrack1, mctrack2)) { - mcDecision |= (uint32_t(1) << isig); - } - } - } - - }// end of loop MC signals - - //dileptonFilterMap = twoTrackFilter; - //dileptonMcDecision = mcDecision; - //dileptonList(event, VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), dileptonFilterMap, dileptonMcDecision); - - - for (unsigned int i = 0; i < fMCSignals.size(); i++) { - if (!(mcDecision & (uint32_t(1) << i))) { - continue; - } - for (unsigned int j = 0; j < fTrackCuts.size(); j++) { - if (twoTrackFilter & (uint8_t(1) << j)) { - dynamic_cast(fHistRecPair.at(j * fMCSignals.size() + i))->Fill(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt]); - } - } - } - } - - } - - void processToEESkimmed(soa::Filtered const& events, - soa::Filtered const& tracks, - ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) - { - - runPairing< gkEventFillMap, gkMCEventFillMap, gkTrackFillMap>(events, tracks, eventsMC, tracksMC); - - } - - void processDummy(MyEvents&) - { - // do nothing - } - - PROCESS_SWITCH(AnalysisSameEventPairing, processToEESkimmed, "Run barrel barrel pairing on DQ skimmed tracks", false); - PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy process function", false); -}; - - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc)}; -} - -void DefineHistograms(HistogramManager* histMan, TString histClasses) -{ - // - // Define here the histograms for all the classes required in analysis. - // The histogram classes are provided in the histClasses string, separated by semicolon ";" - // The histogram classes and their components histograms are defined below depending on the name of the histogram class - // - std::unique_ptr objArray(histClasses.Tokenize(";")); - for (Int_t iclass = 0; iclass < objArray->GetEntries(); ++iclass) { - TString classStr = objArray->At(iclass)->GetName(); - histMan->AddHistClass(classStr.Data()); - - // NOTE: The level of detail for histogramming can be controlled via configurables - if (classStr.Contains("Event")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", "trigger,cent,mc"); - } - - if (classStr.Contains("Track")) { - if (classStr.Contains("Barrel")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", "its,tpcpid,dca,tofpid,mc"); - } - } - - if (classStr.Contains("Pairs")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel", "vertexing-barrel"); - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_dimuon", "vertexing-forward"); - } - - if (classStr.Contains("MCTruthGenPair")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair"); - histMan->AddHistogram(objArray->At(iclass)->GetName(), "Pt", "MC generator p_{T} distribution", false, 200, 0.0, 20.0, VarManager::kMCPt); - histMan->AddHistogram(objArray->At(iclass)->GetName(), "Eta", "MC generator #eta distribution", false, 500, -5.0, 5.0, VarManager::kMCEta); - histMan->AddHistogram(objArray->At(iclass)->GetName(), "Phi", "MC generator #varphi distribution", false, 500, -6.3, 6.3, VarManager::kMCPhi); - } - if (classStr.Contains("MCTruthGen")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth"); - } - if (classStr.Contains("DileptonsSelected")) { - dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel"); - } - - - } -}// end loop over histogram classes - -void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps){ - fBins.clear(); - const double stepSize = (max - min) / steps; - for (unsigned int i = 0; i < steps+1; ++i){ - fBins.push_back(i * stepSize + min); - } -} From 4b57647c1a944040d19934294208f872990ec330 Mon Sep 17 00:00:00 2001 From: Raphaelle Bailhache Date: Fri, 15 Jul 2022 01:43:58 +0200 Subject: [PATCH 3/4] Fixes --- PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx | 854 +++++++++++------------- 1 file changed, 392 insertions(+), 462 deletions(-) diff --git a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx index 5064efa6869..5bbd7cd3c45 100644 --- a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx +++ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx @@ -48,34 +48,31 @@ namespace o2::aod namespace emanalysisflags { - DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); - DECLARE_SOA_COLUMN(IsBarrelSelected, isBarrelSelected, int); -} // namespace reducedevent +DECLARE_SOA_COLUMN(IsEventSelected, isEventSelected, int); +DECLARE_SOA_COLUMN(IsBarrelSelected, isBarrelSelected, int); +} // namespace emanalysisflags - DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", emanalysisflags::IsEventSelected); - DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "BARRELTRACKCUTS", emanalysisflags::IsBarrelSelected); +DECLARE_SOA_TABLE(EventCuts, "AOD", "EVENTCUTS", emanalysisflags::IsEventSelected); +DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "BARRELTRACKCUTS", emanalysisflags::IsBarrelSelected); } // namespace o2::aod // using MyEvents = soa::Join; using MyEventsSelected = soa::Join; using MyBarrelTracks = soa::Join; -//using MyBarrelTracksWithCov = soa::Join; +// using MyBarrelTracksWithCov = soa::Join; using MyBarrelTracksSelected = soa::Join; -//using MyBarrelTracksSelectedWithCov = soa::Join; +// using MyBarrelTracksSelectedWithCov = soa::Join; - -// +// constexpr static uint32_t gkEventFillMap = VarManager::ObjTypes::ReducedEvent | VarManager::ObjTypes::ReducedEventExtended; constexpr static uint32_t gkMCEventFillMap = VarManager::ObjTypes::ReducedEventMC; constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; -//constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; +// constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkParticleMCFillMap = VarManager::ObjTypes::ParticleMC; - - void DefineHistograms(HistogramManager* histMan, TString histClasses); -void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps); +void SetBinsLinear(std::vector& fBins, const double min, const double max, const unsigned int steps); struct AnalysisEventSelection { @@ -144,121 +141,104 @@ struct AnalysisEventQa { Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; OutputObj fOutputList{"output"}; - THashList *fMainList; - TH1D *fNbRecCollisionPerMCCollision; - TProfile *fNbmctrack; - TH1D *fNbMcEvent; - TH1D *fNbRecEvent; + THashList* fMainList; + TH1D* fNbRecCollisionPerMCCollision; + TProfile* fNbmctrack; + TH1D* fNbMcEvent; + TH1D* fNbRecEvent; void init(o2::framework::InitContext&) { - - VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill + + VarManager::SetUseVars(AnalysisCut::fgUsedVars); VarManager::SetDefaultVarNames(); fMainList = new THashList; fMainList->SetOwner(kTRUE); fMainList->SetName("OAnalysisEventQA"); - fNbRecCollisionPerMCCollision = new TH1D("MCEvRecEv","",20,0.,20.); - fNbmctrack = new TProfile("mctrack","",20,0.,20.); - fNbMcEvent = new TH1D("MCEvent","",1,0.,1.); - fNbRecEvent = new TH1D("RecEvent","",1,0.,1.); + fNbRecCollisionPerMCCollision = new TH1D("MCEvRecEv", "", 20, 0., 20.); + fNbmctrack = new TProfile("mctrack", "", 20, 0., 20.); + fNbMcEvent = new TH1D("MCEvent", "", 1, 0., 1.); + fNbRecEvent = new TH1D("RecEvent", "", 1, 0., 1.); fMainList->Add(fNbRecCollisionPerMCCollision); fMainList->Add(fNbmctrack); fMainList->Add(fNbMcEvent); fMainList->Add(fNbRecEvent); fOutputList.setObject(fMainList); - } Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; - + template void runSelection(TEvents const& events, TEventsMC const& eventsMC, TTracksMC const& tracksMC) { - + uint8_t eventFilter = 0; std::map fMCEventNbmctrack; std::map fMCEventNbReco; std::map fMCEventLabels; - - - int fMCCounters = 0; //! - mc event counter - int fEvCounters = 0; // ! - rec event counter + int fMCCounters = 0; + int fEvCounters = 0; // First loop for (auto& event : events) { - - - // Reset the values array - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); - VarManager::FillEvent(event); - eventFilter = uint32_t(event.isEventSelected()); - if(!eventFilter) continue; - fEvCounters++; - fNbRecEvent->Fill(0.5); - - - Int_t midrap = 0; - Int_t globalindexmc = -1; + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event); + eventFilter = uint32_t(event.isEventSelected()); + if (!eventFilter) + continue; + fEvCounters++; + fNbRecEvent->Fill(0.5); - // skimmed data - if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - - auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); - midrap = dNdetach(groupedMCTracks); + Int_t midrap = 0; + Int_t globalindexmc = -1; - auto mcEvent = event.reducedMCevent(); - globalindexmc = mcEvent.globalIndex(); + // skimmed data + if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + midrap = dNdetach(groupedMCTracks); - } + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); + } - if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { - fMCEventLabels[globalindexmc] = fMCCounters; - fMCEventNbReco[globalindexmc] = 1; - fMCEventNbmctrack[globalindexmc] = midrap; - fNbMcEvent->Fill(0.5); - fMCCounters++; - - } - else { - fMCEventNbReco[globalindexmc] = fMCEventNbReco.find(globalindexmc)->second + 1; - } + if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { + fMCEventLabels[globalindexmc] = fMCCounters; + fMCEventNbReco[globalindexmc] = 1; + fMCEventNbmctrack[globalindexmc] = midrap; + fNbMcEvent->Fill(0.5); + fMCCounters++; + } else { + fMCEventNbReco[globalindexmc] = fMCEventNbReco.find(globalindexmc)->second + 1; + } - } // end loop over events - + for (const auto& [mcEv, NbRecEv] : fMCEventNbReco) { - fNbRecCollisionPerMCCollision->Fill(fMCEventNbReco.find(mcEv)->second); - } for (const auto& [mcEv, NbRecEv] : fMCEventNbmctrack) { - - fNbmctrack->Fill(fMCEventNbReco.find(mcEv)->second,fMCEventNbmctrack.find(mcEv)->second); - + fNbmctrack->Fill(fMCEventNbReco.find(mcEv)->second, fMCEventNbmctrack.find(mcEv)->second); } - } template - Int_t dNdetach(TTracksMC const& groupedMCTracks){ - + Int_t dNdetach(TTracksMC const& groupedMCTracks) + { + Int_t midrap = 0; - for(auto mctrack : groupedMCTracks){ - if(TMath::Abs(mctrack.eta())<0.5 && mctrack.isPhysicalPrimary() && (TMath::Abs(mctrack.pdgCode())==211 || mctrack.pdgCode()==111)) { - midrap++; + for (auto mctrack : groupedMCTracks) { + if (TMath::Abs(mctrack.eta()) < 0.5 && mctrack.isPhysicalPrimary() && (TMath::Abs(mctrack.pdgCode()) == 211 || mctrack.pdgCode() == 111)) { + midrap++; } } return midrap; - } - - void processSkimmed(soa::Filtered const&events, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + void processSkimmed(soa::Filtered const& events, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) { runSelection(events, eventsMC, tracksMC); } @@ -269,19 +249,18 @@ struct AnalysisEventQa { PROCESS_SWITCH(AnalysisEventQa, processSkimmed, "Run event selection on DQ skimmed events", false); PROCESS_SWITCH(AnalysisEventQa, processDummy, "Dummy process function", false); - }; struct AnalysisTrackSelection { - + Produces trackSel; Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; - + // configurables - //Configurable fConfigEventCuts{"cfgEventCuts", "eventStandardNoINT7", "Event selection"}; + // Configurable fConfigEventCuts{"cfgEventCuts", "eventStandardNoINT7", "Event selection"}; Configurable fConfigCuts{"cfgTrackCuts", "jpsiPID1", "Comma separated list of barrel track cuts"}; Configurable fConfigMCSignals{"cfgTrackMCSignals", "", "Comma separated list of MC signals"}; - + // 3D histos Configurable fConfigUsePtVec{"cfgUsePtVec", true, "If true, non-linear pt bins predefined"}; Configurable fConfigMinPt{"cfgMinPt", 0., "min Pt in 3D histos"}; @@ -308,31 +287,30 @@ struct AnalysisTrackSelection { Configurable fConfigMaxDeltaPhi{"cfgMaxDeltaPhi", 0.5, "max delta Phi in res histos"}; Configurable fConfigStepDeltaPhi{"cfgStepDeltaPhi", 500, "Nb of steps in delta Phi in res histos"}; - Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; // output lists OutputObj fOutputList{"output"}; - THashList *fMainList; // Main list - THashList *fSingleElectronList; // 3D histos for MC and reconstructed signals - THashList *fResolutionList; // Resolution histograms - THashList *fQASingleElectronList; // QA in case on with histo manager outputs + THashList* fMainList; // Main list + THashList* fSingleElectronList; // 3D histos for MC and reconstructed signals + THashList* fResolutionList; // Resolution histograms + THashList* fQASingleElectronList; // QA in case on with histo manager outputs // Cuts and signals - //AnalysisCompositeCut* fEventCut; // Taken from event selection part + // AnalysisCompositeCut* fEventCut; // Taken from event selection part std::vector fTrackCuts; // list of track cuts - AnalysisCompositeCut* fTrackCutsRes; // track cut for resolution map - std::vector fMCSignals; // list of signals to be checked - MCSignal* fMCSignalRes; // signal for res + AnalysisCompositeCut* fTrackCutsRes; // track cut for resolution map + std::vector fMCSignals; // list of signals to be checked + MCSignal* fMCSignalRes; // signal for res - //3D histos + // 3D histos std::vector fHistGenPosPart; std::vector fHistGenNegPart; std::vector fHistGenSmearedPosPart; std::vector fHistGenSmearedNegPart; std::vector fHistRecPosPart; std::vector fHistRecNegPart; - // Binning + // Binning std::vector fPtBins; std::vector fEtaBins; std::vector fPhiBins; @@ -343,47 +321,46 @@ struct AnalysisTrackSelection { std::vector fPtResBins; std::vector fDeltaEtaBins; std::vector fDeltaPhiBins; - + // QA - HistogramManager* fHistManQA; // histo manager - std::vector fHistNamesRecoQA; // list of histo names for all reconstructed tracks in histo manager + HistogramManager* fHistManQA; // histo manager + std::vector fHistNamesRecoQA; // list of histo names for all reconstructed tracks in histo manager std::vector> fHistNamesMCMatchedQA; // list of histo names for reconstructed signals in histo manager - std::vector> fHistNamesMCQA; // list of histo names for generated signals in histo manager - + std::vector> fHistNamesMCQA; // list of histo names for generated signals in histo manager + void init(o2::framework::InitContext&) { - + // Create list output fMainList = new THashList; fMainList->SetOwner(kTRUE); fMainList->SetName("trackselection"); - + // Create list output for 3D eta,phi,pt fSingleElectronList = new THashList; fSingleElectronList->SetOwner(kTRUE); fSingleElectronList->SetName("SingleElectron"); // Binning 3D histos - if(fConfigUsePtVec){ + if (fConfigUsePtVec) { const Int_t Npt = 68; - Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.}; - std::vector v_pte(pte,std::end(pte)); + Double_t pte[Npt] = {0.00, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.155, 0.16, 0.165, 0.17, 0.175, 0.18, 0.185, 0.19, 0.195, 0.20, 0.205, 0.21, 0.215, 0.22, 0.225, 0.23, 0.235, 0.24, 0.245, 0.25, 0.255, 0.26, 0.265, 0.27, 0.275, 0.28, 0.285, 0.29, 0.295, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.43, 0.46, 0.49, 0.52, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.90, 1.00, 1.10, 1.20, 1.40, 1.60, 1.80, 2.00, 2.40, 2.80, 3.20, 3.70, 4.50, 6.00, 8.00, 10.}; + std::vector v_pte(pte, std::end(pte)); fPtBins = v_pte; - } - else { + } else { SetBinsLinear(fPtBins, fConfigMinPt, fConfigMaxPt, fConfigStepPt); } SetBinsLinear(fEtaBins, fConfigMinEta, fConfigMaxEta, fConfigStepEta); SetBinsLinear(fPhiBins, fConfigMinPhi, fConfigMaxPhi, fConfigStepPhi); - const int fNptBins = fPtBins.size()-1; - const int fNetaBins = fEtaBins.size()-1; - const int fNphiBins = fPhiBins.size()-1; + const int fNptBins = fPtBins.size() - 1; + const int fNetaBins = fEtaBins.size() - 1; + const int fNphiBins = fPhiBins.size() - 1; // Event cut: taken from event selection part - //fEventCut = new AnalysisCompositeCut(true); - //TString eventCutStr = fConfigEventCuts.value; - //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); - + // fEventCut = new AnalysisCompositeCut(true); + // TString eventCutStr = fConfigEventCuts.value; + // fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + // List of track cuts TString cutNamesStr = fConfigCuts.value; if (!cutNamesStr.IsNull()) { @@ -394,32 +371,32 @@ struct AnalysisTrackSelection { } VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill VarManager::SetDefaultVarNames(); - + // List of MC signals TString configSigNamesStr = fConfigMCSignals.value; std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); if (sig) { - if (sig->GetNProngs() != 1) { // NOTE: only 1 prong signals - continue; - } - // List of signal to be checked - fMCSignals.push_back(*sig); + if (sig->GetNProngs() != 1) { // NOTE: only 1 prong signals + continue; + } + // List of signal to be checked + fMCSignals.push_back(*sig); } } - //Configure 3D histograms + // Configure 3D histograms // Create List with generated particles TList* Generated = new TList(); Generated->SetOwner(); Generated->SetName("Generated"); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + for (unsigned int i = 0; i < fMCSignals.size(); ++i) { + TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()), ";p_{T};#eta;#varphi", fNptBins, fPtBins.data(), fNetaBins, fEtaBins.data(), fNphiBins, fPhiBins.data()); th3_tmp_pos->Sumw2(); fHistGenPosPart.push_back(th3_tmp_pos); Generated->Add(th3_tmp_pos); - TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()), ";p_{T};#eta;#varphi", fNptBins, fPtBins.data(), fNetaBins, fEtaBins.data(), fNphiBins, fPhiBins.data()); th3_tmp_neg->Sumw2(); fHistGenNegPart.push_back(th3_tmp_neg); Generated->Add(th3_tmp_neg); @@ -428,151 +405,140 @@ struct AnalysisTrackSelection { TList* GeneratedSmeared = new TList(); GeneratedSmeared->SetName("GeneratedSmeared"); GeneratedSmeared->SetOwner(); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + for (unsigned int i = 0; i < fMCSignals.size(); ++i) { + TH3D* th3_tmp_pos = new TH3D(Form("Ngen_Pos_%s", fMCSignals.at(i).GetName()), ";p_{T};#eta;#varphi", fNptBins, fPtBins.data(), fNetaBins, fEtaBins.data(), fNphiBins, fPhiBins.data()); th3_tmp_pos->Sumw2(); fHistGenSmearedPosPart.push_back(th3_tmp_pos); GeneratedSmeared->Add(th3_tmp_pos); - TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); + TH3D* th3_tmp_neg = new TH3D(Form("Ngen_Neg_%s", fMCSignals.at(i).GetName()), ";p_{T};#eta;#varphi", fNptBins, fPtBins.data(), fNetaBins, fEtaBins.data(), fNphiBins, fPhiBins.data()); th3_tmp_neg->Sumw2(); fHistGenSmearedNegPart.push_back(th3_tmp_neg); GeneratedSmeared->Add(th3_tmp_neg); } - + fSingleElectronList->Add(Generated); fSingleElectronList->Add(GeneratedSmeared); // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge - for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ + for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i) { TList* list = new TList(); list->SetName(fTrackCuts.at(list_i).GetName()); list->SetOwner(); - - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH3D* th3_tmp_pos = new TH3D(Form("Nrec_Pos_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_pos->Sumw2(); - th3_tmp_pos->SetDirectory(0x0); - fHistRecPosPart.push_back(th3_tmp_pos); - list->Add(th3_tmp_pos); - TH3D* th3_tmp_neg = new TH3D(Form("Nrec_Neg_%s", fMCSignals.at(i).GetName()),";p_{T};#eta;#varphi",fNptBins,fPtBins.data(),fNetaBins,fEtaBins.data(),fNphiBins,fPhiBins.data()); - th3_tmp_neg->Sumw2(); - th3_tmp_neg->SetDirectory(0x0); - fHistRecNegPart.push_back(th3_tmp_neg); - list->Add(th3_tmp_neg); - + + for (unsigned int i = 0; i < fMCSignals.size(); ++i) { + TH3D* th3_tmp_pos = new TH3D(Form("Nrec_Pos_%s", fMCSignals.at(i).GetName()), ";p_{T};#eta;#varphi", fNptBins, fPtBins.data(), fNetaBins, fEtaBins.data(), fNphiBins, fPhiBins.data()); + th3_tmp_pos->Sumw2(); + th3_tmp_pos->SetDirectory(0x0); + fHistRecPosPart.push_back(th3_tmp_pos); + list->Add(th3_tmp_pos); + TH3D* th3_tmp_neg = new TH3D(Form("Nrec_Neg_%s", fMCSignals.at(i).GetName()), ";p_{T};#eta;#varphi", fNptBins, fPtBins.data(), fNetaBins, fEtaBins.data(), fNphiBins, fPhiBins.data()); + th3_tmp_neg->Sumw2(); + th3_tmp_neg->SetDirectory(0x0); + fHistRecNegPart.push_back(th3_tmp_neg); + list->Add(th3_tmp_neg); } fSingleElectronList->Add(list); } fMainList->Add(fSingleElectronList); // Resolution histogramms - if(fConfigResolutionOn) { + if (fConfigResolutionOn) { // Binning 3D histos - if(fConfigUsePtVecRes){ - const Int_t Npt = 73; - Double_t pte[Npt] = {0.00,0.10,0.11,0.12,0.13,0.14,0.15,0.155,0.16,0.165,0.17,0.175,0.18,0.185,0.19,0.195,0.20,0.205,0.21,0.215,0.22,0.225,0.23,0.235,0.24,0.245,0.25,0.255,0.26,0.265,0.27,0.275,0.28,0.285,0.29,0.295,0.30,0.32,0.34,0.36,0.38,0.40,0.43,0.46,0.49,0.52,0.55,0.60,0.65,0.70,0.75,0.80,0.90,1.00,1.10,1.20,1.40,1.60,1.80,2.00,2.40,2.80,3.20,3.70,4.50,6.00,8.00,10.,12.0,14.,16.,18.,20.}; - std::vector v_pte(pte,std::end(pte)); - fPtResBins = v_pte; - } - else { - SetBinsLinear(fPtResBins, fConfigMinPtRes, fConfigMaxPtRes, fConfigStepPtRes); + if (fConfigUsePtVecRes) { + const Int_t Npt = 73; + Double_t pte[Npt] = {0.00, 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.155, 0.16, 0.165, 0.17, 0.175, 0.18, 0.185, 0.19, 0.195, 0.20, 0.205, 0.21, 0.215, 0.22, 0.225, 0.23, 0.235, 0.24, 0.245, 0.25, 0.255, 0.26, 0.265, 0.27, 0.275, 0.28, 0.285, 0.29, 0.295, 0.30, 0.32, 0.34, 0.36, 0.38, 0.40, 0.43, 0.46, 0.49, 0.52, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80, 0.90, 1.00, 1.10, 1.20, 1.40, 1.60, 1.80, 2.00, 2.40, 2.80, 3.20, 3.70, 4.50, 6.00, 8.00, 10., 12.0, 14., 16., 18., 20.}; + std::vector v_pte(pte, std::end(pte)); + fPtResBins = v_pte; + } else { + SetBinsLinear(fPtResBins, fConfigMinPtRes, fConfigMaxPtRes, fConfigStepPtRes); } SetBinsLinear(fDeltaEtaBins, fConfigMinDeltaEta, fConfigMaxDeltaEta, fConfigStepDeltaEta); SetBinsLinear(fDeltaPhiBins, fConfigMinDeltaPhi, fConfigMaxDeltaPhi, fConfigStepDeltaPhi); - const int fNptresBins = fPtResBins.size()-1; - const int fNDeltaetaBins = fDeltaEtaBins.size()-1; - const int fNDeltaphiBins = fDeltaPhiBins.size()-1; - + const int fNptresBins = fPtResBins.size() - 1; + const int fNDeltaetaBins = fDeltaEtaBins.size() - 1; + const int fNDeltaphiBins = fDeltaPhiBins.size() - 1; fResolutionList = new THashList; fResolutionList->SetOwner(kTRUE); fResolutionList->SetName("Resolution"); printf("Histos\n"); - TH2D *thPtGen_DeltaPtOverPtGen = new TH2D("PtGen_DeltaPtOverPtGen", "", fNptresBins, fPtResBins.data(), fConfigStepDeltaPt, -1., +1.); - TH2D *thPtGen_DeltaEta = new TH2D("PtGen_DeltaEta", "", fNptresBins, fPtResBins.data(), fNDeltaetaBins, fDeltaEtaBins.data()); - TH2D *thPtGen_DeltaPhi_Ele = new TH2D("PtGen_DeltaPhi_Ele", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); - TH2D *thPtGen_DeltaPhi_Pos = new TH2D("PtGen_DeltaPhi_Pos", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); - - thPtGen_DeltaPtOverPtGen ->Sumw2(); - thPtGen_DeltaEta ->Sumw2(); - thPtGen_DeltaPhi_Ele ->Sumw2(); - thPtGen_DeltaPhi_Pos ->Sumw2(); - - thPtGen_DeltaPtOverPtGen ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaPtOverPtGen ->GetYaxis()->SetTitle("(p^{gen}_{T} - p^{rec}_{T}) / p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaEta ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaEta ->GetYaxis()->SetTitle("#eta^{gen} - #eta^{rec}"); - thPtGen_DeltaPhi_Ele ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaPhi_Ele ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); - thPtGen_DeltaPhi_Pos ->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); - thPtGen_DeltaPhi_Pos ->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); - - printf("Add\n"); + TH2D* thPtGen_DeltaPtOverPtGen = new TH2D("PtGen_DeltaPtOverPtGen", "", fNptresBins, fPtResBins.data(), fConfigStepDeltaPt, -1., +1.); + TH2D* thPtGen_DeltaEta = new TH2D("PtGen_DeltaEta", "", fNptresBins, fPtResBins.data(), fNDeltaetaBins, fDeltaEtaBins.data()); + TH2D* thPtGen_DeltaPhi_Ele = new TH2D("PtGen_DeltaPhi_Ele", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); + TH2D* thPtGen_DeltaPhi_Pos = new TH2D("PtGen_DeltaPhi_Pos", "", fNptresBins, fPtResBins.data(), fNDeltaphiBins, fDeltaPhiBins.data()); + + thPtGen_DeltaPtOverPtGen->Sumw2(); + thPtGen_DeltaEta->Sumw2(); + thPtGen_DeltaPhi_Ele->Sumw2(); + thPtGen_DeltaPhi_Pos->Sumw2(); + + thPtGen_DeltaPtOverPtGen->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPtOverPtGen->GetYaxis()->SetTitle("(p^{gen}_{T} - p^{rec}_{T}) / p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaEta->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaEta->GetYaxis()->SetTitle("#eta^{gen} - #eta^{rec}"); + thPtGen_DeltaPhi_Ele->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPhi_Ele->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); + thPtGen_DeltaPhi_Pos->GetXaxis()->SetTitle("p^{gen}_{T} (GeV/c)"); + thPtGen_DeltaPhi_Pos->GetYaxis()->SetTitle("#varphi^{gen} - #varphi^{rec} (rad)"); + fHistRes.push_back(thPtGen_DeltaPtOverPtGen); fHistRes.push_back(thPtGen_DeltaEta); fHistRes.push_back(thPtGen_DeltaPhi_Ele); fHistRes.push_back(thPtGen_DeltaPhi_Pos); - + fResolutionList->Add(thPtGen_DeltaPtOverPtGen); fResolutionList->Add(thPtGen_DeltaEta); fResolutionList->Add(thPtGen_DeltaPhi_Ele); fResolutionList->Add(thPtGen_DeltaPhi_Pos); - - fMainList->Add(fResolutionList); + fMainList->Add(fResolutionList); } - // Configure QA histogram classes - if (fConfigQA){ + if (fConfigQA) { // Create list output for QA fQASingleElectronList = new THashList; fQASingleElectronList->SetOwner(kTRUE); fQASingleElectronList->SetName("SingleElectronQA"); TString histClassesQA = ""; for (auto& cut : fTrackCuts) { - - // All reconstructed leptons - TString nameStr = Form("TrackBarrel_%s", cut.GetName()); - fHistNamesRecoQA.push_back(nameStr); - histClassesQA += Form("%s;", nameStr.Data()); - - // All reconstructed leptons matched to a 1 prong signal or MC 1 prong signal directly - std::vector mcnamesreco; - for (int isig = 0; isig < fMCSignals.size(); ++isig) { + + // All reconstructed leptons + TString nameStr = Form("TrackBarrel_%s", cut.GetName()); + fHistNamesRecoQA.push_back(nameStr); + histClassesQA += Form("%s;", nameStr.Data()); + + // All reconstructed leptons matched to a 1 prong signal or MC 1 prong signal directly + std::vector mcnamesreco; + for (int isig = 0; isig < fMCSignals.size(); ++isig) { TString nameStr2 = Form("TrackBarrel_%s_%s", cut.GetName(), fMCSignals.at(isig).GetName()); mcnamesreco.push_back(nameStr2); - histClassesQA += Form("%s;", nameStr2.Data()); - } - fHistNamesMCMatchedQA.push_back(mcnamesreco); + histClassesQA += Form("%s;", nameStr2.Data()); + } + fHistNamesMCMatchedQA.push_back(mcnamesreco); } - + // Add histogram classes for each MC signal at generated level std::vector mcnamesgen; for (int isig = 0; isig < fMCSignals.size(); ++isig) { - TString nameStr2 = Form("MCTruthGen_%s", fMCSignals.at(isig).GetName()); - mcnamesgen.push_back(nameStr2); - histClassesQA += Form("%s;", nameStr2.Data()); - + TString nameStr2 = Form("MCTruthGen_%s", fMCSignals.at(isig).GetName()); + mcnamesgen.push_back(nameStr2); + histClassesQA += Form("%s;", nameStr2.Data()); } fHistNamesMCQA.push_back(mcnamesgen); - - + fHistManQA = new HistogramManager("SingleElectronQA", "aa", VarManager::kNVars); fHistManQA->SetUseDefaultVariableNames(kTRUE); fHistManQA->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); - DefineHistograms(fHistManQA, histClassesQA.Data()); // define all histograms - VarManager::SetUseVars(fHistManQA->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + DefineHistograms(fHistManQA, histClassesQA.Data()); // define all histograms + VarManager::SetUseVars(fHistManQA->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill fQASingleElectronList = fHistManQA->GetMainHistogramList(); fMainList->Add(fQASingleElectronList); - } - fOutputList.setObject(fMainList); - -} + } Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; @@ -584,57 +550,53 @@ struct AnalysisTrackSelection { uint8_t eventFilter = 0; std::map fMCEventLabels; int fCounters = 0; //! [0] - particle counter, [1] - event counter - - + for (auto& event : events) { - - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); VarManager::ResetValues(0, VarManager::kNMCParticleVariables); // fill event information which might be needed in histograms that combine track and event properties VarManager::FillEvent(event); - //if(!fEventCut->IsSelected(VarManager::fgValues)) continue; + // if(!fEventCut->IsSelected(VarManager::fgValues)) continue; eventFilter = uint32_t(event.isEventSelected()); - if(!eventFilter) continue; + if (!eventFilter) + continue; if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - VarManager::FillEvent(event.reducedMCevent()); + VarManager::FillEvent(event.reducedMCevent()); } // Look if we did not already saw the collision and fill the denominator of the single electron efficiency Int_t globalindexmc = -1; if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - auto mcEvent = event.reducedMCevent(); - globalindexmc = mcEvent.globalIndex(); + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); } if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { fMCEventLabels[globalindexmc] = fCounters; fCounters++; - // skimmed data - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); - groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenTrack(groupedMCTracks); - } + // skimmed data + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGenTrack(groupedMCTracks); + } } - - + // Loop over reconstructed tracks belonging to the event and fill the numerator of the efficiency as well as the resolution map if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); - runRecTrack(groupedTracks,tracksMC); + auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); + runRecTrack(groupedTracks, tracksMC); } - }// end loop over events - + } // end loop over events } - + template void runMCGenTrack(TTracksMC const& groupedMCTracks) { - // loop over mc stack and fill histograms for pure MC truth signals // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event - //auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); + // auto groupedMCTracks = tracksMC.sliceBy(aod::reducedtrackMC::reducedMCeventId, event.reducedMCevent().globalIndex()); for (auto& mctrack : groupedMCTracks) { VarManager::FillTrack(mctrack); // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. @@ -642,14 +604,14 @@ struct AnalysisTrackSelection { // TODO: Use the mcReducedFlags to select signals int isig = 0; for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - if ((*sig).CheckSignal(true, groupedMCTracks, mctrack)) { - if(mctrack.pdgCode() > 0) { - dynamic_cast(fHistGenNegPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); - } - else { - dynamic_cast(fHistGenPosPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); - } - if (fConfigQA) fHistManQA->FillHistClass(Form("MCTruthGen_%s", (*sig).GetName()), VarManager::fgValues); + if ((*sig).CheckSignal(true, groupedMCTracks, mctrack)) { + if (mctrack.pdgCode() > 0) { + dynamic_cast(fHistGenNegPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); + } else { + dynamic_cast(fHistGenPosPart.at(isig))->Fill(mctrack.pt(), mctrack.eta(), mctrack.phi()); + } + if (fConfigQA) + fHistManQA->FillHistClass(Form("MCTruthGen_%s", (*sig).GetName()), VarManager::fgValues); } } } @@ -663,101 +625,90 @@ struct AnalysisTrackSelection { trackSel.reserve(groupedTracks.size()); for (auto& track : groupedTracks) { filterMap = 0; - + VarManager::FillTrack(track); // compute track quantities // compute MC matched quantities if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - VarManager::FillTrack(track.reducedMCTrack()); + VarManager::FillTrack(track.reducedMCTrack()); } - - + // compute track selection and publish the bit map int i = 0; for (auto cut = fTrackCuts.begin(); cut != fTrackCuts.end(); cut++, i++) { - if ((*cut).IsSelected(VarManager::fgValues)) { - filterMap |= (uint32_t(1) << i); - if (fConfigQA) { - fHistManQA->FillHistClass(fHistNamesRecoQA[i].Data(), VarManager::fgValues); - } - } + if ((*cut).IsSelected(VarManager::fgValues)) { + filterMap |= (uint32_t(1) << i); + if (fConfigQA) { + fHistManQA->FillHistClass(fHistNamesRecoQA[i].Data(), VarManager::fgValues); + } + } } trackSel(static_cast(filterMap)); if (!filterMap) { - continue; + continue; } - + // compute MC matching decisions uint32_t mcDecision = 0; int isig = 0; for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - if ((*sig).CheckSignal(true, tracksMC, track.reducedMCTrack())) { - mcDecision |= (uint32_t(1) << isig); - } - } - + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + if ((*sig).CheckSignal(true, tracksMC, track.reducedMCTrack())) { + mcDecision |= (uint32_t(1) << isig); + } + } } - // fill histograms for (unsigned int i = 0; i < fMCSignals.size(); i++) { - if (!(mcDecision & (uint32_t(1) << i))) { - continue; - } - for (unsigned int j = 0; j < fTrackCuts.size(); j++) { - if (filterMap & (uint8_t(1) << j)) { - if(track.sign() < 0) { - dynamic_cast(fHistRecNegPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); - } - else { - dynamic_cast(fHistRecPosPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); - } - if(fConfigResolutionOn && (i==0) && (j==0)){ - - Double_t mcpt = -10000.; - Double_t mceta = -10000.; - Double_t mcphi = -1000.; - Int_t mcpdg = -10000.; - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto mctrack = track.reducedMCTrack(); - mcpt = mctrack.pt(); - mceta = mctrack.eta(); - mcphi = mctrack.phi(); - mcpdg = mctrack.pdgCode(); - } - - Double_t deltaptoverpt = -1000.; - if(mcpt>0.) deltaptoverpt = (mcpt-track.pt())/mcpt; - Double_t deltaeta = mceta-track.eta(); - Double_t deltaphi = mcphi-track.phi(); - dynamic_cast(fHistRes.at(0))->Fill(mcpt,deltaptoverpt); - dynamic_cast(fHistRes.at(1))->Fill(mcpt,deltaeta); - if(mcpdg<0){ - dynamic_cast(fHistRes.at(2))->Fill(mcpt,deltaphi); - }else { - dynamic_cast(fHistRes.at(3))->Fill(mcpt,deltaphi); - } - } - if(fConfigQA) fHistManQA->FillHistClass(fHistNamesMCMatchedQA[j][i].Data(), VarManager::fgValues); - } - } // end loop over cuts + if (!(mcDecision & (uint32_t(1) << i))) { + continue; + } + for (unsigned int j = 0; j < fTrackCuts.size(); j++) { + if (filterMap & (uint8_t(1) << j)) { + if (track.sign() < 0) { + dynamic_cast(fHistRecNegPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); + } else { + dynamic_cast(fHistRecPosPart.at(j * fMCSignals.size() + i))->Fill(track.pt(), track.eta(), track.phi()); + } + if (fConfigResolutionOn && (i == 0) && (j == 0)) { + + Double_t mcpt = -10000.; + Double_t mceta = -10000.; + Double_t mcphi = -1000.; + Int_t mcpdg = -10000.; + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto mctrack = track.reducedMCTrack(); + mcpt = mctrack.pt(); + mceta = mctrack.eta(); + mcphi = mctrack.phi(); + mcpdg = mctrack.pdgCode(); + } + Double_t deltaptoverpt = -1000.; + if (mcpt > 0.) + deltaptoverpt = (mcpt - track.pt()) / mcpt; + Double_t deltaeta = mceta - track.eta(); + Double_t deltaphi = mcphi - track.phi(); + dynamic_cast(fHistRes.at(0))->Fill(mcpt, deltaptoverpt); + dynamic_cast(fHistRes.at(1))->Fill(mcpt, deltaeta); + if (mcpdg < 0) { + dynamic_cast(fHistRes.at(2))->Fill(mcpt, deltaphi); + } else { + dynamic_cast(fHistRes.at(3))->Fill(mcpt, deltaphi); + } + } + if (fConfigQA) + fHistManQA->FillHistClass(fHistNamesMCMatchedQA[j][i].Data(), VarManager::fgValues); + } + } // end loop over cuts } // end loop over MC signals - - } // end loop over reconstructed track belonging to the events - } - - - void processSkimmed(soa::Filtered const& events, MyBarrelTracks const& tracks, ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) { runSelection(events, tracks, eventsMC, tracksMC); - - } void processDummy(MyEvents&) { @@ -768,16 +719,14 @@ struct AnalysisTrackSelection { PROCESS_SWITCH(AnalysisTrackSelection, processDummy, "Dummy process function", false); }; - struct AnalysisSameEventPairing { - - //Produces dileptonList; + // Produces dileptonList; // Filter based on previous components in the task Filter filterEventSelected = aod::emanalysisflags::isEventSelected == 1; Filter filterBarrelTrackSelected = aod::emanalysisflags::isBarrelSelected > 0; - + Configurable fConfigTrackCuts{"cfgTrackCuts", "", "Comma separated list of barrel track cuts"}; Configurable fConfigMCSignals{"cfgBarrelMCSignals", "", "Comma separated list of MC signals"}; Configurable fConfigMinPt{"cfgMinPt", 0., "Fiducial min Pt for MC signal"}; @@ -791,7 +740,6 @@ struct AnalysisSameEventPairing { // TODO: special AnalysisCuts to be prepared in this direction // TODO: cuts on the MC truth information to be added if needed - // 2D histos: mee and ptee Configurable fConfigMinPtee{"cfgMinPtee", 0., "min Ptee in 2D histos"}; Configurable fConfigMaxPtee{"cfgMaxPtee", 10., "max Ptee in 2D histos"}; @@ -799,39 +747,38 @@ struct AnalysisSameEventPairing { Configurable fConfigMinMee{"cfgMinMee", 0., "min Mee in 2D histos"}; Configurable fConfigMaxMee{"cfgMaxMee", 3.5, "max Mee in 2D histos"}; Configurable fConfigStepMee{"cfgStepMee", 600, "Nb of steps in Mee in 2D histos"}; - + // output lists OutputObj fOutputList{"output"}; - THashList *fMainList; // Main list - THashList *fPairList; // 2D histos for MC and reconstructed signals - THashList *fQAPairList; // QA in case on with histo manager outputs + THashList* fMainList; // Main list + THashList* fPairList; // 2D histos for MC and reconstructed signals + THashList* fQAPairList; // QA in case on with histo manager outputs // Cuts and signals - //AnalysisCompositeCut* fEventCut; // Taken from event selection part + // AnalysisCompositeCut* fEventCut; // Taken from event selection part std::vector fTrackCuts; // list of track cuts - std::vector fMCSignals; // list of signals with one prong to be checked: ULS 2D histos - - //2D histo vectors + std::vector fMCSignals; // list of signals with one prong to be checked: ULS 2D histos + + // 2D histo vectors std::vector fHistGenPair; std::vector fHistGenSmearedPair; std::vector fHistRecPair; - // Binning + // Binning std::vector fPteeBins; std::vector fMeeBins; // QA: to be defined - - + void init(o2::framework::InitContext& context) { printf("Start init\n"); - + // Create list output fMainList = new THashList; fMainList->SetOwner(kTRUE); fMainList->SetName("pairselection"); - + // Create list output for 3D eta,phi,pt fPairList = new THashList; fPairList->SetOwner(kTRUE); @@ -840,14 +787,14 @@ struct AnalysisSameEventPairing { // Binning 2D histos SetBinsLinear(fPteeBins, fConfigMinPtee, fConfigMaxPtee, fConfigStepPtee); SetBinsLinear(fMeeBins, fConfigMinMee, fConfigMaxMee, fConfigStepMee); - const int fNpteeBins = fPteeBins.size()-1; - const int fNmeeBins = fMeeBins.size()-1; - + const int fNpteeBins = fPteeBins.size() - 1; + const int fNmeeBins = fMeeBins.size() - 1; + // Event cut: taken from event selection part - //fEventCut = new AnalysisCompositeCut(true); - //TString eventCutStr = fConfigEventCuts.value; - //fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); - + // fEventCut = new AnalysisCompositeCut(true); + // TString eventCutStr = fConfigEventCuts.value; + // fEventCut->AddCut(dqcuts::GetAnalysisCut(eventCutStr.Data())); + // List of track cuts TString cutNamesStr = fConfigTrackCuts.value; if (!cutNamesStr.IsNull()) { @@ -858,126 +805,118 @@ struct AnalysisSameEventPairing { } VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill VarManager::SetDefaultVarNames(); - + // List of MC signals TString configSigNamesStr = fConfigMCSignals.value; std::unique_ptr sigNamesArray(configSigNamesStr.Tokenize(",")); for (int isig = 0; isig < sigNamesArray->GetEntries(); ++isig) { MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(sigNamesArray->At(isig)->GetName()); if (sig) { - if (sig->GetNProngs() == 2) { // only 2 prong signals - fMCSignals.push_back(*sig); - } - // List of signal to be checked + if (sig->GetNProngs() == 2) { // only 2 prong signals + fMCSignals.push_back(*sig); + } + // List of signal to be checked } } - //Configure 2D histograms + // Configure 2D histograms // Create List with generated particles TList* Generated = new TList(); Generated->SetOwner(); Generated->SetName("Generated"); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH2D* th2_tmp = new TH2D(Form("Ngen_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + for (unsigned int i = 0; i < fMCSignals.size(); ++i) { + TH2D* th2_tmp = new TH2D(Form("Ngen_Pair_%s", fMCSignals.at(i).GetName()), ";m_{ee};p_{T,ee}", fNmeeBins, fMeeBins.data(), fNpteeBins, fPteeBins.data()); th2_tmp->Sumw2(); fHistGenPair.push_back(th2_tmp); Generated->Add(th2_tmp); } - + // Create List with generated+smeared particles TList* GeneratedSmeared = new TList(); GeneratedSmeared->SetName("GeneratedSmeared"); GeneratedSmeared->SetOwner(); - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH2D* th2_tmp = new TH2D(Form("NgenSmeared_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); + for (unsigned int i = 0; i < fMCSignals.size(); ++i) { + TH2D* th2_tmp = new TH2D(Form("NgenSmeared_Pair_%s", fMCSignals.at(i).GetName()), ";m_{ee};p_{T,ee}", fNmeeBins, fMeeBins.data(), fNpteeBins, fPteeBins.data()); th2_tmp->Sumw2(); fHistGenSmearedPair.push_back(th2_tmp); GeneratedSmeared->Add(th2_tmp); } - - + fPairList->Add(Generated); fPairList->Add(GeneratedSmeared); - + // Generated reconstructed lists for every cutsetting one list and every MCsignal 2 histograms with pos and neg charge - for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i){ + for (unsigned int list_i = 0; list_i < fTrackCuts.size(); ++list_i) { TList* list = new TList(); list->SetName(fTrackCuts.at(list_i).GetName()); list->SetOwner(); - - for (unsigned int i = 0; i < fMCSignals.size(); ++i){ - TH2D* th2_tmp = new TH2D(Form("Nrec_Pair_%s", fMCSignals.at(i).GetName()),";m_{ee};p_{T,ee}",fNmeeBins,fMeeBins.data(),fNpteeBins,fPteeBins.data()); - th2_tmp->Sumw2(); - th2_tmp->SetDirectory(0x0); - fHistRecPair.push_back(th2_tmp); - list->Add(th2_tmp); + + for (unsigned int i = 0; i < fMCSignals.size(); ++i) { + TH2D* th2_tmp = new TH2D(Form("Nrec_Pair_%s", fMCSignals.at(i).GetName()), ";m_{ee};p_{T,ee}", fNmeeBins, fMeeBins.data(), fNpteeBins, fPteeBins.data()); + th2_tmp->Sumw2(); + th2_tmp->SetDirectory(0x0); + fHistRecPair.push_back(th2_tmp); + list->Add(th2_tmp); } fPairList->Add(list); } fMainList->Add(fPairList); - fOutputList.setObject(fMainList); - } Preslice perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; Preslice perReducedEventTracks = aod::reducedtrack::reducedeventId; - - template < uint32_t TEventFillMap, uint32_t TEventMCFillMap, uint32_t TTrackFillMap, typename TEvents, typename TTracks, typename TEventsMC, typename TTracksMC> + + template void runPairing(TEvents const& events, TTracks const& tracks, TEventsMC const& eventsMC, TTracksMC const& tracksMC) { - std::map fMCEventLabels; int fCounters = 0; //! [0] - particle counter, [1] - event counter - - + for (auto& event : events) { if (!event.isEventSelected()) { - return; + return; } - - VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); VarManager::ResetValues(0, VarManager::kNMCParticleVariables); // fill event information which might be needed in histograms that combine track and event properties VarManager::FillEvent(event); if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - VarManager::FillEvent(event.reducedMCevent()); + VarManager::FillEvent(event.reducedMCevent()); } - + // Look if we did not already saw the collision and fill the denominator of the single electron efficiency Int_t globalindexmc = -1; if constexpr ((TEventMCFillMap & VarManager::ObjTypes::ReducedEventMC) > 0) { - auto mcEvent = event.reducedMCevent(); - globalindexmc = mcEvent.globalIndex(); + auto mcEvent = event.reducedMCevent(); + globalindexmc = mcEvent.globalIndex(); } if (!(fMCEventLabels.find(globalindexmc) != fMCEventLabels.end())) { fMCEventLabels[globalindexmc] = fCounters; fCounters++; - // skimmed data - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); - groupedMCTracks.bindInternalIndicesTo(&tracksMC); - runMCGenPair(groupedMCTracks); - } + // skimmed data + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { + auto groupedMCTracks = tracksMC.sliceBy(perReducedMcEvent, event.reducedMCevent().globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&tracksMC); + runMCGenPair(groupedMCTracks); + } } - - - //auto groupedTrackspos = posTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); - //groupedTrackspos.bindInternalIndicesTo(&tracks); - //auto groupedTracksneg = negTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); - //groupedTracksneg.bindInternalIndicesTo(&tracks); - + + // auto groupedTrackspos = posTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); + // groupedTrackspos.bindInternalIndicesTo(&tracks); + // auto groupedTracksneg = negTracks->sliceByCached(aod::reducedtrack::reducedeventId, event.globalIndex()); + // groupedTracksneg.bindInternalIndicesTo(&tracks); + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { - auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); - groupedTracks.bindInternalIndicesTo(&tracks); - runRecPair(groupedTracks,tracksMC); - } - + auto groupedTracks = tracks.sliceBy(perReducedEventTracks, event.globalIndex()); + groupedTracks.bindInternalIndicesTo(&tracks); + runRecPair(groupedTracks, tracksMC); + } } // end loop over reconstructed event } // end loop pairing function - template void runMCGenPair(TTracksMC const& groupedMCTracks) @@ -987,100 +926,93 @@ struct AnalysisSameEventPairing { for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { - if((abs(t1.pdgCode()) != 11) || (abs(t2.pdgCode()) != 11)) continue; - if(t1.pdgCode()*t2.pdgCode() > 0) continue; // ULS only - + if ((abs(t1.pdgCode()) != 11) || (abs(t2.pdgCode()) != 11)) + continue; + if (t1.pdgCode() * t2.pdgCode() > 0) + continue; // ULS only + TLorentzVector Lvec1; TLorentzVector Lvec2; - Lvec1.SetPtEtaPhiM(t1.pt(),t1.eta(),t1.phi(),masse); - Lvec2.SetPtEtaPhiM(t2.pt(),t2.eta(),t2.phi(),masse); + Lvec1.SetPtEtaPhiM(t1.pt(), t1.eta(), t1.phi(), masse); + Lvec2.SetPtEtaPhiM(t2.pt(), t2.eta(), t2.phi(), masse); TLorentzVector LvecM = Lvec1 + Lvec2; double mass = LvecM.M(); double pairpt = LvecM.Pt(); - //double opangle = Lvec1.Angle(Lvec2.Vect()); - - //printf("Check before\n"); + // double opangle = Lvec1.Angle(Lvec2.Vect()); + + // printf("Check before\n"); // Fiducial cut Bool_t genfidcut = kTRUE; - if((t1.eta()> fConfigMaxEta) || (t2.eta() > fConfigMaxEta) || (t1.eta() < fConfigMinEta) || (t2.eta() < fConfigMinEta) || (t1.pt()> fConfigMaxPt) || (t2.pt() > fConfigMaxPt) || (t1.pt() < fConfigMinPt) || (t2.pt() < fConfigMinPt)) genfidcut = kFALSE; + if ((t1.eta() > fConfigMaxEta) || (t2.eta() > fConfigMaxEta) || (t1.eta() < fConfigMinEta) || (t2.eta() < fConfigMinEta) || (t1.pt() > fConfigMaxPt) || (t2.pt() > fConfigMaxPt) || (t1.pt() < fConfigMinPt) || (t2.pt() < fConfigMinPt)) + genfidcut = kFALSE; int isig = 0; for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - if ((*sig).CheckSignal(true, groupedMCTracks, t1, t2)) { - - // not smeared - if(genfidcut) dynamic_cast(fHistGenPair.at(isig))->Fill(mass,pairpt); - // need to implement smeared - + if ((*sig).CheckSignal(true, groupedMCTracks, t1, t2)) { + + // not smeared + if (genfidcut) + dynamic_cast(fHistGenPair.at(isig))->Fill(mass, pairpt); + // need to implement smeared } } - } //end of true pairing loop + } // end of true pairing loop } // end runMCGen - - - - template < uint32_t TTrackFillMap, typename TTracks, typename TTracksMC> + template void runRecPair(TTracks const& tracks, TTracksMC const& tracksMC) { - // Loop over two track combinations uint8_t twoTrackFilter = 0; - //uint32_t dileptonFilterMap = 0; + // uint32_t dileptonFilterMap = 0; // uint32_t dileptonMcDecision = 0; - //dileptonList.reserve(1); - - - - for (auto& [t1, t2] : combinations(tracks, tracks)) { - - twoTrackFilter = uint32_t(t1.isBarrelSelected()) & uint32_t(t2.isBarrelSelected()); - - if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue - continue; - } - VarManager::FillPair(t1, t2); - - // run MC matching for this pair - uint32_t mcDecision = 0; - int isig = 0; - for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { - - if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { // for skimmed DQ model - if ((*sig).CheckSignal(true, tracksMC, t1.reducedMCTrack(), t2.reducedMCTrack())) { - mcDecision |= (uint32_t(1) << isig); - } - } - - }// end of loop MC signals - - //dileptonFilterMap = twoTrackFilter; - //dileptonMcDecision = mcDecision; - //dileptonList(event, VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), dileptonFilterMap, dileptonMcDecision); - - - for (unsigned int i = 0; i < fMCSignals.size(); i++) { - if (!(mcDecision & (uint32_t(1) << i))) { - continue; - } - for (unsigned int j = 0; j < fTrackCuts.size(); j++) { - if (twoTrackFilter & (uint8_t(1) << j)) { - dynamic_cast(fHistRecPair.at(j * fMCSignals.size() + i))->Fill(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt]); - } - } - } + // dileptonList.reserve(1); + + for (auto& [t1, t2] : combinations(tracks, tracks)) { + + twoTrackFilter = uint32_t(t1.isBarrelSelected()) & uint32_t(t2.isBarrelSelected()); + + if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + continue; + } + VarManager::FillPair(t1, t2); + + // run MC matching for this pair + uint32_t mcDecision = 0; + int isig = 0; + for (auto sig = fMCSignals.begin(); sig != fMCSignals.end(); sig++, isig++) { + + if constexpr ((TTrackFillMap & VarManager::ObjTypes::ReducedTrack) > 0) { // for skimmed DQ model + if ((*sig).CheckSignal(true, tracksMC, t1.reducedMCTrack(), t2.reducedMCTrack())) { + mcDecision |= (uint32_t(1) << isig); + } + } + + } // end of loop MC signals + + // dileptonFilterMap = twoTrackFilter; + // dileptonMcDecision = mcDecision; + // dileptonList(event, VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), dileptonFilterMap, dileptonMcDecision); + + for (unsigned int i = 0; i < fMCSignals.size(); i++) { + if (!(mcDecision & (uint32_t(1) << i))) { + continue; + } + for (unsigned int j = 0; j < fTrackCuts.size(); j++) { + if (twoTrackFilter & (uint8_t(1) << j)) { + dynamic_cast(fHistRecPair.at(j * fMCSignals.size() + i))->Fill(VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt]); + } + } + } } - } void processToEESkimmed(soa::Filtered const& events, - soa::Filtered const& tracks, - ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) + soa::Filtered const& tracks, + ReducedMCEvents const& eventsMC, ReducedMCTracks const& tracksMC) { - - runPairing< gkEventFillMap, gkMCEventFillMap, gkTrackFillMap>(events, tracks, eventsMC, tracksMC); - + runPairing(events, tracks, eventsMC, tracksMC); } - + void processDummy(MyEvents&) { // do nothing @@ -1090,7 +1022,6 @@ struct AnalysisSameEventPairing { PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy process function", false); }; - WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ @@ -1140,15 +1071,14 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses) if (classStr.Contains("DileptonsSelected")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair_barrel"); } - - } -}// end loop over histogram classes - -void SetBinsLinear(std::vector &fBins, const double min, const double max, const unsigned int steps){ +} // end loop over histogram classes + +void SetBinsLinear(std::vector& fBins, const double min, const double max, const unsigned int steps) +{ fBins.clear(); const double stepSize = (max - min) / steps; - for (unsigned int i = 0; i < steps+1; ++i){ + for (unsigned int i = 0; i < steps + 1; ++i) { fBins.push_back(i * stepSize + min); } } From 743fd947df1d112a05e64861636cab5428e0f32d Mon Sep 17 00:00:00 2001 From: Raphaelle Bailhache Date: Fri, 15 Jul 2022 01:59:37 +0200 Subject: [PATCH 4/4] Further Fixe --- PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx index 5bbd7cd3c45..a07a00cb2c8 100644 --- a/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx +++ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx @@ -512,7 +512,7 @@ struct AnalysisTrackSelection { // All reconstructed leptons matched to a 1 prong signal or MC 1 prong signal directly std::vector mcnamesreco; - for (int isig = 0; isig < fMCSignals.size(); ++isig) { + for (unsigned int isig = 0; isig < fMCSignals.size(); ++isig) { TString nameStr2 = Form("TrackBarrel_%s_%s", cut.GetName(), fMCSignals.at(isig).GetName()); mcnamesreco.push_back(nameStr2); histClassesQA += Form("%s;", nameStr2.Data()); @@ -522,7 +522,7 @@ struct AnalysisTrackSelection { // Add histogram classes for each MC signal at generated level std::vector mcnamesgen; - for (int isig = 0; isig < fMCSignals.size(); ++isig) { + for (unsigned int isig = 0; isig < fMCSignals.size(); ++isig) { TString nameStr2 = Form("MCTruthGen_%s", fMCSignals.at(isig).GetName()); mcnamesgen.push_back(nameStr2); histClassesQA += Form("%s;", nameStr2.Data());