diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 389050980e2..a2ce21c7ad9 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1362,6 +1362,10 @@ class VarManager : public TObject static void FillDileptonTrackVertexing(C const& collision, T1 const& lepton1, T1 const& lepton2, T1 const& track, float* values); template static void FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float* values = nullptr, float hadronMass = 0.0f); + template + static void FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2, T3 const& hadron, float* values = nullptr, float Translow = 1. / 3, float Transhigh = 2. / 3, bool applyFitMass = false, float sidebandMass = 0.0f); + template + static void FillEnergyCorrelatorsUnfoldingTriple(T1 const& lepton1, T2 const& lepton2, T3 const& hadron, T4 const& track, T5 const& t1, float* values = nullptr); template static void FillEnergyCorrelator(T1 const& dilepton, T2 const& hadron, float* values = nullptr, float Translow = 1. / 3, float Transhigh = 2. / 3, bool applyFitMass = false, float sidebandMass = 0.0f); template @@ -5872,6 +5876,114 @@ void VarManager::FillEnergyCorrelator(T1 const& dilepton, T2 const& hadron, floa } } } + +template +void VarManager::FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2, T3 const& hadron, float* values, float Translow, float Transhigh, bool applyFitMass, float sidebandMass) +{ + float m1 = o2::constants::physics::MassElectron; + float m2 = o2::constants::physics::MassElectron; + + ROOT::Math::PtEtaPhiMVector v_lepton1(lepton1.pt(), lepton1.eta(), lepton1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v_lepton2(lepton2.pt(), lepton2.eta(), lepton2.phi(), m2); + ROOT::Math::PtEtaPhiMVector dilepton = v_lepton1 + v_lepton2; + + float dileptonmass = o2::constants::physics::MassJPsi; + if (applyFitMass) { + dileptonmass = dilepton.mass(); + } + if (applyFitMass && sidebandMass > 0) { + dileptonmass = sidebandMass; + } + + if (fgUsedVars[kCosChi] || fgUsedVars[kECWeight] || fgUsedVars[kCosTheta] || fgUsedVars[kEWeight_before] || fgUsedVars[kPtDau] || fgUsedVars[kEtaDau] || fgUsedVars[kPhiDau] || fgUsedVars[kCosChi_randomPhi_trans] || fgUsedVars[kCosChi_randomPhi_toward] || fgUsedVars[kCosChi_randomPhi_away]) { + values[kdileptonmass] = dileptonmass; + ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dileptonmass); + ROOT::Math::PtEtaPhiMVector v2(hadron.pt(), hadron.eta(), hadron.phi(), o2::constants::physics::MassPionCharged); + values[kCosChi] = LorentzTransformJpsihadroncosChi("coschi", v1, v2); + float E_boost = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2); + values[kECWeight] = E_boost / v1.M(); + values[kCosTheta] = LorentzTransformJpsihadroncosChi("costheta", v1, v2); + values[kEWeight_before] = v2.Pt() / v1.M(); + values[kPtDau] = v2.pt(); + values[kEtaDau] = v2.eta(); + values[kPhiDau] = RecoDecay::constrainAngle(v2.phi(), -o2::constants::math::PIHalf); + + float deltaphi = RecoDecay::constrainAngle(v1.phi() - v2.phi(), -o2::constants::math::PI); + values[kDeltaPhi] = RecoDecay::constrainAngle(v1.phi() - v2.phi(), -o2::constants::math::PIHalf); + values[kDeltaEta] = v1.eta() - v2.eta(); + values[kCosChi_randomPhi_trans] = -999.9f; + values[kCosChi_randomPhi_toward] = -999.9f; + values[kCosChi_randomPhi_away] = -999.9f; + + values[kdeltaphi_randomPhi_trans] = -999.9f; + values[kdeltaphi_randomPhi_toward] = -999.9f; + values[kdeltaphi_randomPhi_away] = -999.9f; + + float randomPhi_trans = -o2::constants::math::PIHalf; + float randomPhi_toward = -o2::constants::math::PIHalf; + float randomPhi_away = -o2::constants::math::PIHalf; + + if ((deltaphi > -Transhigh * TMath::Pi() && deltaphi < -Translow * TMath::Pi()) || (deltaphi > Translow * TMath::Pi() && deltaphi < Transhigh * TMath::Pi())) { + randomPhi_trans = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); + randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); + randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf); + + ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, o2::constants::physics::MassPionCharged); + values[kCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans); + values[kWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M(); + + ROOT::Math::PtEtaPhiMVector v2_randomPhi_toward(v2.pt(), v2.eta(), randomPhi_toward, o2::constants::physics::MassPionCharged); + values[kCosChi_randomPhi_toward] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_toward); + values[kWeight_randomPhi_toward] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_toward) / v1.M(); + + ROOT::Math::PtEtaPhiMVector v2_randomPhi_away(v2.pt(), v2.eta(), randomPhi_away, o2::constants::physics::MassPionCharged); + values[kCosChi_randomPhi_away] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_away); + values[kWeight_randomPhi_away] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_away) / v1.M(); + + values[kdeltaphi_randomPhi_trans] = RecoDecay::constrainAngle(v1.phi() - randomPhi_trans, -o2::constants::math::PIHalf); + values[kdeltaphi_randomPhi_toward] = RecoDecay::constrainAngle(v1.phi() - randomPhi_toward, -o2::constants::math::PIHalf); + values[kdeltaphi_randomPhi_away] = RecoDecay::constrainAngle(v1.phi() - randomPhi_away, -o2::constants::math::PIHalf); + } + } +} + +template +void VarManager::FillEnergyCorrelatorsUnfoldingTriple(T1 const& lepton1, T2 const& lepton2, T3 const& hadron, T4 const& track, T5 const& t1, float* values) +{ + if (fgUsedVars[kMCCosChi_gen] || fgUsedVars[kMCWeight_gen] || fgUsedVars[kMCdeltaeta_gen] || fgUsedVars[kMCCosChi_rec] || fgUsedVars[kMCWeight_rec] || fgUsedVars[kMCdeltaeta_rec]) { + // energy correlators + + float m1 = o2::constants::physics::MassElectron; + float m2 = o2::constants::physics::MassElectron; + + ROOT::Math::PtEtaPhiMVector v_lepton1(lepton1.pt(), lepton1.eta(), lepton1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v_lepton2(lepton2.pt(), lepton2.eta(), lepton2.phi(), m2); + ROOT::Math::PtEtaPhiMVector dilepton = v_lepton1 + v_lepton2; + + float MassHadron; + if constexpr (pairType == kJpsiHadronMass) { + MassHadron = TMath::Sqrt(t1.e() * t1.e() - t1.p() * t1.p()); + } + if constexpr (pairType == kJpsiPionMass) { + MassHadron = o2::constants::physics::MassPionCharged; + } + ROOT::Math::PtEtaPhiMVector v1_gen(track.pt(), track.eta(), track.phi(), o2::constants::physics::MassJPsi); + ROOT::Math::PtEtaPhiMVector v2_gen(t1.pt(), t1.eta(), t1.phi(), MassHadron); + float E_boost_gen = LorentzTransformJpsihadroncosChi("weight_boost", v1_gen, v2_gen); + float CosChi_gen = LorentzTransformJpsihadroncosChi("coschi", v1_gen, v2_gen); + values[kMCCosChi_gen] = CosChi_gen; + values[kMCWeight_gen] = E_boost_gen / o2::constants::physics::MassJPsi; + values[kMCdeltaeta_gen] = track.eta() - t1.eta(); + + ROOT::Math::PtEtaPhiMVector v1_rec(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); + ROOT::Math::PtEtaPhiMVector v2_rec(hadron.pt(), hadron.eta(), hadron.phi(), o2::constants::physics::MassPionCharged); + values[kMCCosChi_rec] = LorentzTransformJpsihadroncosChi("coschi", v1_rec, v2_rec); + float E_boost_rec = LorentzTransformJpsihadroncosChi("weight_boost", v1_rec, v2_rec); + values[kMCWeight_rec] = E_boost_rec / v1_rec.M(); + values[kMCdeltaeta_rec] = dilepton.eta() - hadron.eta(); + } +} + template void VarManager::FillDileptonPhoton(T1 const& dilepton, T2 const& photon, float* values) { diff --git a/PWGDQ/Tasks/CMakeLists.txt b/PWGDQ/Tasks/CMakeLists.txt index cc660b53202..4bb0cf6d714 100644 --- a/PWGDQ/Tasks/CMakeLists.txt +++ b/PWGDQ/Tasks/CMakeLists.txt @@ -39,6 +39,11 @@ o2physics_add_dpl_workflow(efficiency-with-assoc-direct PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(energy-correlator-direct + SOURCES dqEnergyCorrelator_direct.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(filter-pp SOURCES filterPP.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore diff --git a/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx b/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx new file mode 100644 index 00000000000..0f62c150ae0 --- /dev/null +++ b/PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx @@ -0,0 +1,1013 @@ +// 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. +// +// Contact: iarsene@cern.ch, i.c.arsene@fys.uio.no +// Configurable workflow for running several DQ or other PWG analyses + +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "PWGDQ/Core/MCSignal.h" +#include "PWGDQ/Core/MCSignalLibrary.h" +#include "PWGDQ/Core/MixingHandler.h" +#include "PWGDQ/Core/MixingLibrary.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" + +#include "Common/Core/PID/PIDTOFParamService.h" +#include "Common/Core/TableHelper.h" +#include "Common/DataModel/CollisionAssociationTables.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/McCollisionExtra.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "Field/MagneticField.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +#include "TGeoGlobalMagField.h" +#include +#include +#include +#include +#include +#include + +#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; + +using MyEvents = soa::Join; +using MyBarrelTracksWithCov = soa::Join; + +// bit maps used for the Fill functions of the VarManager +constexpr static uint32_t gkEventFillMapWithMults = VarManager::ObjTypes::BC | VarManager::ObjTypes::Collision | VarManager::ObjTypes::CollisionMult | VarManager::ObjTypes::CollisionMultExtra; +constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::Track | VarManager::ObjTypes::TrackExtra | VarManager::ObjTypes::TrackDCA | VarManager::ObjTypes::TrackCov | VarManager::ObjTypes::TrackPID; + +// Forward declarations +void DefineHistograms(HistogramManager* histMan, TString histClasses, const char* histGroups); // defines histograms for all tasks + +struct AnalysisEnergyCorrelator { + OutputObj fOutputList{"output"}; + + struct : ConfigurableGroup { // Event selection configurables + Configurable> fConfigZBins{"cfgZBins", std::vector{-10.0, -5.0, 0.0, 5.0, 10.0}, "Z vertex bins for mixing"}; + Configurable> fConfigMultBins{"cfgMultBins", std::vector{0.0, 20.0, 40.0, 60.0, 100.0}, "Multiplicity bins for mixing"}; + Configurable fConfigEventCuts{"cfgEventCuts", "eventStandard", "Event selection"}; + Configurable fConfigEventCutsJSON{"cfgEventCutsJSON", "", "Additional event cuts in JSON"}; + Configurable fConfigAddEventHistogram{"cfgAddEventHistogram", "", "Event histograms"}; + Configurable fConfigAddEventMCHistogram{"cfgAddEventMCHistogram", "generator", "MC Event histograms"}; + Configurable fConfigMixingDepth{"cfgMixingDepth", 5, "Event mixing pool depth"}; + Configurable fConfigEventfilterVtz{"cfgEventfilterVtz", 10.0, "Event filter Vtz"}; + Configurable fConfigEventQA{"cfgEventQA", false, "If true, fill Event QA histograms"}; + } fConfigEventOptions; + + struct : ConfigurableGroup { // Track selection configurables + Configurable fConfigTrackCuts{"cfgTrackCuts", "electronSelection1_ionut", "Comma separated list of barrel track cuts for electrons"}; + Configurable fConfigTrackCutsJSON{"cfgTrackCutsJSON", "", "Additional track cuts in JSON"}; + Configurable fConfigAddTrackHistogram{"cfgAddTrackHistogram", "", "Track histograms"}; + Configurable fConfigTrackQA{"cfgTrackQA", false, "If true, fill Track QA histograms"}; + } fConfigTrackOptions; + + struct : ConfigurableGroup { // Pair selection configurables + Configurable fConfigJpsiMassMin{"cfgJpsiMassMin", 2.8, "J/psi mass minimum"}; + Configurable fConfigJpsiMassMax{"cfgJpsiMassMax", 3.3, "J/psi mass maximum"}; + Configurable fConfigJpsiPtMin{"cfgJpsiPtMin", 0.0, "J/psi pt minimum"}; + Configurable fConfigJpsiPtMax{"cfgJpsiPtMax", 100.0, "J/psi pt maximum"}; + Configurable fConfigJpsiRapMax{"cfgJpsiRapMax", 0.9, "J/psi rapidity maximum"}; + Configurable fConfigAddSEPHistogram{"cfgAddSEPHistogram", "", "Comma separated list of histograms"}; + Configurable recSignals{"cfgMCRecSignals", "eeFromJpsi", "Comma separated list of MC signals (reconstructed)"}; + Configurable recSignalsJSON{"cfgMCRecSignalsJSON", "", "Comma separated list of MC signals (reconstructed) via JSON"}; + } fConfigPairOptions; + + struct : ConfigurableGroup { // Dilepton selection configurables + Configurable fConfigHadronCuts{"cfgHadronCuts", "NoPID", "Comma separated list of hadron track cuts"}; + Configurable fConfigHadronCutsJSON{"cfgHadronCutsJSON", "", "Additional hadron cuts in JSON"}; + Configurable fConfigApplyMassEC{"cfgApplyMassEC", false, "Apply fit mass for sideband for the energy correlator study"}; + Configurable> fConfigSavelessevents{"cfgSavelessevents", std::vector{1, 0}, "Save less events for the energy correlator study"}; + Configurable> fConfigTransRange{"cfgTransRange", std::vector{0.333333, 0.666667}, "Transverse region for the energy correlstor analysis"}; + Configurable fConfigAddDileptonHadronHistogram{"cfgAddDileptonHadronHistogram", "", "Dilepton-hadron histograms"}; + Configurable fConfigMCRecSignals{"cfgMCRecDileptonHadronSignals", "", "Comma separated list of MC signals (reconstructed)"}; + Configurable fConfigMCGenSignals{"cfgMCGenDileptonHadronSignals", "", "Comma separated list of MC signals (generated)"}; + Configurable fConfigMCRecSignalsJSON{"cfgMCRecDileptonHadronSignalsJSON", "", "Additional list of MC signals (reconstructed) via JSON"}; + Configurable fConfigMCGenSignalsJSON{"cfgMCGenDileptonHadronSignalsJSON", "", "Comma separated list of MC signals (generated) via JSON"}; + Configurable fConfigMCGenHadronEtaAbs{"cfgMCGenHadronEtaAbs", 0.9f, "eta abs range for the hadron"}; + Configurable fConfigMCGenHadronPtMin{"cfgMCGenHadronPtMin", 0.1f, "minimum pt for the hadron"}; + } fConfigDileptonHadronOptions; + + // Histogram configurables + Configurable fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON"}; + + // CCDB configurables + Configurable fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "CCDB url"}; + Configurable fConfigNoLaterThan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "CCDB timestamp"}; + + // Member variables + HistogramManager* fHistMan = nullptr; + MixingHandler* fMixHandler = nullptr; + + AnalysisCompositeCut* fEventCut = nullptr; + std::vector fTrackCuts; // Electron cuts + std::vector fHadronCuts; // Hadron cuts + + std::vector fTrackCutNames; + std::vector fHadronCutNames; + std::vector fHistNamesReco; + + std::map> fTrackHistNames; + std::map> fBarrelHistNamesMCmatched; + std::map fSelMap; + + std::vector fRecMCSignals; // MC signals for reconstructed pairs + std::vector fGenMCSignals; + std::vector fRecMCTripleSignals; // MC signals for reconstructed triples + + Service fCCDB; + int fCurrentRun = -1; + + // Preslice for association table + Preslice preslice = aod::track_association::collisionId; + + using MixingBinning = ColumnBinningPolicy; + std::unique_ptr fMixingBinning; + + void init(o2::framework::InitContext& context) + { + std::vector zBins = fConfigEventOptions.fConfigZBins.value; + zBins.insert(zBins.begin(), VARIABLE_WIDTH); + + std::vector multBins = fConfigEventOptions.fConfigMultBins.value; + multBins.insert(multBins.begin(), VARIABLE_WIDTH); + + fMixingBinning = std::make_unique(std::array, 2>{zBins, multBins}, true); + + bool isBarrelME = context.mOptions.get("processBarrelMixedEvent"); + bool isMCGen_energycorrelators = context.mOptions.get("processMCGenEnergyCorrelators") || context.mOptions.get("processMCGenEnergyCorrelatorsPion"); + bool isMCGen_energycorrelatorsME = context.mOptions.get("processMCGenEnergyCorrelatorsME") || context.mOptions.get("processMCGenEnergyCorrelatorsPionME"); + + if (context.mOptions.get("processDummy")) { + return; + } + VarManager::SetDefaultVarNames(); + + // Setup Event Cuts + fEventCut = new AnalysisCompositeCut(true); + TString eventCutStr = fConfigEventOptions.fConfigEventCuts.value; + if (eventCutStr != "") { + AnalysisCut* cut = dqcuts::GetAnalysisCut(eventCutStr.Data()); + if (cut != nullptr) { + fEventCut->AddCut(cut); + } + } + // Additional cuts via JSON + TString eventCutJSONStr = fConfigEventOptions.fConfigEventCutsJSON.value; + if (eventCutJSONStr != "") { + std::vector jsonCuts = dqcuts::GetCutsFromJSON(eventCutJSONStr.Data()); + for (auto& cutIt : jsonCuts) { + fEventCut->AddCut(cutIt); + } + } + + // Setup Electron Track Cuts + TString trackCutStr = fConfigTrackOptions.fConfigTrackCuts.value; + if (!trackCutStr.IsNull()) { + std::unique_ptr objArrayTrack(trackCutStr.Tokenize(",")); + for (int icut = 0; icut < objArrayTrack->GetEntries(); ++icut) { + fTrackCuts.push_back(dqcuts::GetCompositeCut(objArrayTrack->At(icut)->GetName())); + fTrackCutNames.push_back(objArrayTrack->At(icut)->GetName()); + } + } + + TString trackCutsJSON = fConfigTrackOptions.fConfigTrackCutsJSON.value; + if (trackCutsJSON != "") { + std::vector addTrackCuts = dqcuts::GetCutsFromJSON(trackCutsJSON.Data()); + for (auto& t : addTrackCuts) { + fTrackCuts.push_back(reinterpret_cast(t)); + fTrackCutNames.push_back(t->GetName()); + } + } + + // Setting the MC rec signal names + TString sigNamesStr = fConfigPairOptions.recSignals.value; + std::unique_ptr objRecSigArray(sigNamesStr.Tokenize(",")); + for (int isig = 0; isig < objRecSigArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objRecSigArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() != 2) { // NOTE: 2-prong signals required + continue; + } + fRecMCSignals.push_back(sig); + } + } + + TString sigNamesHadronStr = fConfigDileptonHadronOptions.fConfigMCRecSignals.value; + std::unique_ptr objRecSigTripleArray(sigNamesHadronStr.Tokenize(",")); + if (!sigNamesHadronStr.IsNull()) { + for (int isig = 0; isig < objRecSigTripleArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objRecSigTripleArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() != 3) { + LOG(fatal) << "Signal at reconstructed level requested (" << sig->GetName() << ") " << "does not have 3 prongs! Fix it"; + } + fRecMCTripleSignals.push_back(sig); + } else { + LOG(fatal) << "Signal at reconstructed level requested (" << objRecSigTripleArray->At(isig)->GetName() << ") " << "could not be retrieved from the library! -> skipped"; + } + } + } + + // Add the MCSignals from the JSON config + TString addMCSignalsStr = fConfigPairOptions.recSignalsJSON.value; + if (addMCSignalsStr != "") { + std::vector addMCSignals = dqmcsignals::GetMCSignalsFromJSON(addMCSignalsStr.Data()); + for (auto& mcIt : addMCSignals) { + if (mcIt->GetNProngs() != 2) { // NOTE: only 2 prong signals + continue; + } + fRecMCSignals.push_back(mcIt); + } + } + + // Add the reco MCSignals from the JSON config + TString addMCTripleSignalsStr = fConfigDileptonHadronOptions.fConfigMCRecSignalsJSON.value; + if (addMCTripleSignalsStr != "") { + std::vector addMCTripleSignals = dqmcsignals::GetMCSignalsFromJSON(addMCTripleSignalsStr.Data()); + for (auto& mcIt : addMCTripleSignals) { + if (mcIt->GetNProngs() != 3) { + LOG(fatal) << "Signal at reconstructed level requested (" << mcIt->GetName() << ") " << "does not have 3 prongs! Fix it"; + } + fRecMCTripleSignals.push_back(mcIt); + } + } + + // Setup Hadron Cuts + TString hadronCutStr = fConfigDileptonHadronOptions.fConfigHadronCuts.value; + if (!hadronCutStr.IsNull()) { + std::unique_ptr objArrayHadron(hadronCutStr.Tokenize(",")); + for (int icut = 0; icut < objArrayHadron->GetEntries(); ++icut) { + TString cutName = objArrayHadron->At(icut)->GetName(); + fHadronCuts.push_back(dqcuts::GetCompositeCut(cutName.Data())); + fHadronCutNames.push_back(cutName); + } + } + TString hadronCutsJSON = fConfigDileptonHadronOptions.fConfigHadronCutsJSON.value; + if (hadronCutsJSON != "") { + std::vector addHadronCuts = dqcuts::GetCutsFromJSON(hadronCutsJSON.Data()); + for (auto& t : addHadronCuts) { + fHadronCuts.push_back(reinterpret_cast(t)); + fHadronCutNames.push_back(t->GetName()); + } + } + + // Add histogram classes for each specified MCsignal at the generator level + // TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function + TString sigGenNamesStr = fConfigDileptonHadronOptions.fConfigMCGenSignals.value; + std::unique_ptr objGenSigArray(sigGenNamesStr.Tokenize(",")); + for (int isig = 0; isig < objGenSigArray->GetEntries(); isig++) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objGenSigArray->At(isig)->GetName()); + if (sig) { + fGenMCSignals.push_back(sig); + } + } + + // Add the MCSignals from the JSON config + TString addMCSignalsGenStr = fConfigDileptonHadronOptions.fConfigMCGenSignalsJSON.value; + if (addMCSignalsGenStr != "") { + std::vector addMCSignals = dqmcsignals::GetMCSignalsFromJSON(addMCSignalsGenStr.Data()); + for (auto& mcIt : addMCSignals) { + if (mcIt->GetNProngs() > 2) { // NOTE: only 2 prong signals + continue; + } + fGenMCSignals.push_back(mcIt); + } + } + + VarManager::SetUseVars(AnalysisCut::fgUsedVars); + + fHistMan = new HistogramManager("analysisHistos", "", VarManager::kNVars); + fHistMan->SetUseDefaultVariableNames(true); + fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + + // Setup Histograms + if (fConfigEventOptions.fConfigEventQA) { + DefineHistograms(fHistMan, "TimeFrameStats;Event_BeforeCuts;Event_AfterCuts;", fConfigEventOptions.fConfigAddEventHistogram.value.data()); + DefineHistograms(fHistMan, "EventsMC", fConfigEventOptions.fConfigAddEventMCHistogram.value.data()); // mc + } + + if (fConfigTrackOptions.fConfigTrackQA) { + TString histClasses = "AssocsBarrel_BeforeCuts;"; + // Configure histogram classes for each track cut; + // Add histogram classes for each track cut and for each requested MC signal (reconstructed tracks with MC truth) + for (auto& cut : fTrackCuts) { + TString nameStr = Form("AssocsBarrel_%s", cut->GetName()); + fHistNamesReco.push_back(nameStr); + histClasses += Form("%s;", nameStr.Data()); + } + DefineHistograms(fHistMan, histClasses.Data(), fConfigTrackOptions.fConfigAddTrackHistogram.value.data()); + } + TString histNames = ""; + // check that the barrel track cuts array required in this task is not empty + if (!trackCutStr.IsNull()) { + // tokenize and loop over the barrel cuts produced by the barrel track selection task + std::unique_ptr objArray(trackCutStr.Tokenize(",")); + for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + TString tempStr = objArray->At(icut)->GetName(); + // if the current barrel selection cut is required in this task, then switch on the corresponding bit in the mask + // and assign histogram directories + + // assign the pair hist directories for the current cut + std::vector names = { + Form("PairsBarrelSEPM_%s", objArray->At(icut)->GetName()), + Form("PairsBarrelSEPP_%s", objArray->At(icut)->GetName()), + Form("PairsBarrelSEMM_%s", objArray->At(icut)->GetName())}; + for (auto& n : names) { + histNames += Form("%s;", n.Data()); + } + fTrackHistNames[icut] = names; + // assign hist directories for the MC matched pairs for each (track cut,MCsignal) combination + if (!sigNamesStr.IsNull()) { + for (size_t isig = 0; isig < fRecMCSignals.size(); isig++) { + auto sig = fRecMCSignals.at(isig); + names = { + Form("PairsBarrelSEPM_%s_%s", objArray->At(icut)->GetName(), sig->GetName()), + Form("PairsBarrelSEPP_%s_%s", objArray->At(icut)->GetName(), sig->GetName()), + Form("PairsBarrelSEMM_%s_%s", objArray->At(icut)->GetName(), sig->GetName())}; + for (auto& n : names) { + histNames += Form("%s;", n.Data()); + } + fBarrelHistNamesMCmatched.try_emplace(icut * fRecMCSignals.size() + isig, names); + } // end loop over MC signals + } + } + DefineHistograms(fHistMan, histNames.Data(), fConfigPairOptions.fConfigAddSEPHistogram.value.data()); + } + + for (size_t iCutTrack = 0; iCutTrack < fTrackCutNames.size(); iCutTrack++) { + for (size_t iCutHadron = 0; iCutHadron < fHadronCutNames.size(); iCutHadron++) { + DefineHistograms(fHistMan, Form("DileptonTrack_%s_%s", fTrackCutNames[iCutTrack].Data(), fHadronCutNames[iCutHadron].Data()), fConfigDileptonHadronOptions.fConfigAddDileptonHadronHistogram.value.data()); + for (auto& sig : fRecMCTripleSignals) { + DefineHistograms(fHistMan, Form("DileptonTrackMCMatched_%s_%s_%s", fTrackCutNames[iCutTrack].Data(), fHadronCutNames[iCutHadron].Data(), sig->GetName()), fConfigDileptonHadronOptions.fConfigAddDileptonHadronHistogram.value.data()); + if (isBarrelME) { + DefineHistograms(fHistMan, Form("DileptonTrackMCMatchedME_%s_%s_%s", fTrackCutNames[iCutTrack].Data(), fHadronCutNames[iCutHadron].Data(), sig->GetName()), fConfigDileptonHadronOptions.fConfigAddDileptonHadronHistogram.value.data()); + } + } + } + } + + for (auto& sig : fGenMCSignals) { + if (sig->GetNProngs() == 1) { + if (isMCGen_energycorrelators) { + DefineHistograms(fHistMan, Form("MCTruthGenSel_%s", sig->GetName()), ""); + DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelators_%s", sig->GetName()), ""); + } + if (isMCGen_energycorrelatorsME) { + DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelatorsME_%s", sig->GetName()), ""); + } + } + } + + dqhistograms::AddHistogramsFromJSON(fHistMan, fConfigAddJSONHistograms.value.c_str()); // aditional histograms via JSON + + VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + fOutputList.setObject(fHistMan->GetMainHistogramList()); + + fCCDB->setURL(fConfigCcdbUrl.value); + fCCDB->setCaching(true); + fCCDB->setLocalObjectValidityChecking(); + fCCDB->setCreatedNotAfter(fConfigNoLaterThan.value); + } + + template + void runDileptonHadron(TTrack1 const& track1, TTrack2 const& track2, int iEleCut, + THadron const& hadron, TEvent const& event, aod::McParticles const& /*mcParticles*/) + { + VarManager::ResetValues(0, VarManager::kNVars); // reset variables before filling + VarManager::FillEvent(event); + VarManager::FillTrack(hadron); + VarManager::FillTrackCollision(hadron, event); + + // Check that hadron is not one of the dilepton legs + if (hadron.globalIndex() == track1.globalIndex() || hadron.globalIndex() == track2.globalIndex()) { + return; + } + + if (!track1.has_mcParticle() || !track2.has_mcParticle() || !hadron.has_mcParticle()) { + return; + } + auto hadronMC = hadron.mcParticle(); + auto lepton1MC = track1.mcParticle(); + auto lepton2MC = track2.mcParticle(); + uint32_t mcDecision = 0; + int isig = 0; + for (auto sig = fRecMCTripleSignals.begin(); sig != fRecMCTripleSignals.end(); sig++, isig++) { + if ((*sig)->CheckSignal(true, lepton1MC, lepton2MC, hadronMC)) { + mcDecision |= (static_cast(1) << isig); + } + } + + auto motherParticle = lepton1MC.template mothers_first_as(); + // Fill dilepton-hadron variables + std::vector fTransRange = fConfigDileptonHadronOptions.fConfigTransRange; + VarManager::FillEnergyCorrelatorTriple(track1, track2, hadron, VarManager::fgValues, fTransRange[0], fTransRange[1], fConfigDileptonHadronOptions.fConfigApplyMassEC.value); + VarManager::FillEnergyCorrelatorsUnfoldingTriple(track1, track2, hadron, motherParticle, hadronMC, VarManager::fgValues); + + int iHadronCut = 0; + for (auto hCut = fHadronCuts.begin(); hCut != fHadronCuts.end(); hCut++, iHadronCut++) { + if (!(*hCut)->IsSelected(VarManager::fgValues)) { + continue; + } + // Fill the corresponding histogram + if (!MixedEvent) { + fHistMan->FillHistClass( + Form("DileptonTrack_%s_%s", fTrackCutNames[iEleCut].Data(), fHadronCutNames[iHadronCut].Data()), + VarManager::fgValues); + } + for (uint32_t isig = 0; isig < fRecMCTripleSignals.size(); isig++) { + if (mcDecision & (static_cast(1) << isig)) { + if (!MixedEvent) { + fHistMan->FillHistClass(Form("DileptonTrackMCMatched_%s_%s_%s", fTrackCutNames[iEleCut].Data(), fHadronCutNames[iHadronCut].Data(), fRecMCTripleSignals[isig]->GetName()), VarManager::fgValues); + } + if (MixedEvent) { + fHistMan->FillHistClass(Form("DileptonTrackMCMatchedME_%s_%s_%s", fTrackCutNames[iEleCut].Data(), fHadronCutNames[iHadronCut].Data(), fRecMCTripleSignals[isig]->GetName()), VarManager::fgValues); + } + } // end loop over MC signals + } + } + } + + // Template function to run same event pairing (barrel-barrel, muon-muon, barrel-muon) + template + void runSameEventPairing(TTrack1 const& t1, TTrack2 const& t2, int iEleCut, aod::McParticles const& /*mcParticles*/) + { + std::map> histNames = fTrackHistNames; + std::map> histNamesMC = fBarrelHistNamesMCmatched; + int sign1 = t1.sign(); + int sign2 = t2.sign(); + uint32_t mcDecision = static_cast(0); + // run MC matching for this pair + int isig = 0; + mcDecision = 0; + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + if (t1.has_mcParticle() && t2.has_mcParticle()) { + if ((*sig)->CheckSignal(true, t1.mcParticle(), t2.mcParticle())) { + mcDecision |= (static_cast(1) << isig); + } + } + } // end loop over MC signals + + VarManager::FillPair(t1, t2); + + if (sign1 * sign2 < 0) { // +- pairs + fHistMan->FillHistClass(histNames[iEleCut][0].Data(), VarManager::fgValues); // reconstructed, unmatched + for (size_t isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[iEleCut * fRecMCSignals.size() + isig][0].Data(), VarManager::fgValues); // matched signal + } + } + } else { + if (sign1 > 0) { // ++ pairs + fHistMan->FillHistClass(histNames[iEleCut][1].Data(), VarManager::fgValues); + for (size_t isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[iEleCut * fRecMCSignals.size() + isig][1].Data(), VarManager::fgValues); + } + } + } else { // -- pairs + fHistMan->FillHistClass(histNames[iEleCut][2].Data(), VarManager::fgValues); + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[iEleCut * fRecMCSignals.size() + isig][2].Data(), VarManager::fgValues); + } + } + } + } + } + + void processBarrel(MyEvents const& events, aod::TrackAssoc const& assocs, MyBarrelTracksWithCov const& /*tracks*/, soa::Join const& mcEvents, aod::McParticles const& mcParticles, BCsWithTimestamps const& bcs) + { + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillTimeFrame(bcs); + VarManager::FillTimeFrame(events); + VarManager::FillTimeFrame(mcEvents); + + if (events.size() == 0) + return; + + // CCDB initialization + if (fCurrentRun != bcs.begin().runNumber()) { + fCurrentRun = bcs.begin().runNumber(); + } + + if (fConfigEventOptions.fConfigEventQA) { + fHistMan->FillHistClass("TimeFrameStats", VarManager::fgValues); + } + + for (auto& event : mcEvents) { + // Reset the fValues array and fill event observables + VarManager::FillEvent(event); + if (fConfigEventOptions.fConfigEventQA) { + fHistMan->FillHistClass("EventsMC", VarManager::fgValues); + } + } + + // Event loop + for (auto& event : events) { + // Fill event variables first + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::FillEvent(event); + + if (fConfigEventOptions.fConfigEventQA) { + fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); + } + + // Event selection + if (!fEventCut->IsSelected(VarManager::fgValues)) + continue; + + if (fConfigEventOptions.fConfigEventQA) { + fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + } + + // saveless events for the energy correlator analysis + std::vector fSavelessevents = fConfigDileptonHadronOptions.fConfigSavelessevents.value; + if (fSavelessevents[0] > 1 && event.globalIndex() % fSavelessevents[0] == fSavelessevents[1]) { + continue; + } + + // Get associated tracks for this event + auto groupedAssocs = assocs.sliceBy(preslice, event.globalIndex()); + + // Triple loop: track1 (electron) x track2 (electron) x hadron + for (auto& a1 : groupedAssocs) { + auto t1 = a1.template track_as(); + + uint32_t filter1 = 0; + // Fill track variables + VarManager::FillTrack(t1); + VarManager::FillTrackCollision(t1, event); + if (t1.has_mcParticle()) { + VarManager::FillTrackMC(mcParticles, t1.mcParticle()); + } + + if (fConfigTrackOptions.fConfigTrackQA) { + fHistMan->FillHistClass("AssocsBarrel_BeforeCuts", VarManager::fgValues); + } + + // Apply electron cuts and fill histograms + int iCut1 = 0; + for (auto cut1 = fTrackCuts.begin(); cut1 != fTrackCuts.end(); cut1++, iCut1++) { + if ((*cut1)->IsSelected(VarManager::fgValues)) { + filter1 |= (static_cast(1) << iCut1); + if (fConfigTrackOptions.fConfigTrackQA) { + fHistMan->FillHistClass(fHistNamesReco[iCut1], VarManager::fgValues); + } + } + } + + // Check opposite charge with t2 + for (auto& a2 : groupedAssocs) { + auto t2 = a2.template track_as(); + + // Avoid double counting: use track globalIndex + if (t2.globalIndex() <= t1.globalIndex()) { + continue; + } + + // Fill track variables for t2 (only once per t2) + VarManager::FillTrack(t2); + VarManager::FillTrackCollision(t2, event); + + // Compute filter2: which cuts t2 passes + uint32_t filter2 = 0; + int iCut2 = 0; + for (auto cut2 = fTrackCuts.begin(); cut2 != fTrackCuts.end(); cut2++, iCut2++) { + if ((*cut2)->IsSelected(VarManager::fgValues)) { + filter2 |= (static_cast(1) << iCut2); + } + } + + // Both tracks must pass at least one common cut + uint32_t twoTrackFilter = filter1 & filter2; + if (!twoTrackFilter) { + continue; + } + + // Fill pair histograms for all cuts that both tracks pass + for (size_t iCut = 0; iCut < fTrackCuts.size(); iCut++) { + if (twoTrackFilter & (static_cast(1) << iCut)) { + runSameEventPairing(t1, t2, iCut, mcParticles); + } + } + + float mass = VarManager::fgValues[VarManager::kMass]; + float pt = VarManager::fgValues[VarManager::kPt]; + float rap = VarManager::fgValues[VarManager::kRap]; + + // Apply J/psi cuts + if (mass < fConfigPairOptions.fConfigJpsiMassMin.value || mass > fConfigPairOptions.fConfigJpsiMassMax.value || + pt < fConfigPairOptions.fConfigJpsiPtMin.value || pt > fConfigPairOptions.fConfigJpsiPtMax.value || + std::abs(rap) > fConfigPairOptions.fConfigJpsiRapMax.value) { + continue; + } + + if (t1.sign() * t2.sign() >= 0) { + continue; // Must be opposite charge + } + + // correlate J/psi with hadrons + for (auto& aHadron : groupedAssocs) { + auto hadron = aHadron.template track_as(); + // Process dilepton-hadron correlation for each common cut + for (size_t iCut = 0; iCut < fTrackCuts.size(); iCut++) { + if (twoTrackFilter & (static_cast(1) << iCut)) { + runDileptonHadron(t1, t2, iCut, hadron, event, mcParticles); + } + } + } // end hadron loop + } // end track2 loop + } // end track1 loop + } // end event loop + } + + Filter eventFilter = nabs(aod::collision::posZ) < fConfigEventOptions.fConfigEventfilterVtz && aod::evsel::sel8 == true; + void processBarrelMixedEvent(soa::Filtered& events, aod::TrackAssoc const& assocs, MyBarrelTracksWithCov const& /*tracks*/, aod::McCollisions const& /*mcCollisions*/, aod::McParticles const& mcParticles, BCsWithTimestamps const& bcs) + { + if (events.size() == 0) { + return; + } + + // CCDB initialization + if (fCurrentRun != bcs.begin().runNumber()) { + fCurrentRun = bcs.begin().runNumber(); + } + + fSelMap.clear(); + // Event loop + for (auto& event : events) { + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event); + if (event.has_mcCollision()) { + VarManager::FillEvent(event.mcCollision()); + } + bool decision = false; + if (fEventCut->IsSelected(VarManager::fgValues)) { + decision = true; + } + fSelMap[event.globalIndex()] = decision; + } + + for (auto& [event1, event2] : selfCombinations(*fMixingBinning, fConfigEventOptions.fConfigMixingDepth.value, -1, events, events)) { + VarManager::ResetValues(0, VarManager::kNVars); + if (!fSelMap[event1.globalIndex()] || !fSelMap[event2.globalIndex()]) { + continue; + } + + // save less events if configured + std::vector fSavelessevents = fConfigDileptonHadronOptions.fConfigSavelessevents.value; + if (fSavelessevents[0] > 1 && event1.globalIndex() % fSavelessevents[0] == fSavelessevents[1]) { + continue; + } + + // Get associated tracks for this event + auto groupedAssocs1 = assocs.sliceBy(preslice, event1.globalIndex()); + if (groupedAssocs1.size() < 2) { + continue; // Need at least 2 tracks for pairing + } + auto groupedAssocs2 = assocs.sliceBy(preslice, event2.globalIndex()); + + // Triple loop: track1 (electron) x track2 (electron) x hadron + for (auto& a1 : groupedAssocs1) { + auto t1 = a1.template track_as(); + + uint32_t filter1 = 0; + // Fill track variables + VarManager::FillTrack(t1); + VarManager::FillTrackCollision(t1, event1); + + // Apply electron cuts and fill histograms + int iCut1 = 0; + for (auto cut1 = fTrackCuts.begin(); cut1 != fTrackCuts.end(); cut1++, iCut1++) { + if ((*cut1)->IsSelected(VarManager::fgValues)) { + filter1 |= (static_cast(1) << iCut1); + } + } + + // Check opposite charge with t2 + for (auto& a2 : groupedAssocs1) { + auto t2 = a2.template track_as(); + + // Avoid double counting: use track globalIndex + if (t2.globalIndex() <= t1.globalIndex()) + continue; + + // Fill track variables for t2 (only once per t2) + VarManager::FillTrack(t2); + VarManager::FillTrackCollision(t2, event1); + + // Compute filter2: which cuts t2 passes + uint32_t filter2 = 0; + int iCut2 = 0; + for (auto cut2 = fTrackCuts.begin(); cut2 != fTrackCuts.end(); cut2++, iCut2++) { + if ((*cut2)->IsSelected(VarManager::fgValues)) { + filter2 |= (static_cast(1) << iCut2); + } + } + + // Both tracks must pass at least one common cut + uint32_t twoTrackFilter = filter1 & filter2; + if (!twoTrackFilter) { + continue; + } + // Fill pair variables for cut + VarManager::FillPair(t1, t2, VarManager::fgValues); + float mass = VarManager::fgValues[VarManager::kMass]; + float pt = VarManager::fgValues[VarManager::kPt]; + float rap = VarManager::fgValues[VarManager::kRap]; + // Apply J/psi cuts + if (mass < fConfigPairOptions.fConfigJpsiMassMin.value || mass > fConfigPairOptions.fConfigJpsiMassMax.value || + pt < fConfigPairOptions.fConfigJpsiPtMin.value || pt > fConfigPairOptions.fConfigJpsiPtMax.value || + std::abs(rap) > fConfigPairOptions.fConfigJpsiRapMax.value) { + continue; + } + if (t1.sign() * t2.sign() >= 0) { + continue; // Must be opposite charge + } + // correlate J/psi with hadrons from different events + for (auto& aHadron : groupedAssocs2) { + auto hadron = aHadron.template track_as(); + // Process dilepton-hadron correlation for each common cut + for (size_t iCut = 0; iCut < fTrackCuts.size(); iCut++) { + if (twoTrackFilter & (static_cast(1) << iCut)) { + runDileptonHadron(t1, t2, iCut, hadron, event2, mcParticles); + } + } + } // end hadron loop + } // end track2 loop + } // end track1 loop + } // end event loop + } + + PresliceUnsorted perReducedMcEvent = aod::mcparticle::mcCollisionId; + template + void runEnergyCorrelators(TEvent const& event1, TEvent const& event2, McParticles const& mcTracks) + { + auto groupedMCTracks1 = mcTracks.sliceBy(perReducedMcEvent, event1.mcCollisionId()); + auto groupedMCTracks2 = mcTracks.sliceBy(perReducedMcEvent, event2.mcCollisionId()); + groupedMCTracks1.bindInternalIndicesTo(&mcTracks); + groupedMCTracks2.bindInternalIndicesTo(&mcTracks); + for (auto& t1 : groupedMCTracks1) { + auto t1_raw = groupedMCTracks1.rawIteratorAt(t1.globalIndex()); + for (auto& sig : fGenMCSignals) { + if (sig->CheckSignal(true, t1_raw)) { + if (t1.mcCollisionId() != event1.mcCollisionId()) { // check that the mc track belongs to the same mc collision as the reconstructed event + continue; + } + VarManager::FillTrackMC(groupedMCTracks1, t1_raw); + if (!MixedEvent) { + fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues); + } + } + } + // apply kinematic cuts for signal + if ((t1_raw.pt() < fConfigPairOptions.fConfigJpsiPtMin || t1_raw.pt() > fConfigPairOptions.fConfigJpsiPtMax)) + continue; + if (abs(t1_raw.y()) > fConfigPairOptions.fConfigJpsiRapMax) + continue; + // for the energy correlators + for (auto& t2 : groupedMCTracks2) { + auto t2_raw = groupedMCTracks2.rawIteratorAt(t2.globalIndex()); + if (t2.mcCollisionId() != event2.mcCollisionId()) { // check that the mc track belongs to the same mc collision as the reconstructed event + continue; + } + if (!t2_raw.isPhysicalPrimary()) { + continue; + } + if (t2_raw.has_mothers()) { + auto mother_raw = t2_raw.template mothers_first_as(); + if (mother_raw.globalIndex() == t1_raw.globalIndex()) { + continue; + } + } + if (t2_raw.pt() < fConfigDileptonHadronOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigDileptonHadronOptions.fConfigMCGenHadronEtaAbs.value) { + continue; + } + std::vector fTransRange = fConfigDileptonHadronOptions.fConfigTransRange; + VarManager::FillEnergyCorrelatorsMC(t1_raw, t2_raw, VarManager::fgValues, fTransRange[0], fTransRange[1]); + for (auto& sig : fGenMCSignals) { + if (sig->CheckSignal(true, t1_raw)) { + if (!MixedEvent) { + fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues); + } + if (MixedEvent) { + fHistMan->FillHistClass(Form("MCTruthEenergyCorrelatorsME_%s", sig->GetName()), VarManager::fgValues); + } + } + } + } + } + } + + void processMCGenEnergyCorrelators(soa::Filtered& events, + McCollisions const& /*mcEvents*/, McParticles const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + for (auto& event : events) { + // Fill event variables first + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event); + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + if (!event.has_mcCollision()) { + continue; + } + // saveless events for the energy correlator analysis + std::vector fSavelessevents = fConfigDileptonHadronOptions.fConfigSavelessevents.value; + if (fSavelessevents[0] > 1 && event.globalIndex() % fSavelessevents[0] == fSavelessevents[1]) { + continue; + } + runEnergyCorrelators(event, event, mcTracks); + } + } + + void processMCGenEnergyCorrelatorsME(soa::Filtered& events, + McCollisions const& /*mcEvents*/, McParticles const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + // loop over two event comibnations + for (auto& [event1, event2] : selfCombinations(*fMixingBinning, fConfigEventOptions.fConfigMixingDepth.value, -1, events, events)) { + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event1); + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + VarManager::FillEvent(event2); + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + if (!event1.has_mcCollision() || !event2.has_mcCollision()) { + continue; + } + // saveless events for the energy correlator analysis + std::vector fSavelessevents = fConfigDileptonHadronOptions.fConfigSavelessevents.value; + if (fSavelessevents[0] > 1 && event1.globalIndex() % fSavelessevents[0] == fSavelessevents[1]) { + continue; + } + runEnergyCorrelators(event1, event2, mcTracks); + } + } + + void processMCGenEnergyCorrelatorsPion(soa::Filtered& events, + McCollisions const& /*mcEvents*/, McParticles const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + for (auto& event : events) { + // Fill event variables first + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event); + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + if (!event.has_mcCollision()) { + continue; + } + // saveless events for the energy correlator analysis + std::vector fSavelessevents = fConfigDileptonHadronOptions.fConfigSavelessevents.value; + if (fSavelessevents[0] > 1 && event.globalIndex() % fSavelessevents[0] == fSavelessevents[1]) { + continue; + } + runEnergyCorrelators(event, event, mcTracks); + } + } + + void processMCGenEnergyCorrelatorsPionME(soa::Filtered& events, + McCollisions const& /*mcEvents*/, McParticles const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + // loop over two event comibnations + for (auto& [event1, event2] : selfCombinations(*fMixingBinning, fConfigEventOptions.fConfigMixingDepth.value, -1, events, events)) { + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event1); + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + VarManager::FillEvent(event2); + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + if (!event1.has_mcCollision() || !event2.has_mcCollision()) { + continue; + } + // saveless events for the energy correlator analysis + std::vector fSavelessevents = fConfigDileptonHadronOptions.fConfigSavelessevents.value; + if (fSavelessevents[0] > 1 && event1.globalIndex() % fSavelessevents[0] == fSavelessevents[1]) { + continue; + } + runEnergyCorrelators(event1, event2, mcTracks); + } + } + + void processDummy(aod::Collisions const&) + { + // Do nothing + } + + PROCESS_SWITCH(AnalysisEnergyCorrelator, processBarrel, "Process barrel analysis", false); + PROCESS_SWITCH(AnalysisEnergyCorrelator, processBarrelMixedEvent, "Run barrel dilepton-hadron mixed event pairing", false); + PROCESS_SWITCH(AnalysisEnergyCorrelator, processMCGenEnergyCorrelators, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); + PROCESS_SWITCH(AnalysisEnergyCorrelator, processMCGenEnergyCorrelatorsPion, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); + PROCESS_SWITCH(AnalysisEnergyCorrelator, processMCGenEnergyCorrelatorsME, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); + PROCESS_SWITCH(AnalysisEnergyCorrelator, processMCGenEnergyCorrelatorsPionME, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); + PROCESS_SWITCH(AnalysisEnergyCorrelator, processDummy, "Dummy process function", true); +}; + +// Histogram definitions +void DefineHistograms(HistogramManager* histMan, TString histClasses, const char* histGroups) +{ + 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()); + TString histName = histGroups; + // NOTE: The level of detail for histogramming can be controlled via configurables + if (classStr.Contains("TimeFrameStats")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "timeframe"); + } + if (classStr.Contains("Event")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "event", histName); + } + if ((classStr.Contains("Track") || classStr.Contains("Assoc")) && !classStr.Contains("Pairs")) { + if (classStr.Contains("Barrel")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "track", histName); + } + } + if (classStr.Contains("Pairs")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair", histName); + } + if (classStr.Contains("DileptonTrack")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-track", histName); + } + if (classStr.Contains("MCTruthEenergyCorrelators")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "energy-correlator-gen"); + } + if (classStr.Contains("MCTruthGenSel")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_track"); + } + } // end loop over histogram classes +} + +// Workflow definition +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}