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/emEfficiencyEE.cxx b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx new file mode 100644 index 00000000000..a07a00cb2c8 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/emEfficiencyEE.cxx @@ -0,0 +1,1084 @@ +// 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 emanalysisflags + +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); + 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; + int fEvCounters = 0; + + // First loop + + for (auto& event : events) { + 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)"); + + 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 (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()); + } + fHistNamesMCMatchedQA.push_back(mcnamesreco); + } + + // Add histogram classes for each MC signal at generated level + std::vector mcnamesgen; + 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()); + } + 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 + 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 + 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(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); + } +}