diff --git a/PWGHF/HFC/TableProducer/correlatorFlowCharmHadronsReduced.cxx b/PWGHF/HFC/TableProducer/correlatorFlowCharmHadronsReduced.cxx index 50e99a602e1..9740b182ae1 100644 --- a/PWGHF/HFC/TableProducer/correlatorFlowCharmHadronsReduced.cxx +++ b/PWGHF/HFC/TableProducer/correlatorFlowCharmHadronsReduced.cxx @@ -17,6 +17,7 @@ #include "Common/Core/RecoDecay.h" +#include #include #include #include @@ -33,9 +34,17 @@ #include #include +#include +#include +#include +#include +#include + +#include #include #include #include +#include #include #include @@ -82,20 +91,18 @@ double getPhi(const TTrack& track) } } +enum HistoEffType { + kTH1 = 0, + kTH2, + kTH3, + kTHnSparse +}; + struct HfCorrelatorFlowCharmHadronsReduced { - // Produces rowPairSECharmHads; //! Correlation pairs information Same Event - // Produces rowPairMECharmHads; //! Correlation pairs information Mixed Event - // Produces rowPairSEHadHads; //! Correlation pairs information Same Event - // Produces rowPairMEHadHads; //! Correlation pairs information Mixed Event - // Produces rowCollInfos; //! Collision info - - Configurable fillSparses{"fillSparses", true, "Fill sparse histograms"}; - Configurable fillTables{"fillTables", false, "Fill tables"}; + Configurable numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"}; - Configurable> binsPtTrig{"binsPtTrig", std::vector{0., 3., 5., 8., 16., 36.}, "pT bin limits for trigger candidates"}; Configurable> bdtScore0PtMaxs{"bdtScore0PtMaxs", std::vector{0.1, 0.1, 0.1, 0.1, 0.1}, "pT-differential maximum score 0 for charm candidates"}; Configurable> bdtScore1PtMins{"bdtScore1PtMins", std::vector{0.1, 0.1, 0.1, 0.1, 0.1}, "pT-differential minimum score 1 for charm candidates"}; - Configurable> binsPtAssoc{"binsPtAssoc", std::vector{0.3, 1., 2., 50.}, "pT bin limits for associated particles"}; Configurable selCentrality{"selCentrality", false, "Use centrality to select events"}; Configurable selMultiplicity{"selMultiplicity", false, "Use multiplicity to select events"}; Configurable centMin{"centMin", 0, "min. centrality"}; @@ -112,10 +119,26 @@ struct HfCorrelatorFlowCharmHadronsReduced { Configurable ptMaxForDownSample{"ptMaxForDownSample", 10., "Maximum pt for the application of the downsampling factor"}; Configurable centMultMaxForDownSample{"centMultMaxForDownSample", 10., "Maximum centrality for the application of the downsampling factor"}; + // Efficiency correction + Configurable doEfficiencyCorrection{"doEfficiencyCorrection", false, "Flag for applying efficiency corrections"}; + Configurable effHistTypeCharm{"effHistTypeCharm", 0, "Histogram type for charm efficiency (0=TH1, 1=TH2, 2=TH3, 3=THnSparse)"}; + Configurable effHistTypeAssoc{"effHistTypeAssoc", 0, "Histogram type for associated efficiency (0=TH1, 1=TH2, 2=TH3, 3=THnSparse)"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable assocEffCcdbPath{"assocEffCcdbPath", "", "CCDB path for associated efficiency"}; + Configurable promptEffCcdbPath{"promptEffCcdbPath", "", "CCDB path for trigger efficiency"}; + Configurable timestampCcdb{"timestampCcdb", -1, "timestamp of the efficiency files used to query in CCDB"}; + Configurable ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; + SliceCache cache; + Service ccdb{}; + std::shared_ptr hEffPrompt = nullptr; + std::shared_ptr hEffAssoc = nullptr; + int poolBins{0}; + std::vector mPtEdgesTrig; + using SameEvtPairsChHadWEffInfo = soa::Filtered>; using SameEvtPairsChHad = soa::Filtered>; using SameEvtPairsHadHad = soa::Filtered>; using AssocTracks = soa::Filtered>; @@ -128,9 +151,11 @@ struct HfCorrelatorFlowCharmHadronsReduced { Preslice assocTracksPerCol = aod::hf_correl_charm_had_reduced::hfcRedCorrCollId; Preslice trigCharmCandsPerCol = aod::hf_correl_charm_had_reduced::hfcRedCorrCollId; - ConfigurableAxis zPoolBins{"zPoolBins", {VARIABLE_WIDTH, -10.0, -2.5, 2.5, 10.0}, "Z vertex position pools"}; - ConfigurableAxis multPoolBins{"multPoolBins", {VARIABLE_WIDTH, 0., 900., 1800., 6000.}, "Event multiplicity pools (FT0M)"}; - ConfigurableAxis centPoolBins{"centPoolBins", {VARIABLE_WIDTH, 0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100}, "Event centrality pools"}; + ConfigurableAxis binsPtTrig{"binsPtTrig", {VARIABLE_WIDTH, 0., 1., 2., 3., 4., 5., 6., 8., 10., 12., 16., 24.}, "Trigger candidate pT bins"}; + ConfigurableAxis binsPtAssoc{"binsPtAssoc", {VARIABLE_WIDTH, 0.3, 1., 2., 50.}, "Associated particle pT bins"}; + ConfigurableAxis binsZPool{"binsZPool", {VARIABLE_WIDTH, -10.0, -2.5, 2.5, 10.0}, "Z vertex position pools"}; + ConfigurableAxis binsMultPool{"binsMultPool", {VARIABLE_WIDTH, 0., 900., 1800., 6000.}, "Event multiplicity pools (FT0M)"}; + ConfigurableAxis binsCentPool{"binsCentPool", {VARIABLE_WIDTH, 0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100}, "Event centrality pools"}; ConfigurableAxis binsInvMass{"binsInvMass", {300, 1.6, 2.2}, "Invariant mass bins"}; ConfigurableAxis binsMultFT0M{"binsMultFT0M", {100, 0., 10000.}, "Multiplicity as FT0M signal amplitude"}; ConfigurableAxis binsCent{"binsCent", {100, 0., 100.}, "Centrality bins"}; @@ -142,6 +167,34 @@ struct HfCorrelatorFlowCharmHadronsReduced { HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; + std::shared_ptr loadEfficiencyHistogram(const std::string& path, int histoType) + { + TObject* object = nullptr; + + switch (histoType) { + case kTH1: + object = ccdb->getForTimeStamp(Form("%s/TH1/", path.c_str()), timestampCcdb); + break; + case kTH2: + object = ccdb->getForTimeStamp(Form("%s/TH2/", path.c_str()), timestampCcdb); + break; + case kTH3: + object = ccdb->getForTimeStamp(Form("%s/TH3/", path.c_str()), timestampCcdb); + break; + case kTHnSparse: + object = ccdb->getForTimeStamp(Form("%s/THnSparse/", path.c_str()), timestampCcdb); + break; + default: + LOGF(fatal, "Invalid efficiency histogram type: %d", histoType); + } + + if (!object) { + LOGF(fatal, "Could not load efficiency histogram from %s", path.c_str()); + } + + return std::shared_ptr(object); + } + void init(InitContext&) { if (selCentrality && selMultiplicity) { @@ -151,31 +204,26 @@ struct HfCorrelatorFlowCharmHadronsReduced { LOGP(fatal, "You must select events based on either centrality or multiplicity! Please check your configuration!"); } - // Fatal if selCentrality and process functions using multiplicity mixing, or selMultiplicity and process functions using centrality mixing - if (selCentrality && (doprocessSameEventCharmHadWMultMix || doprocessMixedEventCharmHadWMultMix || doprocessSameEventHadHadWMultMix || doprocessMixedEventHadHadWMultMix)) { - LOGP(fatal, "You cannot select events based on centrality and process functions using multiplicity mixing! Please check your configuration!"); - } - if (selMultiplicity && (doprocessSameEventCharmHadWCentMix || doprocessMixedEventCharmHadWCentMix || doprocessSameEventHadHadWCentMix || doprocessMixedEventHadHadWCentMix)) { - LOGP(fatal, "You cannot select events based on multiplicity and process functions using centrality mixing! Please check your configuration!"); - } - - if ((doprocessSameEventCharmHadWCentMix && doprocessMixedEventCharmHadWMultMix) || - (doprocessSameEventCharmHadWMultMix && doprocessMixedEventCharmHadWCentMix) || - (doprocessSameEventHadHadWCentMix && doprocessMixedEventHadHadWMultMix) || - (doprocessSameEventHadHadWMultMix && doprocessMixedEventHadHadWCentMix)) { + if ((doprocessSECharmHadWCentMix && doprocessMECharmHadWMultMix) || + (doprocessSECharmHadWMultMix && doprocessMECharmHadWCentMix) || + (doprocessSEHadHadWCentMix && doprocessMEHadHadWMultMix) || + (doprocessSEHadHadWMultMix && doprocessMEHadHadWCentMix)) { LOGP(fatal, "You cannot mix centrality and multiplicity mixing in the same processing! Please check your configuration!"); } - if (!fillSparses && !fillTables) { - LOGP(fatal, "At least one of fillSparses or fillTables must be true!"); - } - if ((binsPtTrig.value.size() != (bdtScore0PtMaxs.value.size() + 1) || binsPtTrig.value.size() != (bdtScore1PtMins.value.size() + 1))) { + + mPtEdgesTrig.assign(binsPtTrig.value.begin() + 1, binsPtTrig.value.end()); + if ((mPtEdgesTrig.size() != (bdtScore0PtMaxs.value.size() + 1) || + mPtEdgesTrig.size() != (bdtScore1PtMins.value.size() + 1))) { + LOG(info) << "mPtEdgesTrig.size(): " << mPtEdgesTrig.size(); + LOG(info) << "bdtScore0PtMaxs.value.size(): " << bdtScore0PtMaxs.value.size(); + LOG(info) << "bdtScore1PtMins.value.size(): " << bdtScore1PtMins.value.size(); LOGP(fatal, "The size of bdtScore0PtMaxs and bdtScore1PtMins must be the one of binsPtTrig minus one!"); } - if (doprocessSameEventCharmHadWCentMix || doprocessSameEventHadHadWCentMix || doprocessMixedEventCharmHadWCentMix || doprocessMixedEventHadHadWCentMix || doprocessSameEventCharmHadWCentMixBase || doprocessMixedEventCharmHadWCentMixBase) { - poolBins = (centPoolBins->size() - 2) * (zPoolBins->size() - 2); + if (doprocessSECharmHadWCentMix || doprocessSEHadHadWCentMix || doprocessMECharmHadWCentMix || doprocessMEHadHadWCentMix) { + poolBins = (binsCentPool->size() - 2) * (binsZPool->size() - 2); } else { - poolBins = (multPoolBins->size() - 2) * (zPoolBins->size() - 2); + poolBins = (binsMultPool->size() - 2) * (binsZPool->size() - 2); } const AxisSpec axisInvMass{binsInvMass, "Inv. mass (GeV/#it{c}^{2})"}; @@ -185,19 +233,19 @@ struct HfCorrelatorFlowCharmHadronsReduced { const AxisSpec axisPoolBin = {poolBins, 0., static_cast(poolBins), "PoolBin"}; const AxisSpec axisDeltaEta = {binsDeltaEta, "#Delta#it{#eta}"}; const AxisSpec axisDeltaPhi = {binsDeltaPhi, "#Delta#it{#varphi}"}; - const AxisSpec axisPtTrig = {std::vector(binsPtTrig), "#it{p}_{T} Trig (GeV/#it{c})"}; - const AxisSpec axisPtAssoc = {std::vector(binsPtAssoc), "#it{p}_{T} Assoc (GeV/#it{c})"}; + const AxisSpec axisPtTrig = {binsPtTrig, "#it{p}_{T} Trig (GeV/#it{c})"}; + const AxisSpec axisPtAssoc = {binsPtAssoc, "#it{p}_{T} Assoc (GeV/#it{c})"}; const AxisSpec axisMlOne{binsMlOne, "bdtScore0"}; const AxisSpec axisMlTwo{binsMlTwo, "bdtScore1"}; // Histograms for data analysis - if (doprocessSameEventCharmHadWCentMix || doprocessSameEventHadHadWCentMix) { + if (doprocessSECharmHadWCentMix || doprocessSECharmHadWCentMixWEffInfo || doprocessSEHadHadWCentMix) { registry.add("hCentPoolBinSE", "Centrality SE", {HistType::kTH2F, {{axisCent}, {axisPoolBin}}}); - } else if (doprocessSameEventCharmHadWMultMix || doprocessSameEventHadHadWMultMix) { + } else if (doprocessSECharmHadWMultMix || doprocessSECharmHadWMultMixWEffInfo || doprocessSEHadHadWMultMix) { registry.add("hMultFT0MPoolBinSE", "Multiplicity FT0M SE", {HistType::kTH2F, {{axisMultFT0M}, {axisPoolBin}}}); - } else if (doprocessMixedEventCharmHadWCentMix || doprocessMixedEventHadHadWCentMix) { + } else if (doprocessMECharmHadWCentMix || doprocessMECharmHadWCentMixWEffInfo || doprocessMEHadHadWCentMix) { registry.add("hCentPoolBinME", "Centrality ME", {HistType::kTH2F, {{axisCent}, {axisPoolBin}}}); - } else if (doprocessMixedEventCharmHadWMultMix || doprocessMixedEventHadHadWMultMix) { + } else if (doprocessMECharmHadWMultMix || doprocessMECharmHadWMultMixWEffInfo || doprocessMEHadHadWMultMix) { registry.add("hMultFT0MPoolBinME", "Multiplicity FT0M ME", {HistType::kTH2F, {{axisMultFT0M}, {axisPoolBin}}}); } registry.add("hZVtxPoolBinSE", "z vertex SE", {HistType::kTH2F, {{axisPosZ}, {axisPoolBin}}}); @@ -208,28 +256,116 @@ struct HfCorrelatorFlowCharmHadronsReduced { registry.add("hPoolBinTrigME", "Trigger candidates pool bin ME", {HistType::kTH1F, {axisPoolBin}}); registry.add("hPoolBinAssocSE", "Associated particles pool bin SE", {HistType::kTH1F, {axisPoolBin}}); registry.add("hPoolBinAssocME", "Associated particles pool bin ME", {HistType::kTH1F, {axisPoolBin}}); - if (fillSparses) { - std::vector axesTrigger = {axisInvMass, axisPtTrig, axisMlOne, axisMlTwo}; - std::vector axes = {axisPoolBin, axisPtTrig, axisPtAssoc, axisDeltaEta, axisDeltaPhi}; - if (doprocessSameEventHadHadWCentMix || doprocessSameEventHadHadWMultMix) { - registry.add("hSparseCorrelationsSEHadHad", "THn for SE Had-Had correlations", HistType::kTHnSparseF, axes); - } else if (doprocessMixedEventHadHadWCentMix || doprocessMixedEventHadHadWMultMix) { - registry.add("hSparseCorrelationsMEHadHad", "THn for ME Had-Had correlations", HistType::kTHnSparseF, axes); - } else { - // axes.insert(axes.end(), {axisInvMass}); - axes.insert(axes.end(), {axisInvMass, axisMlOne, axisMlTwo}); - if (doprocessSameEventCharmHadWCentMix || doprocessSameEventCharmHadWMultMix || doprocessSameEventCharmHadWCentMixBase) { - registry.add("hSparseCorrelationsSECharmHad", "THn for SE Charm-Had correlations", HistType::kTHnSparseF, axes); - } else if (doprocessMixedEventCharmHadWCentMix || doprocessMixedEventCharmHadWMultMix || doprocessMixedEventCharmHadWCentMixBase) { - registry.add("hSparseCorrelationsMECharmHad", "THn for ME Charm-Had correlations", HistType::kTHnSparseF, axes); + + std::vector axesTrigger = {axisInvMass, axisPtTrig, axisMlOne, axisMlTwo}; + std::vector axes = {axisPoolBin, axisPtTrig, axisPtAssoc, axisDeltaEta, axisDeltaPhi}; + if (doprocessSEHadHadWCentMix || doprocessSEHadHadWMultMix) { + registry.add("hSparseCorrelationsSEHadHad", "THn for SE Had-Had correlations", HistType::kTHnSparseF, axes); + if (doEfficiencyCorrection) { + registry.get(HIST("hSparseCorrelationsSEHadHad"))->Sumw2(); + } + } else if (doprocessMEHadHadWCentMix || doprocessMEHadHadWMultMix) { + registry.add("hSparseCorrelationsMEHadHad", "THn for ME Had-Had correlations", HistType::kTHnSparseF, axes); + if (doEfficiencyCorrection) { + registry.get(HIST("hSparseCorrelationsMEHadHad"))->Sumw2(); + } + } else { + axes.insert(axes.end(), {axisInvMass, axisMlOne, axisMlTwo}); + if (doprocessSECharmHadWCentMix || doprocessSECharmHadWCentMixWEffInfo || doprocessSECharmHadWMultMix || doprocessSECharmHadWMultMixWEffInfo) { + registry.add("hSparseCorrelationsSECharmHad", "THn for SE Charm-Had correlations", HistType::kTHnSparseF, axes); + if (doEfficiencyCorrection) { + registry.get(HIST("hSparseCorrelationsSECharmHad"))->Sumw2(); + } + } else if (doprocessMECharmHadWCentMix || doprocessMECharmHadWMultMixWEffInfo || doprocessMECharmHadWMultMix || doprocessMECharmHadWCentMixWEffInfo) { + registry.add("hSparseCorrelationsMECharmHad", "THn for ME Charm-Had correlations", HistType::kTHnSparseF, axes); + if (doEfficiencyCorrection) { + registry.get(HIST("hSparseCorrelationsMECharmHad"))->Sumw2(); + } + } + if (doprocessCharmTriggers || doprocessCharmTriggersWEffInfo) { + registry.add("hSparseTrigCandsCharm", "THn for Charm trigger candidates", HistType::kTHnSparseF, axesTrigger); + if (doEfficiencyCorrection) { + registry.get(HIST("hSparseTrigCandsCharm"))->Sumw2(); + } + } + } + + if (doEfficiencyCorrection) { + hEffPrompt = loadEfficiencyHistogram(promptEffCcdbPath, effHistTypeCharm); + hEffAssoc = loadEfficiencyHistogram(assocEffCcdbPath, effHistTypeAssoc); + } + } + + template + double getEfficiency(TObject* hist, int histoType, const TCand& cand) //, const std::vector& values) + { + if (!hist) { + LOG(debug) << "Candidate does not have required efficiency histogram."; + return 0.; + } + + std::vector kineFeatures{}; + + // Charm triggers + if constexpr (IsCharm) { + kineFeatures = {cand.ptTrig(), cand.bdtScore0Trig(), cand.bdtScore1Trig(), cand.etaTrig(), cand.phiTrig()}; + } else { + auto collision = cand.template hfcRedCorrColl_as(); + kineFeatures = {cand.ptAssoc(), cand.etaAssoc(), cand.phiAssoc(), collision.posZ()}; + } + + switch (histoType) { + + case kTH1: { + auto* h = dynamic_cast(hist); + if (!h || kineFeatures.size() < 1) { + return 0.; + } + return h->GetBinContent(h->FindBin(kineFeatures[0])); + } + + case kTH2: { + auto* h = dynamic_cast(hist); + if (!h || kineFeatures.size() < 2) { + return 0.; + } + return h->GetBinContent(h->FindBin(kineFeatures[0], kineFeatures[1])); + } + + case kTH3: { + auto* h = dynamic_cast(hist); + if (!h || kineFeatures.size() < 3) { + return 0.; } - if (doprocessCharmTriggers) { - registry.add("hSparseTrigCandsCharm", "THn for Charm trigger candidates", HistType::kTHnSparseF, axesTrigger); + return h->GetBinContent( + h->FindBin(kineFeatures[0], kineFeatures[1], kineFeatures[2])); + } + + case kTHnSparse: { + auto* h = dynamic_cast(hist); + if (!h || kineFeatures.size() < static_cast(h->GetNdimensions())) { + return 0.; } + + std::vector bins(h->GetNdimensions()); + + for (int i = 0; i < h->GetNdimensions(); ++i) { + bins[i] = h->GetAxis(i)->FindBin(kineFeatures[i]); + } + + return h->GetBinContent(bins.data()); } + + default: + return 0.; } } + bool isValidEfficiency(double efficiency) + { + return std::isfinite(efficiency) && efficiency > 0.; + } + template bool isSelCollision(const TColl& collision) { @@ -277,8 +413,8 @@ struct HfCorrelatorFlowCharmHadronsReduced { bool isSelBdtScoreCut(TCand const& cand, double ptTrig) { - for (size_t iPt = 0; iPt < binsPtTrig.value.size() - 1; iPt++) { - if (ptTrig >= binsPtTrig.value[iPt] && ptTrig < binsPtTrig.value[iPt + 1]) { + for (size_t iPt = 0; iPt < mPtEdgesTrig.size() - 1; iPt++) { + if (ptTrig >= mPtEdgesTrig[iPt] && ptTrig < mPtEdgesTrig[iPt + 1]) { return (cand.bdtScore0Trig() < bdtScore0PtMaxs.value[iPt]) && (cand.bdtScore1Trig() > bdtScore1PtMins.value[iPt]); } } @@ -290,7 +426,7 @@ struct HfCorrelatorFlowCharmHadronsReduced { /// \param trigCands are the selected trigger candidates /// \param assocTracks are the selected associated tracks /// \param binPolicy is the binning policy for the correlation - template + template void fillSameEvent(TPair const& pair, TTrigCand const& trigCand, const TBinningType& binPolicy) @@ -318,33 +454,61 @@ struct HfCorrelatorFlowCharmHadronsReduced { int const poolBin = getPoolBin(collision, binPolicy); registry.fill(HIST("hPoolBinTrigSE"), poolBin); registry.fill(HIST("hPoolBinAssocSE"), poolBin); - // if constexpr (FillTables) { - // if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - // rowPairSECharmHads(poolBin, ptTrig, pair.ptAssoc(), pair.deltaEta(), pair.deltaPhi(), - // trigCand.invMassTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig(), - // pair.nTpcCrossedRowsAssoc(), pair.itsClsMapAssoc(), pair.itsNClsAssoc(), pair.dcaXYAssoc(), pair.dcaZAssoc()); - // } else { - // rowPairSEHadHads(poolBin, ptTrig, pair.ptAssoc(), pair.deltaEta(), pair.deltaPhi(), - // trigCand.nTpcCrossedRowsTrig(), trigCand.itsClsMapTrig(), trigCand.itsNClsTrig(), trigCand.dcaXYTrig(), trigCand.dcaZTrig(), - // pair.nTpcCrossedRowsAssoc(), pair.itsClsMapAssoc(), pair.itsNClsAssoc(), pair.dcaXYAssoc(), pair.dcaZAssoc()); - // } - // rowCollInfos(collision.multiplicity(), collision.numPvContrib(), collision.centrality()); - // } - if constexpr (FillSparses) { - if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - registry.fill(HIST("hSparseCorrelationsSECharmHad"), poolBin, ptTrig, pair.ptAssoc(), pair.deltaEta(), - pair.deltaPhi(), trigCand.invMassTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()); - } else { - registry.fill(HIST("hSparseCorrelationsSEHadHad"), poolBin, ptTrig, pair.ptAssoc(), pair.deltaEta(), pair.deltaPhi()); + + double effAssoc{1.f}; + if constexpr (CorrectEff) { + if (doEfficiencyCorrection) { + effAssoc = getEfficiency(hEffAssoc.get(), effHistTypeAssoc.value, pair); + if (!isValidEfficiency(effAssoc)) { + LOGF(warn, "Invalid efficiency value for associated track: effAssoc = %f", effAssoc); + return; + } } } + + if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases + + double effTrig{1.f}; + if constexpr (CorrectEff) { + if (doEfficiencyCorrection) { + effTrig = getEfficiency(hEffPrompt.get(), effHistTypeCharm.value, trigCand); + if (!isValidEfficiency(effTrig)) { + LOGF(warn, "Invalid efficiency value for trigger candidate: effTrig = %f", effTrig); + return; + } + } + } + + double weight = doEfficiencyCorrection ? 1 / (effTrig * effAssoc) : 1.f; + + std::vector values = {static_cast(poolBin), ptTrig, pair.ptAssoc(), pair.deltaEta(), pair.deltaPhi(), trigCand.invMassTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()}; + + registry.get(HIST("hSparseCorrelationsSECharmHad"))->Fill(values.data(), weight); + } else { + + double effTrig{1.f}; + if constexpr (CorrectEff) { + if (doEfficiencyCorrection) { + effTrig = getEfficiency(hEffAssoc.get(), effHistTypeAssoc.value, trigCand); + if (!isValidEfficiency(effTrig)) { + LOGF(warn, "Invalid efficiency value for trigger candidate: effTrig = %f", effTrig); + return; + } + } + } + + double weight = doEfficiencyCorrection ? 1 / (effTrig * effAssoc) : 1.f; + + std::vector valuesHadHad = {static_cast(poolBin), ptTrig, pair.ptAssoc(), pair.deltaEta(), pair.deltaPhi()}; + registry.get(HIST("hSparseCorrelationsSEHadHad"))->Fill(valuesHadHad.data(), weight); + } } /// Save info for Mixed Event pairs /// \param collisions are the selected collisions /// \param pairs are the mixed event pairs of trigger candidates and associated tracks /// \param binPolicy is the binning policy for the correlation - template + template void fillMixedEvent(TPairs const& pairs, const TBinningType& binPolicy) { @@ -392,248 +556,174 @@ struct HfCorrelatorFlowCharmHadronsReduced { } } double const deltaPhi = RecoDecay::constrainAngle(getPhi(assocTrack) - getPhi(trigCand), -o2::constants::math::PIHalf); - // if constexpr (FillTables) { - // if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - // rowPairMECharmHads(poolBinTrig, ptTrig, ptAssoc, deltaEta, deltaPhi, - // trigCand.invMassTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig(), - // assocTrack.nTpcCrossedRowsAssoc(), assocTrack.itsClsMapAssoc(), assocTrack.itsNClsAssoc(), assocTrack.dcaXYAssoc(), assocTrack.dcaZAssoc()); - // } else { - // rowPairMEHadHads(poolBinTrig, ptTrig, ptAssoc, deltaEta, deltaPhi, - // trigCand.nTpcCrossedRowsAssoc(), trigCand.itsClsMapAssoc(), trigCand.itsNClsAssoc(), trigCand.dcaXYAssoc(), trigCand.dcaZAssoc(), - // assocTrack.nTpcCrossedRowsAssoc(), assocTrack.itsClsMapAssoc(), assocTrack.itsNClsAssoc(), assocTrack.dcaXYAssoc(), assocTrack.dcaZAssoc()); - // } - // rowCollInfos(trigColl.multiplicity(), trigColl.numPvContrib(), trigColl.centrality()); - // } - if constexpr (FillSparses) { - if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - registry.fill(HIST("hSparseCorrelationsMECharmHad"), poolBinTrig, ptTrig, ptAssoc, deltaEta, - deltaPhi, trigCand.invMassTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()); - } else { - registry.fill(HIST("hSparseCorrelationsMEHadHad"), poolBinTrig, ptTrig, ptAssoc, deltaEta, deltaPhi); - } - } - } - } - } - /// Correlations for Same Event pairs - /// \param poolBin collision pool bin based on multiplicity and z-vertex position - /// \param trigCandsThisColl are the selected trigger candidates in the collision - /// \param assocTracksThisColl are the selected associated tracks in the collision - template - void doCorrelationsSameEvent(int poolBin, - const TTrigCand& trigCandsThisColl, - const TTrackAssoc& assocTracksThisColl) - { - for (const auto& trigCand : trigCandsThisColl) { - double const ptTrig = trigCand.ptTrig(); - if constexpr (requires { trigCand.bdtScore0Trig(); }) { // ML selection on bkg score for Charm-Had case - if (!isSelBdtScoreCut(trigCand, ptTrig)) { - continue; - } - } - - for (const auto& assTrk : assocTracksThisColl) { - // TODO: Remove Ds daughters - /*if (assTrk.originTrackId() == candidate.prong0Id() || - assTrk.originTrackId() == candidate.prong1Id() || - assTrk.originTrackId() == candidate.prong2Id()) { - continue; - }*/ - // TODO: DCA cut - double deltaPhi = RecoDecay::constrainAngle(assTrk.phiAssoc() - trigCand.phiTrig(), -o2::constants::math::PIHalf); - double deltaEta = assTrk.etaAssoc() - trigCand.etaTrig(); - if constexpr (FillSparses) { - if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - registry.fill(HIST("hSparseCorrelationsSECharmHad"), poolBin, ptTrig, assTrk.ptAssoc(), deltaEta, - deltaPhi, trigCand.invMassTrig()); - } else { - registry.fill(HIST("hSparseCorrelationsSEHadHad"), poolBin, ptTrig, assTrk.ptAssoc(), deltaEta, deltaPhi); + double effAssoc{1.f}; + if (doEfficiencyCorrection) { + effAssoc = getEfficiency(hEffAssoc.get(), effHistTypeAssoc.value, assocTrack); + if (!isValidEfficiency(effAssoc)) { + LOGF(warn, "Invalid efficiency value for associated track: effAssoc = %f", effAssoc); + continue; } } - } - } - } - /// Correlations for Mixed Event pairs - /// \param collisions are the selected collisions - /// \param trigCands are the trigger candidates - /// \param assocTracks are the associated tracks - /// \param binPolicy is the binning policy for the correlation - template - void doCorrelationsMixedEvent(const TCollision& collisions, - const TTrigCand& trigCands, - const TTrackAssoc& assocTracks, - const TBinningType& binPolicy) - { - auto tracksTuple = std::make_tuple(trigCands, assocTracks); + if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - Pair pairData{binPolicy, numberEventsMixed, -1, collisions, tracksTuple, &cache}; + double effTrig{1.f}; + if (doEfficiencyCorrection) { + effTrig = getEfficiency(hEffPrompt.get(), effHistTypeCharm.value, trigCand); + if (!isValidEfficiency(effTrig)) { + LOGF(warn, "Invalid efficiency value for trigger candidate: effTrig = %f", effTrig); + continue; + } + } - for (const auto& [c1, tracks1, c2, tracks2] : pairData) { - if (!isSelCollision(c1) || !isSelCollision(c2)) { - continue; - } - if (tracks1.size() == 0) { - continue; - } + double weight = doEfficiencyCorrection ? 1 / (effTrig * effAssoc) : 1.f; - int poolBin = binPolicy.getBin({c2.posZ(), selCentrality ? c2.centrality() : c2.multiplicity()}); - int poolBinTrigCand = binPolicy.getBin({c1.posZ(), selCentrality ? c1.centrality() : c1.multiplicity()}); + std::vector values = {static_cast(poolBinTrig), ptTrig, ptAssoc, deltaEta, deltaPhi, trigCand.invMassTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()}; + registry.get(HIST("hSparseCorrelationsMECharmHad"))->Fill(values.data(), weight); + } else { - if (poolBin != poolBinTrigCand) { - LOGF(info, "Error, poolBins are different"); - continue; - } + double effTrig{1.f}; + if (doEfficiencyCorrection) { + effTrig = getEfficiency(hEffAssoc.get(), effHistTypeAssoc.value, trigCand); + if (!isValidEfficiency(effTrig)) { + LOGF(warn, "Invalid efficiency value for trigger candidate: effTrig = %f", effTrig); + continue; + } + } - for (const auto& [trigCand, assTrk] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (!isSelBdtScoreCut(trigCand, trigCand.ptTrig())) { - continue; - } + double weight = doEfficiencyCorrection ? 1 / (effTrig * effAssoc) : 1.f; - double deltaPhi = RecoDecay::constrainAngle(assTrk.phiAssoc() - trigCand.phiTrig(), -o2::constants::math::PIHalf); - double deltaEta = assTrk.etaAssoc() - trigCand.etaTrig(); - if constexpr (FillSparses) { - if constexpr (requires { trigCand.bdtScore0Trig(); }) { // Separate Charm-Had and Had-Had cases - registry.fill(HIST("hSparseCorrelationsMECharmHad"), poolBin, trigCand.ptTrig(), assTrk.ptAssoc(), deltaEta, - deltaPhi, trigCand.invMassTrig()); - } else { - registry.fill(HIST("hSparseCorrelationsMEHadHad"), poolBin, trigCand.ptTrig(), assTrk.ptAssoc(), deltaEta, deltaPhi); - } + std::vector valuesHadHad = {static_cast(poolBinTrig), ptTrig, ptAssoc, deltaEta, deltaPhi}; + registry.get(HIST("hSparseCorrelationsMEHadHad"))->Fill(valuesHadHad.data(), weight); } } } } - void processSameEventCharmHadWMultMix(SameEvtPairsChHad::iterator const& pair, - aod::HfcRedTrigCharms const&, - aod::HfcRedCorrColls const&) + void processSECharmHadWMultMix(SameEvtPairsChHad::iterator const& pair, + aod::HfcRedTrigCharms const&, + aod::HfcRedCorrColls const&) { - BinningMultPosZ binPolicyPosZMult{{zPoolBins, multPoolBins}}; + BinningMultPosZ binPolicyPosZMult{{binsZPool, binsMultPool}}; auto trigCand = pair.template hfcRedTrigCharm_as(); - if (fillSparses && fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZMult); - } else if (fillSparses) { - fillSameEvent(pair, trigCand, binPolicyPosZMult); - } else if (fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZMult); - } + fillSameEvent(pair, trigCand, binPolicyPosZMult); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSameEventCharmHadWMultMix, "Process Same Event for Charm-Had with multiplicity pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSECharmHadWMultMix, "Process Same Event for Charm-Had with multiplicity pools", false); - void processSameEventHadHadWMultMix(SameEvtPairsHadHad::iterator const& pair, - aod::HfcRedTrigTracks const&, - aod::HfcRedCorrColls const&) + void processSECharmHadWMultMixWEffInfo(SameEvtPairsChHadWEffInfo::iterator const& pair, + TrigCharmCands const&, + aod::HfcRedCorrColls const&) { - BinningMultPosZ binPolicyPosZMult{{zPoolBins, multPoolBins}}; + BinningMultPosZ binPolicyPosZMult{{binsZPool, binsMultPool}}; + auto trigCand = pair.template hfcRedTrigCharm_as(); + fillSameEvent(pair, trigCand, binPolicyPosZMult); + } + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSECharmHadWMultMixWEffInfo, "Process Same Event for Charm-Had with multiplicity pools", false); + + void processSEHadHadWMultMix(SameEvtPairsHadHad::iterator const& pair, + aod::HfcRedTrigTracks const&, + aod::HfcRedCorrColls const&) + { + BinningMultPosZ binPolicyPosZMult{{binsZPool, binsMultPool}}; auto trigCand = pair.template hfcRedTrigTrack_as(); - if (fillSparses && fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZMult); - } else if (fillSparses) { - fillSameEvent(pair, trigCand, binPolicyPosZMult); - } else if (fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZMult); - } + fillSameEvent(pair, trigCand, binPolicyPosZMult); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSameEventHadHadWMultMix, "Process Same Event for Had-Had with multiplicity pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSEHadHadWMultMix, "Process Same Event for Had-Had with multiplicity pools", false); - void processSameEventCharmHadWCentMix(SameEvtPairsChHad::iterator const& pair, - aod::HfcRedTrigCharms const&, - aod::HfcRedCorrColls const&) + void processSECharmHadWCentMix(SameEvtPairsChHad::iterator const& pair, + aod::HfcRedTrigCharms const&, + aod::HfcRedCorrColls const&) { - BinningCentPosZ binPolicyPosZCent{{zPoolBins, centPoolBins}}; + BinningCentPosZ binPolicyPosZCent{{binsZPool, binsCentPool}}; auto trigCand = pair.template hfcRedTrigCharm_as(); - if (fillSparses && fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZCent); - } else if (fillSparses) { - fillSameEvent(pair, trigCand, binPolicyPosZCent); - } else if (fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZCent); - } + fillSameEvent(pair, trigCand, binPolicyPosZCent); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSameEventCharmHadWCentMix, "Process Same Event for Charm-Had with centrality pools", true); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSECharmHadWCentMix, "Process Same Event for Charm-Had with centrality pools", true); - void processSameEventHadHadWCentMix(SameEvtPairsHadHad::iterator const& pair, - aod::HfcRedTrigTracks const&, - aod::HfcRedCorrColls const&) + void processSECharmHadWCentMixWEffInfo(SameEvtPairsChHadWEffInfo::iterator const& pair, + TrigCharmCands const&, + aod::HfcRedCorrColls const&) { - BinningCentPosZ binPolicyPosZCent{{zPoolBins, centPoolBins}}; + BinningCentPosZ binPolicyPosZCent{{binsZPool, binsCentPool}}; + auto trigCand = pair.template hfcRedTrigCharm_as(); + fillSameEvent(pair, trigCand, binPolicyPosZCent); + } + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSECharmHadWCentMixWEffInfo, "Process Same Event for Charm-Had with centrality pools", false); + + void processSEHadHadWCentMix(SameEvtPairsHadHad::iterator const& pair, + aod::HfcRedTrigTracks const&, + aod::HfcRedCorrColls const&) + { + BinningCentPosZ binPolicyPosZCent{{binsZPool, binsCentPool}}; auto trigCand = pair.template hfcRedTrigTrack_as(); - if (fillSparses && fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZCent); - } else if (fillSparses) { - fillSameEvent(pair, trigCand, binPolicyPosZCent); - } else if (fillTables) { - fillSameEvent(pair, trigCand, binPolicyPosZCent); - } + fillSameEvent(pair, trigCand, binPolicyPosZCent); + } + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSEHadHadWCentMix, "Process Same Event for Had-Had with centrality pools", false); + + void processMECharmHadWCentMix(aod::HfcRedCorrColls const& collisions, + TrigCharmCands const& candidates, + AssocTracks const& tracks) + { + BinningCentPosZ binPolicyPosZCent{{binsZPool, binsCentPool}}; + auto pairsTuple = std::make_tuple(candidates, tracks); + Pair const pairs{binPolicyPosZCent, numberEventsMixed, -1, collisions, pairsTuple, &cache}; + fillMixedEvent(pairs, binPolicyPosZCent); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSameEventHadHadWCentMix, "Process Same Event for Had-Had with centrality pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMECharmHadWCentMix, "Process Mixed Event for Charm-Had with centrality pools", false); - void processMixedEventCharmHadWCentMix(aod::HfcRedCorrColls const& collisions, + void processMECharmHadWCentMixWEffInfo(aod::HfcRedCorrColls const& collisions, TrigCharmCands const& candidates, AssocTracks const& tracks) { - BinningCentPosZ binPolicyPosZCent{{zPoolBins, centPoolBins}}; + BinningCentPosZ binPolicyPosZCent{{binsZPool, binsCentPool}}; auto pairsTuple = std::make_tuple(candidates, tracks); Pair const pairs{binPolicyPosZCent, numberEventsMixed, -1, collisions, pairsTuple, &cache}; - if (fillSparses && fillTables) { - fillMixedEvent(pairs, binPolicyPosZCent); - } else if (fillSparses) { - fillMixedEvent(pairs, binPolicyPosZCent); - } else if (fillTables) { - fillMixedEvent(pairs, binPolicyPosZCent); - } + fillMixedEvent(pairs, binPolicyPosZCent); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMixedEventCharmHadWCentMix, "Process Mixed Event for Charm-Had with centrality pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMECharmHadWCentMixWEffInfo, "Process Mixed Event for Charm-Had with centrality pools and efficiency information", false); - void processMixedEventCharmHadWMultMix(aod::HfcRedCorrColls const& collisions, + void processMECharmHadWMultMix(aod::HfcRedCorrColls const& collisions, + TrigCharmCands const& candidates, + AssocTracks const& tracks) + { + BinningMultPosZ binPolicyPosZMult{{binsZPool, binsMultPool}}; + auto pairsTuple = std::make_tuple(candidates, tracks); + Pair const pairs{binPolicyPosZMult, numberEventsMixed, -1, collisions, pairsTuple, &cache}; + fillMixedEvent(pairs, binPolicyPosZMult); + } + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMECharmHadWMultMix, "Process Mixed Event for Charm-Had with multiplicity pools", false); + + void processMECharmHadWMultMixWEffInfo(aod::HfcRedCorrColls const& collisions, TrigCharmCands const& candidates, AssocTracks const& tracks) { - BinningMultPosZ binPolicyPosZMult{{zPoolBins, multPoolBins}}; + BinningMultPosZ binPolicyPosZMult{{binsZPool, binsMultPool}}; auto pairsTuple = std::make_tuple(candidates, tracks); Pair const pairs{binPolicyPosZMult, numberEventsMixed, -1, collisions, pairsTuple, &cache}; - if (fillSparses && fillTables) { - fillMixedEvent(pairs, binPolicyPosZMult); - } else if (fillSparses) { - fillMixedEvent(pairs, binPolicyPosZMult); - } else if (fillTables) { - fillMixedEvent(pairs, binPolicyPosZMult); - } + fillMixedEvent(pairs, binPolicyPosZMult); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMixedEventCharmHadWMultMix, "Process Mixed Event for Charm-Had with multiplicity pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMECharmHadWMultMixWEffInfo, "Process Mixed Event for Charm-Had with multiplicity pools and efficiency information", false); - void processMixedEventHadHadWCentMix(aod::HfcRedCorrColls const& collisions, - AssocTracks const& tracks) + void processMEHadHadWCentMix(aod::HfcRedCorrColls const& collisions, + AssocTracks const& tracks) { - BinningCentPosZ binPolicyPosZCent{{zPoolBins, centPoolBins}}; + BinningCentPosZ binPolicyPosZCent{{binsZPool, binsCentPool}}; auto tracksTuple = std::make_tuple(tracks); SameKindPair const pairs{binPolicyPosZCent, numberEventsMixed, -1, collisions, tracksTuple, &cache}; - if (fillSparses && fillTables) { - fillMixedEvent(pairs, binPolicyPosZCent); - } else if (fillSparses) { - fillMixedEvent(pairs, binPolicyPosZCent); - } else if (fillTables) { - fillMixedEvent(pairs, binPolicyPosZCent); - } + fillMixedEvent(pairs, binPolicyPosZCent); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMixedEventHadHadWCentMix, "Process Mixed Event for Had-Had with centrality pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMEHadHadWCentMix, "Process Mixed Event for Had-Had with centrality pools", false); - void processMixedEventHadHadWMultMix(aod::HfcRedCorrColls const& collisions, - AssocTracks const& tracks) + void processMEHadHadWMultMix(aod::HfcRedCorrColls const& collisions, + AssocTracks const& tracks) { - BinningMultPosZ binPolicyPosZMult{{zPoolBins, multPoolBins}}; + BinningMultPosZ binPolicyPosZMult{{binsZPool, binsMultPool}}; auto tracksTuple = std::make_tuple(tracks); SameKindPair const pairs{binPolicyPosZMult, numberEventsMixed, -1, collisions, tracksTuple, &cache}; - if (fillSparses && fillTables) { - fillMixedEvent(pairs, binPolicyPosZMult); - } else if (fillSparses) { - fillMixedEvent(pairs, binPolicyPosZMult); - } else if (fillTables) { - fillMixedEvent(pairs, binPolicyPosZMult); - } + fillMixedEvent(pairs, binPolicyPosZMult); } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMixedEventHadHadWMultMix, "Process Mixed Event for Had-Had with multiplicity pools", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMEHadHadWMultMix, "Process Mixed Event for Had-Had with multiplicity pools", false); void processCharmTriggers(aod::HfcRedTrigCharms const& trigCands, aod::HfcRedCorrColls const&) @@ -647,41 +737,37 @@ struct HfCorrelatorFlowCharmHadronsReduced { if (!isSelBdtScoreCut(trigCand, trigCand.ptTrig())) { continue; } - registry.fill(HIST("hSparseTrigCandsCharm"), trigCand.invMassTrig(), trigCand.ptTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()); + + std::vector values = {trigCand.invMassTrig(), trigCand.ptTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()}; + registry.get(HIST("hSparseTrigCandsCharm"))->Fill(values.data()); } } PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processCharmTriggers, "Process charm trigger info", false); - void processSameEventCharmHadWCentMixBase(aod::HfcRedCorrColls const& collisions, - TrigCharmCands const& candidates, - aod::HfcRedAssBases const& tracks) + void processCharmTriggersWEffInfo(TrigCharmCands const& trigCands, + aod::HfcRedCorrColls const&) { - BinningCentPosZ binPolicyPosZCent{{zPoolBins, centPoolBins}}; - for (const auto& collision : collisions) { + for (const auto& trigCand : trigCands) { + auto collision = trigCand.template hfcRedCorrColl_as(); if (!isSelCollision(collision)) { continue; } - int poolBin = binPolicyPosZCent.getBin({collision.posZ(), collision.centrality()}); - - auto thisCollId = collision.globalIndex(); - auto candsThisColl = candidates.sliceBy(trigCharmCandsPerCol, thisCollId); - auto tracksThisColl = tracks.sliceBy(assocTracksPerCol, thisCollId); + if (!isSelBdtScoreCut(trigCand, trigCand.ptTrig())) { + continue; + } - doCorrelationsSameEvent(poolBin, candsThisColl, tracksThisColl); + double effTrig = doEfficiencyCorrection ? getEfficiency(hEffPrompt.get(), effHistTypeCharm.value, trigCand) : 1.; + if (doEfficiencyCorrection && !isValidEfficiency(effTrig)) { + LOGF(warn, "Invalid efficiency value for trigger candidate: effTrig = %f", effTrig); + continue; + } + std::vector values = {trigCand.invMassTrig(), trigCand.ptTrig(), trigCand.bdtScore0Trig(), trigCand.bdtScore1Trig()}; + double weight = doEfficiencyCorrection ? 1. / effTrig : 1.; + registry.get(HIST("hSparseTrigCandsCharm"))->Fill(values.data(), weight); } } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processSameEventCharmHadWCentMixBase, "Process Same Event base", false); - - void processMixedEventCharmHadWCentMixBase(aod::HfcRedCorrColls const& collisions, - TrigCharmCands const& candidates, - aod::HfcRedAssBases const& tracks) - { - BinningCentPosZ binPolicyPosZCent{{zPoolBins, centPoolBins}}; - - doCorrelationsMixedEvent(collisions, candidates, tracks, binPolicyPosZCent); - } - PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processMixedEventCharmHadWCentMixBase, "Process Mixed Event base", false); + PROCESS_SWITCH(HfCorrelatorFlowCharmHadronsReduced, processCharmTriggersWEffInfo, "Process charm trigger info with efficiency", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx b/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx index 06ae7c8a758..3efa36a6d9f 100644 --- a/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx +++ b/PWGHF/HFC/TableProducer/derivedDataCreatorCorrelationsReduced.cxx @@ -16,11 +16,13 @@ /// \author Wu Chuntai , CCNU, INFN Padova, and Padova University #include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/Core/DecayChannels.h" #include "PWGHF/Core/HfHelper.h" #include "PWGHF/DataModel/AliasTables.h" #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" #include "PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h" +#include "PWGHF/HFC/Utils/utilsFlow.h" #include "PWGHF/Utils/utilsEvSelHf.h" #include "Common/Core/RecoDecay.h" @@ -49,6 +51,7 @@ #include #include #include +#include #include #include @@ -57,19 +60,10 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::hf_centrality; using namespace o2::hf_evsel; +using namespace o2::hf_decay; +using namespace o2::hfc_utils_flow; -enum CandidateType { - DplusToPiKPi = 0, - DsToKKPi, - DsToPiKK, - D0ToPiK, - D0ToKPi, - LcToPKPi, - LcToPiKP, - Hadron -}; - -/// Code to select collisions with at least one Ds meson +/// Code to produce self-contained derived datasets for v2 analyses of charm hadrons via 2PC struct HfDerivedDataCreatorCorrelationsReduced { Produces rowCollisions; // Table with reduced collision info Produces rowSECharmHadPairs; // Table with same-event pairs info @@ -83,6 +77,7 @@ struct HfDerivedDataCreatorCorrelationsReduced { Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4)"}; Configurable selectionFlag{"selectionFlag", 15, "Selection Flag for hadron (ML score tables are required to run the task)"}; Configurable forceCharmInCollision{"forceCharmInCollision", true, "Flag to force charm in collision"}; + Configurable acceptOnlyPrompt{"acceptOnlyPrompt", true, "Flag to accept only prompt charm in the MC process functions"}; Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."}; Configurable yCandMax{"yCandMax", 0.8, "max. cand. rapidity"}; @@ -102,84 +97,95 @@ struct HfDerivedDataCreatorCorrelationsReduced { Configurable downSampleTrksFactor{"downSampleTrksFactor", 1., "Fraction of associated tracks to keep"}; Configurable ptMaxForDownSample{"ptMaxForDownSample", 10., "Maximum pt for the application of the downsampling factor"}; Configurable centMaxForDownSample{"centMaxForDownSample", 101., "Maximum centrality for the application of the downsampling factor"}; - Configurable> binsPtTrig{"binsPtTrig", std::vector{0., 1., 2., 3., 5., 8., 12., 24., 36.}, "pT bin limits for trigger candidates"}; - Configurable> binsPtAssoc{"binsPtAssoc", std::vector{0.2, 1., 2., 50.}, "pT bin limits for associated particles"}; HfEventSelection hfEvSel; // event selection and monitoring + HfcFlowHelper hfcFlowHelper; // flow utilities o2::framework::Service ccdb{}; SliceCache cache; - double massCharm{0.}; - TF1* funcDcaXYPtCutPrimTrk = nullptr; + std::unique_ptr funcDcaXYPtCutPrimTrk{nullptr}; - using CollsWithCentMult = soa::Join; - using CandDsData = soa::Filtered>; - using CandDplusData = soa::Filtered>; + using DataCollsWithCentMult = soa::Join; + using McCollsWithCentMult = soa::Join; + + // Data using CandD0Data = soa::Filtered>; + using CandDplusData = soa::Filtered>; + using CandDsData = soa::Filtered>; using CandLcData = soa::Filtered>; - using TracksData = soa::Filtered>; + using TracksData = soa::Filtered>; //, aod::pidTPCFullMu, aod::pidTPCFullEl, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr>>; + // MC + using CandD0Mc = soa::Filtered>; + using CandDplusMc = soa::Filtered>; + using CandDsMc = soa::Filtered>; + using CandLcMc = soa::Filtered>; + using TracksMc = soa::Filtered>; //, aod::pidTPCFullMu, aod::pidTPCFullEl, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr>>; - Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; - Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; Filter filterSelectD0Candidates = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; + Filter filterSelectDplusCandidates = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; + Filter filterSelectDsCandidates = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; Filter filterSelectLcCandidates = aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlag || aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlag; Filter filterSelectTrkData = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < dcaXYTrkMax) && (nabs(aod::track::dcaZ) < dcaZTrkMax); - Preslice candsDsPerColl = aod::hf_cand::collisionId; - Preslice candsDplusPerColl = aod::hf_cand::collisionId; - Preslice candsD0PerColl = aod::hf_cand::collisionId; - Preslice candsLcPerColl = aod::hf_cand::collisionId; - Preslice trackIndicesPerColl = aod::track::collisionId; - - Partition selectedDsToKKPi = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag; - Partition selectedDsToPiKK = aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; - Partition selectedD0ToPiK = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; - Partition selectedD0ToKPi = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; - Partition selectedLcToPKPi = aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlag; - Partition selectedLcToPiKP = aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlag; - - ConfigurableAxis binsInvMass{"binsInvMass", {300, 1.6, 2.2}, ""}; + Partition selD0ToPiK = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag; + Partition selD0ToKPi = aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; + Partition selDsToKKPi = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag; + Partition selDsToPiKK = aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; + Partition selLcToPKPi = aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlag; + Partition selLcToPiKP = aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlag; + + Partition recoSelD0ToPiK = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_2prong::DecayChannelMain::D0ToPiK)) && aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + Partition recoSelD0ToKPi = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_2prong::DecayChannelMain::D0ToPiK)) && aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + Partition recoSelDPlusToPiKPi = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_3prong::DecayChannelMain::DplusToPiKPi)) && aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + Partition recoSelDsToKKPi = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_3prong::DecayChannelMain::DsToPiKK) && nabs(aod::hf_cand_mc_flag::flagMcDecayChanRec) == static_cast(hf_cand_3prong::DecayChannelResonant::DsToPhiPi)) && aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + Partition recoSelDsToPiKK = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_3prong::DecayChannelMain::DsToPiKK) && nabs(aod::hf_cand_mc_flag::flagMcDecayChanRec) == static_cast(hf_cand_3prong::DecayChannelResonant::DsToPhiPi)) && aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + Partition recoSelLcToPKPi = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_3prong::DecayChannelMain::LcToPKPi)) && aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + Partition recoSelLcToPiKP = (nabs(aod::hf_cand_mc_flag::flagMcMatchRec) == static_cast(hf_cand_3prong::DecayChannelMain::LcToPKPi)) && aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlag && ((!acceptOnlyPrompt) || (acceptOnlyPrompt && aod::hf_cand_mc_flag::originMcRec == static_cast(RecoDecay::OriginType::Prompt))); + + ConfigurableAxis binsTpcSignal{"binsTpcSignal", {500, 0., 500.}, "TPC signal bins"}; + ConfigurableAxis binsPtTrig{"binsPtTrig", {24, 0, 24}, ""}; + ConfigurableAxis binsPtAssoc{"binsPtAssoc", {20, 0, 10}, ""}; ConfigurableAxis binsMult{"binsMult", {100, 0., 10000.}, "Multiplicity bins"}; ConfigurableAxis binsCent{"binsCent", {100, 0., 100.}, "Centrality bins"}; ConfigurableAxis binsPosZ{"binsPosZ", {100, -10., 10.}, "Primary vertex z coordinate"}; ConfigurableAxis binsEta{"binsEta", {50, -2., 2.}, "Eta bins"}; ConfigurableAxis binsPhi{"binsPhi", {64, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, "Phi bins"}; - ConfigurableAxis binsDeltaEta{"binsDeltaEta", {100, -2., 2.}, "Delta Eta bins"}; - ConfigurableAxis binsDeltaPhi{"binsDeltaPhi", {64, -3., 3.}, "Delta Phi bins"}; - ConfigurableAxis binsMlOne{"binsMlOne", {100, 0., 1.}, ""}; - ConfigurableAxis binsMlTwo{"binsMlTwo", {100, 0., 1.}, ""}; ConfigurableAxis binsDca{"binsDca", {200, -1., 1.}, ""}; HistogramRegistry registry{"registry", {}}; void init(InitContext&) { - if (doprocessDplusSameEvent || doprocessDplusMixedEvent) { - massCharm = o2::constants::physics::MassDPlus; - } else if (doprocessDsSameEvent || doprocessDsMixedEvent) { - massCharm = o2::constants::physics::MassDS; - } else if (doprocessD0SameEvent || doprocessD0MixedEvent) { - massCharm = o2::constants::physics::MassD0; - } else if (doprocessLcSameEvent || doprocessLcMixedEvent) { - massCharm = o2::constants::physics::MassLambdaCPlus; - } else if (doprocessHadronHadronSameEvent || doprocessHadronHadronMixedEvent) { - LOG(info) << "Charm mass not set, processing Hadron-Hadron case"; - } else { - LOG(fatal) << "No decay channel selected to process"; - } hfEvSel.addHistograms(registry); // collision monitoring ccdb->setURL(ccdbUrl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); + std::vector classesMl; + for (const auto scoreClass : classMl.value) { + if (scoreClass < 0) { + LOG(fatal) << "ML class indexes must be non-negative"; + return; + } + classesMl.push_back(static_cast(scoreClass)); + } + + if (classesMl.size() > 2) { + LOG(fatal) << "At most two ML classes can be selected"; + return; + } + + hfcFlowHelper.init(ccdb, registry, classesMl, yCandMax.value); + + const AxisSpec axisTpcSignal = {binsTpcSignal, "TPC Signal"}; const AxisSpec axisCent = {binsCent, "Centrality"}; const AxisSpec axisMult = {binsMult, "Multiplicity"}; const AxisSpec axisPosZ = {binsPosZ, "PosZ"}; const AxisSpec axisEta = {binsEta, "#it{#eta}"}; const AxisSpec axisPhi = {binsPhi, "#it{#varphi}"}; - const AxisSpec axisPtTrig = {std::vector(binsPtTrig), "#it{p}_{T} Trig (GeV/#it{c})"}; - const AxisSpec axisPtAssoc = {std::vector(binsPtAssoc), "#it{p}_{T} Assoc (GeV/#it{c})"}; + const AxisSpec axisPtTrig = {binsPtTrig, "#it{p}_{T} Trig (GeV/#it{c})"}; + const AxisSpec axisPtAssoc = {binsPtAssoc, "#it{p}_{T} Assoc (GeV/#it{c})"}; const AxisSpec axisDcaXY = {binsDca, "DCA XY (cm)"}; const AxisSpec axisDcaZ = {binsDca, "DCA Z (cm)"}; @@ -197,105 +203,12 @@ struct HfDerivedDataCreatorCorrelationsReduced { // Setup pt-dependent DCAxy cut function if (usePtDiffDcaXYCut) { - funcDcaXYPtCutPrimTrk = new TF1("funcDcaXYPtCutPrimTrk", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data()), 0.001, 100); + funcDcaXYPtCutPrimTrk = std::make_unique("funcDcaXYPtCutPrimTrk", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data()), 0.001, 100.); funcDcaXYPtCutPrimTrk->SetParameter(0, dcaXYTrkNSigmaMax); LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data())); } }; // end init - /// Get charm hadron candidate mass - /// \param candidate is the charm hadron candidate - template - double getCandMass(const TCand& candidate) - { - if constexpr (CandType == CandidateType::DsToKKPi) { - return HfHelper::invMassDsToKKPi(candidate); - } - if constexpr (CandType == CandidateType::DsToPiKK) { - return HfHelper::invMassDsToPiKK(candidate); - } - if constexpr (CandType == CandidateType::DplusToPiKPi) { - return HfHelper::invMassDplusToPiKPi(candidate); - } - if constexpr (CandType == CandidateType::D0ToPiK) { - return HfHelper::invMassD0ToPiK(candidate); - } - if constexpr (CandType == CandidateType::D0ToKPi) { - return HfHelper::invMassD0barToKPi(candidate); - } - if constexpr (CandType == CandidateType::LcToPKPi) { - return HfHelper::invMassLcToPKPi(candidate); - } - if constexpr (CandType == CandidateType::LcToPiKP) { - return HfHelper::invMassLcToPiKP(candidate); - } - return -1.; - } - - /// Get charm hadron bdt scores - /// \param candidate is the charm hadron candidate - template - std::array getCandMlScores(const TCand& candidate) - { - std::array outputMl{-999.f, -999.f}; - if constexpr (CandType == CandidateType::DsToKKPi) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbDsToKKPi()[classMl->at(iclass)]; - } - } - if constexpr (CandType == CandidateType::DsToPiKK) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbDsToPiKK()[classMl->at(iclass)]; - } - } - if constexpr (CandType == CandidateType::DplusToPiKPi) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbDplusToPiKPi()[classMl->at(iclass)]; - } - } - if constexpr (CandType == CandidateType::D0ToPiK) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbD0()[classMl->at(iclass)]; - } - } - if constexpr (CandType == CandidateType::D0ToKPi) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbD0bar()[classMl->at(iclass)]; - } - } - if constexpr (CandType == CandidateType::LcToPKPi) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbLcToPKPi()[classMl->at(iclass)]; - } - } - if constexpr (CandType == CandidateType::LcToPiKP) { - for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) { - outputMl[iclass] = candidate.mlProbLcToPiKP()[classMl->at(iclass)]; - } - } - return outputMl; - } - - /// Cut on rapidity of the candidate - /// \param candidate is the charm hadron candidate - template - bool cutCandRapidity(const TCand& candidate) - { - double y = 0.0; - if constexpr (CandType == CandidateType::DsToKKPi || CandType == CandidateType::DsToPiKK) { - y = HfHelper::yDs(candidate); - } else if constexpr (CandType == CandidateType::DplusToPiKPi) { - y = HfHelper::yDplus(candidate); - } else if constexpr (CandType == CandidateType::D0ToPiK || CandType == CandidateType::D0ToKPi) { - y = HfHelper::yD0(candidate); - } else if constexpr (CandType == CandidateType::LcToPKPi || CandType == CandidateType::LcToPiKP) { - y = HfHelper::yLc(candidate); - } else { - return true; - } - return std::fabs(y) < yCandMax; - } - /// Check event selections for collision and fill the collision table /// \param collision is the collision template @@ -325,7 +238,7 @@ struct HfDerivedDataCreatorCorrelationsReduced { /// Checks if the trigger cand-associated track pair can be accepted for SE correlation /// \param assTrk is the associated track /// \param cand is the trigger candidate - template + template bool acceptSameEvtPair(TAssocTrk const& assTrk, TCand const& cand, double deltaEta) { if (std::abs(deltaEta) <= deltaEtaAbsMin || std::abs(deltaEta) > deltaEtaAbsMax) { @@ -337,7 +250,7 @@ struct HfDerivedDataCreatorCorrelationsReduced { } int const trackGlobalIndex = assTrk.globalIndex(); - if constexpr (CandType == CandidateType::Hadron) { + if constexpr (CandType == TriggerChannels::Hadron) { if (!cand.isGlobalTrackWoDCA() || cand.tpcNClsCrossedRows() < tpcNClsCrossedRowsMin) { return false; } @@ -363,24 +276,25 @@ struct HfDerivedDataCreatorCorrelationsReduced { /// \param assTrks are the associated tracks /// \param cent is the collision centrality /// \param mult is the collision multiplicity - template + template void fillSameEvent(TTrigCands const& trigCands, TAssocTrks const& assTrks, const float cent, const float mult) { for (const auto& trigCand : trigCands) { - if (!cutCandRapidity(trigCand)) { + if (!hfcFlowHelper.cutCandRapidity(trigCand)) { continue; } double trigCandPt = trigCand.pt(); registry.fill(HIST("hPhiVsPtTrig"), RecoDecay::constrainAngle(trigCand.phi(), -o2::constants::math::PIHalf), trigCandPt); registry.fill(HIST("hEtaVsPtTrig"), trigCand.eta(), trigCandPt); - if constexpr (CandType == CandidateType::Hadron) { + if constexpr (CandType == TriggerChannels::Hadron) { rowTrigHads(rowCollisions.lastIndex(), trigCandPt, trigCand.tpcNClsCrossedRows(), trigCand.itsClusterMap(), trigCand.itsNCls(), trigCand.dcaXY(), trigCand.dcaZ()); } else { - std::array outputMl = getCandMlScores(trigCand); - rowTrigCharms(rowCollisions.lastIndex(), trigCandPt, getCandMass(trigCand), outputMl[0], outputMl[1]); + std::array outputMl = hfcFlowHelper.getCandMlScores(trigCand); + rowTrigCharms(rowCollisions.lastIndex(), trigCandPt, hfcFlowHelper.getCandMass(trigCand), outputMl[0], outputMl[1]); + rowTrigBases(trigCand.phi(), trigCand.eta()); } for (const auto& assTrk : assTrks) { @@ -407,13 +321,14 @@ struct HfDerivedDataCreatorCorrelationsReduced { registry.fill(HIST("hPhiVsPtAssoc"), RecoDecay::constrainAngle(assTrk.phi(), -o2::constants::math::PIHalf), assTrkPt); registry.fill(HIST("hEtaVsPtAssoc"), assTrk.eta(), assTrkPt); registry.fill(HIST("hDcaXYVsPtVsCentAssoc"), assTrkPt, assTrk.dcaXY(), cent); - registry.fill(HIST("hDcaZVsPtVsCentAssoc"), assTrkPt, assTrk.dcaZ(), cent); registry.fill(HIST("hDcaXYVsPtVsMultAssoc"), assTrkPt, assTrk.dcaXY(), mult); + registry.fill(HIST("hDcaZVsPtVsCentAssoc"), assTrkPt, assTrk.dcaZ(), cent); registry.fill(HIST("hDcaZVsPtVsMultAssoc"), assTrkPt, assTrk.dcaZ(), mult); double deltaPhi = RecoDecay::constrainAngle(assTrk.phi() - trigCand.phi(), -o2::constants::math::PIHalf); + rowAssocBases(rowCollisions.lastIndex(), assTrk.phi(), assTrk.eta(), assTrkPt); rowAssocTrkSels(assTrk.tpcNClsCrossedRows(), assTrk.itsClusterMap(), assTrk.itsNCls(), assTrk.dcaXY(), assTrk.dcaZ()); - if constexpr (CandType == CandidateType::Hadron) { + if constexpr (CandType == TriggerChannels::Hadron) { rowSEHadHadPairs(rowCollisions.lastIndex(), rowTrigHads.lastIndex(), assTrkPt, deltaEta, deltaPhi); } else { rowSECharmHadPairs(rowCollisions.lastIndex(), rowTrigCharms.lastIndex(), assTrkPt, deltaEta, deltaPhi); @@ -424,19 +339,19 @@ struct HfDerivedDataCreatorCorrelationsReduced { /// Fill charm hadron tables for mixed-event /// \param trigCands are the charm trigger candidates - template + template void fillCharmMixedEvent(TTrigCands const& trigCands) { for (const auto& trigCand : trigCands) { - if (!cutCandRapidity(trigCand)) { + if (!hfcFlowHelper.cutCandRapidity(trigCand)) { continue; } registry.fill(HIST("hPhiVsPtTrig"), RecoDecay::constrainAngle(trigCand.phi(), -o2::constants::math::PIHalf), trigCand.pt()); registry.fill(HIST("hEtaVsPtTrig"), trigCand.eta(), trigCand.pt()); - std::array outputMl = getCandMlScores(trigCand); + std::array outputMl = hfcFlowHelper.getCandMlScores(trigCand); rowTrigBases(trigCand.phi(), trigCand.eta()); - rowTrigCharms(rowCollisions.lastIndex(), trigCand.pt(), getCandMass(trigCand), outputMl[0], outputMl[1]); + rowTrigCharms(rowCollisions.lastIndex(), trigCand.pt(), hfcFlowHelper.getCandMass(trigCand), outputMl[0], outputMl[1]); } } @@ -444,7 +359,7 @@ struct HfDerivedDataCreatorCorrelationsReduced { /// \param assTrks are the associated tracks /// \param cent is the collision centrality /// \param mult is the collision multiplicity - template + template void fillTrkMixedEvent(TAssocTrks const& assTrks, const float cent, const float mult) @@ -479,10 +394,11 @@ struct HfDerivedDataCreatorCorrelationsReduced { } } - // Dplus with ML selections - void processDplusSameEvent(CollsWithCentMult::iterator const& coll, - CandDplusData const& candsDplus, - TracksData const& tracks) + // CHARM-HADRON CORRELATIONS PROCESS FUNCTIONS --- DATA + // Dplus data + void processDataDplusSE(DataCollsWithCentMult::iterator const& coll, + CandDplusData const& candsDplus, + TracksData const& tracks) { if (forceCharmInCollision && candsDplus.size() < 1) { return; @@ -492,14 +408,14 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillSameEvent(candsDplus, tracks, cent, mult); + fillSameEvent(candsDplus, tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDplusSameEvent, "Process Same Event for Dplus candidates", true); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataDplusSE, "Process Same Event for Dplus candidates", true); - // Dplus with ML selections - void processDplusMixedEvent(CollsWithCentMult::iterator const& coll, - CandDplusData const& candsDplus, - TracksData const& tracks) + // Dplus data + void processDataDplusME(DataCollsWithCentMult::iterator const& coll, + CandDplusData const& candsDplus, + TracksData const& tracks) { if (forceCharmInCollision && candsDplus.size() < 1) { return; @@ -509,18 +425,18 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillCharmMixedEvent(candsDplus); - fillTrkMixedEvent(tracks, cent, mult); + fillCharmMixedEvent(candsDplus); + fillTrkMixedEvent(tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDplusMixedEvent, "Process Mixed Event for Dplus candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataDplusME, "Process Mixed Event for Dplus candidates", false); - // Ds with ML selections - void processDsSameEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks, - CandDsData const&) + // Ds data + void processDataDsSE(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks, + CandDsData const&) { - auto candsDsToPiKK = selectedDsToPiKK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); - auto candsDsToKKPi = selectedDsToKKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsDsToPiKK = selDsToPiKK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsDsToKKPi = selDsToKKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); if (forceCharmInCollision && candsDsToPiKK.size() < 1 && candsDsToKKPi.size() < 1) { return; } @@ -529,18 +445,18 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillSameEvent(candsDsToPiKK, tracks, cent, mult); - fillSameEvent(candsDsToKKPi, tracks, cent, mult); + fillSameEvent(candsDsToPiKK, tracks, cent, mult); + fillSameEvent(candsDsToKKPi, tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDsSameEvent, "Process Same Event for Ds candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataDsSE, "Process Same Event for Ds candidates", false); - // Ds with ML selections - void processDsMixedEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks, - CandDsData const&) + // Ds data + void processDataDsME(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks, + CandDsData const&) { - auto candsDsToPiKK = selectedDsToPiKK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); - auto candsDsToKKPi = selectedDsToKKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsDsToPiKK = selDsToPiKK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsDsToKKPi = selDsToKKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); if (forceCharmInCollision && candsDsToPiKK.size() < 1 && candsDsToKKPi.size() < 1) { return; } @@ -549,19 +465,19 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillCharmMixedEvent(candsDsToPiKK); - fillCharmMixedEvent(candsDsToKKPi); - fillTrkMixedEvent(tracks, cent, mult); + fillCharmMixedEvent(candsDsToPiKK); + fillCharmMixedEvent(candsDsToKKPi); + fillTrkMixedEvent(tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDsMixedEvent, "Process Mixed Event for Ds candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataDsME, "Process Mixed Event for Ds candidates", false); - // D0 with ML selections - void processD0SameEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks, - CandD0Data const&) + // D0 data + void processDataD0SE(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks, + CandD0Data const&) { - auto candsD0ToPiK = selectedD0ToPiK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); - auto candsD0ToKPi = selectedD0ToKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsD0ToPiK = selD0ToPiK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsD0ToKPi = selD0ToKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); if (forceCharmInCollision && candsD0ToPiK.size() < 1 && candsD0ToKPi.size() < 1) { return; } @@ -570,18 +486,18 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillSameEvent(candsD0ToPiK, tracks, cent, mult); - fillSameEvent(candsD0ToKPi, tracks, cent, mult); + fillSameEvent(candsD0ToPiK, tracks, cent, mult); + fillSameEvent(candsD0ToKPi, tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processD0SameEvent, "Process Same Event for D0 candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataD0SE, "Process Same Event for D0 candidates", false); - // D0 with ML selections - void processD0MixedEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks, - CandD0Data const&) + // D0 data + void processDataD0ME(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks, + CandD0Data const&) { - auto candsD0ToPiK = selectedD0ToPiK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); - auto candsD0ToKPi = selectedD0ToKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsD0ToPiK = selD0ToPiK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsD0ToKPi = selD0ToKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); if (forceCharmInCollision && candsD0ToPiK.size() < 1 && candsD0ToKPi.size() < 1) { return; } @@ -590,19 +506,19 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillCharmMixedEvent(candsD0ToPiK); - fillCharmMixedEvent(candsD0ToKPi); - fillTrkMixedEvent(tracks, cent, mult); + fillCharmMixedEvent(candsD0ToPiK); + fillCharmMixedEvent(candsD0ToKPi); + fillTrkMixedEvent(tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processD0MixedEvent, "Process Mixed Event for D0 candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataD0ME, "Process Mixed Event for D0 candidates", false); - // Lc with ML selections - void processLcSameEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks, - CandLcData const&) + // Lc data + void processDataLcSE(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks, + CandLcData const&) { - auto candsLcToPKPi = selectedLcToPKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); - auto candsLcToPiKP = selectedLcToPiKP->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsLcToPKPi = selLcToPKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsLcToPiKP = selLcToPiKP->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); if (forceCharmInCollision && candsLcToPKPi.size() < 1 && candsLcToPiKP.size() < 1) { return; } @@ -611,18 +527,18 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillSameEvent(candsLcToPKPi, tracks, cent, mult); - fillSameEvent(candsLcToPiKP, tracks, cent, mult); + fillSameEvent(candsLcToPKPi, tracks, cent, mult); + fillSameEvent(candsLcToPiKP, tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processLcSameEvent, "Process Same Event for Lc candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataLcSE, "Process Same Event for Lc candidates", false); - // Lc with ML selections - void processLcMixedEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks, - CandLcData const&) + // Lc data + void processDataLcME(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks, + CandLcData const&) { - auto candsLcToPKPi = selectedLcToPKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); - auto candsLcToPiKP = selectedLcToPiKP->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsLcToPKPi = selLcToPKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsLcToPiKP = selLcToPiKP->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); if (forceCharmInCollision && candsLcToPKPi.size() < 1 && candsLcToPiKP.size() < 1) { return; } @@ -631,37 +547,226 @@ struct HfDerivedDataCreatorCorrelationsReduced { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillCharmMixedEvent(candsLcToPKPi); - fillCharmMixedEvent(candsLcToPiKP); - fillTrkMixedEvent(tracks, cent, mult); + fillCharmMixedEvent(candsLcToPKPi); + fillCharmMixedEvent(candsLcToPiKP); + fillTrkMixedEvent(tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataLcME, "Process Mixed Event for Lc candidates", false); + + // CHARM-HADRON CORRELATIONS PROCESS FUNCTIONS --- MC + // Dplus MC + void processMcDplusSE(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandDplusMc const&) + { + auto candsRecoDplus = recoSelDPlusToPiKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoDplus.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillSameEvent(candsRecoDplus, tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcDplusSE, "Process Same Event for Dplus candidates", false); + + // Dplus MC + void processMcDplusME(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandDplusMc const&) + { + auto candsRecoDplus = recoSelDPlusToPiKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoDplus.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillCharmMixedEvent(candsRecoDplus); + fillTrkMixedEvent(tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcDplusME, "Process Mixed Event for Dplus candidates", false); + + // Ds MC + void processMcDsSE(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandDsMc const&) + { + auto candsRecoDsToPiKK = recoSelDsToPiKK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsRecoDsToKKPi = recoSelDsToKKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoDsToPiKK.size() < 1 && candsRecoDsToKKPi.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillSameEvent(candsRecoDsToPiKK, tracks, cent, mult); + fillSameEvent(candsRecoDsToKKPi, tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcDsSE, "Process Same Event for Ds candidates", false); + + // Ds MC + void processMcDsME(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandDsMc const&) + { + auto candsRecoDsToPiKK = recoSelDsToPiKK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsRecoDsToKKPi = recoSelDsToKKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoDsToPiKK.size() < 1 && candsRecoDsToKKPi.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillCharmMixedEvent(candsRecoDsToPiKK); + fillCharmMixedEvent(candsRecoDsToKKPi); + fillTrkMixedEvent(tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcDsME, "Process Mixed Event for Ds candidates", false); + + // D0 MC + void processMcD0SE(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandD0Mc const&) + { + auto candsRecoD0ToPiK = recoSelD0ToPiK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsRecoD0ToKPi = recoSelD0ToKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoD0ToPiK.size() < 1 && candsRecoD0ToKPi.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillSameEvent(candsRecoD0ToPiK, tracks, cent, mult); + fillSameEvent(candsRecoD0ToKPi, tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcD0SE, "Process Same Event for D0 candidates", false); + + // D0 MC + void processMcD0ME(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandD0Mc const&) + { + auto candsRecoD0ToPiK = recoSelD0ToPiK->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsRecoD0ToKPi = recoSelD0ToKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoD0ToPiK.size() < 1 && candsRecoD0ToKPi.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillCharmMixedEvent(candsRecoD0ToPiK); + fillCharmMixedEvent(candsRecoD0ToKPi); + fillTrkMixedEvent(tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcD0ME, "Process Mixed Event for D0 candidates", false); + + // Lc MC + void processMcLcSE(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandLcMc const&) + { + auto candsRecoLcToPKPi = recoSelLcToPKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsRecoLcToPiKP = recoSelLcToPiKP->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoLcToPKPi.size() < 1 && candsRecoLcToPiKP.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillSameEvent(candsRecoLcToPKPi, tracks, cent, mult); + fillSameEvent(candsRecoLcToPiKP, tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcLcSE, "Process Same Event for Lc candidates", false); + + // Lc MC + void processMcLcME(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks, + CandLcMc const&) + { + auto candsRecoLcToPKPi = recoSelLcToPKPi->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + auto candsRecoLcToPiKP = recoSelLcToPiKP->sliceByCached(aod::hf_cand::collisionId, coll.globalIndex(), cache); + if (forceCharmInCollision && candsRecoLcToPKPi.size() < 1 && candsRecoLcToPiKP.size() < 1) { + return; + } + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillCharmMixedEvent(candsRecoLcToPKPi); + fillCharmMixedEvent(candsRecoLcToPiKP); + fillTrkMixedEvent(tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcLcME, "Process Mixed Event for Lc candidates", false); + + // HADRON-HADRON CORRELATIONS PROCESS FUNCTIONS --- DATA + // Hadron Hadron Same Event + void processDataHadronHadronSE(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks) + { + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillSameEvent(tracks, tracks, cent, mult); + } + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataHadronHadronSE, "Process Same Event for hadron candidates", false); + + // Hadron Hadron Mixed Event + void processDataHadronHadronME(DataCollsWithCentMult::iterator const& coll, + TracksData const& tracks) + { + float cent{-1.}, mult{-1.}; + if (!checkCollision(coll, cent, mult)) { + return; + } + rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); + fillTrkMixedEvent(tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processLcMixedEvent, "Process Mixed Event for Lc candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processDataHadronHadronME, "Process Mixed Event for hadron candidates", false); + // HADRON-HADRON CORRELATIONS PROCESS FUNCTIONS --- MC // Hadron Hadron Same Event - void processHadronHadronSameEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks) + void processMcHadronHadronSE(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks) { float cent{-1.}, mult{-1.}; if (!checkCollision(coll, cent, mult)) { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillSameEvent(tracks, tracks, cent, mult); + fillSameEvent(tracks, tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processHadronHadronSameEvent, "Process Same Event for hadron candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcHadronHadronSE, "Process Same Event for hadron candidates", false); // Hadron Hadron Mixed Event - void processHadronHadronMixedEvent(CollsWithCentMult::iterator const& coll, - TracksData const& tracks) + void processMcHadronHadronME(McCollsWithCentMult::iterator const& coll, + TracksMc const& tracks) { float cent{-1.}, mult{-1.}; if (!checkCollision(coll, cent, mult)) { return; } rowCollisions(mult, coll.numContrib(), cent, coll.posZ()); - fillTrkMixedEvent(tracks, cent, mult); + fillTrkMixedEvent(tracks, cent, mult); } - PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processHadronHadronMixedEvent, "Process Mixed Event for hadron candidates", false); + PROCESS_SWITCH(HfDerivedDataCreatorCorrelationsReduced, processMcHadronHadronME, "Process Mixed Event for hadron candidates", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGHF/HFC/Tasks/CMakeLists.txt b/PWGHF/HFC/Tasks/CMakeLists.txt index 4551ccd2b0f..8c60c81b29d 100644 --- a/PWGHF/HFC/Tasks/CMakeLists.txt +++ b/PWGHF/HFC/Tasks/CMakeLists.txt @@ -74,6 +74,11 @@ o2physics_add_dpl_workflow(task-flow PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGCFCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(task-correlations-efficiencies + SOURCES taskCorrelationsEfficiencies.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + include_directories(${CMAKE_CURRENT_SOURCE_DIR}/PWGCF/Core) o2physics_add_dpl_workflow(task-charm-hadrons-charm-femto-dream diff --git a/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx b/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx new file mode 100644 index 00000000000..8407dcbde82 --- /dev/null +++ b/PWGHF/HFC/Tasks/taskCorrelationsEfficiencies.cxx @@ -0,0 +1,438 @@ +// 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. + +/// \file taskCorrelationsEfficiencies.cxx +/// \brief Charm hadrons and associated tracks efficiency evaluation for correlations studies +/// \author Marcello Di Costanzo , Politecnico and INFN Torino + +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/DataModel/AliasTables.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "PWGHF/HFC/DataModel/DerivedDataCorrelationTables.h" +#include "PWGHF/HFC/Utils/utilsFlow.h" +#include "PWGHF/Utils/utilsEvSelHf.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::hf_centrality; +using namespace o2::hf_evsel; +using namespace o2::hfc_utils_flow; + +namespace o2::aod +{ +namespace hf_correlations_eff +{ +// Charm-hadron candidate columns +DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c) +DECLARE_SOA_COLUMN(PvContrib, pvContrib, int); //! Number of primary vertex contributors +DECLARE_SOA_COLUMN(Eta, eta, float); //! Pseudorapidity of candidate +DECLARE_SOA_COLUMN(Phi, phi, float); //! Azimuthal angle of candidate +DECLARE_SOA_COLUMN(Origin, origin, int); //! Origin of the candidate (prompt or non-prompt) +DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! ML score of the first configured index +DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! ML score of the second configured index + +// Associated track columns +DECLARE_SOA_COLUMN(PtAssoc, ptAssoc, float); //! Transverse momentum of associated track (GeV/c) +DECLARE_SOA_COLUMN(EtaAssoc, etaAssoc, float); //! Pseudorapidity of associated track +DECLARE_SOA_COLUMN(PhiAssoc, phiAssoc, float); //! Azimuthal angle of associated track +DECLARE_SOA_COLUMN(ZVtx, zVtx, float); //! Primary vertex z position (cm) +DECLARE_SOA_COLUMN(AssTrkSpecie, assTrkSpecie, int); //! Associated track species (0=pi, 1=K, 2=p, 3=el, 4=mu) +DECLARE_SOA_COLUMN(AssTrkOrigin, assTrkOrigin, int); //! Associated track origin (0=primary, 1=secondary) +} // namespace hf_correlations_eff + +DECLARE_SOA_TABLE(HfcMcrigInfo, "AOD", "HFCMCTRIGINFO", + hf_correlations_eff::Pt, + hf_correlations_eff::PvContrib, + hf_correlations_eff::Eta, + hf_correlations_eff::Phi, + hf_correlations_eff::Origin, + hf_correlations_eff::MlScore0, + hf_correlations_eff::MlScore1); + +DECLARE_SOA_TABLE(HfcMcAssocInfo, "AOD", "HFCMCASSOCINFO", + hf_correlations_eff::PtAssoc, + hf_correlations_eff::EtaAssoc, + hf_correlations_eff::PhiAssoc, + hf_correlations_eff::ZVtx, + hf_correlations_eff::AssTrkSpecie, + hf_correlations_eff::AssTrkOrigin); +} // namespace o2::aod + +/// Code to obtain efficiency correction factors for (Delta Eta, Delta Phi) correlation maps +struct HfTaskCorrelationsEfficiencies { + + Produces rowTrigMc; + Produces rowAssocMc; + Produces rowAssocTrkSels; + + Configurable selectionFlag{"selectionFlag", 15, "Selection Flag for hadron (ML score tables are required to run the task)"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> classMl{"classMl", {0, 2}, "Indexes of BDT scores to be stored. Two indexes max."}; + Configurable yCandMax{"yCandMax", 0.8, "max. cand. rapidity"}; + Configurable tpcNClsCrossedRowsMin{"tpcNClsCrossedRowsMin", 70, "min. TPC crossed rows for associated tracks"}; + Configurable etaTrkMax{"etaTrkMax", 1., "max. track eta"}; + Configurable ptTrkMin{"ptTrkMin", 0.3, "min. track pT"}; + Configurable ptTrkMax{"ptTrkMax", 3., "max. track pT"}; + Configurable dcaXYTrkMax{"dcaXYTrkMax", 1., "max. track DCA XY"}; + Configurable dcaZTrkMax{"dcaZTrkMax", 1., "max. track DCA Z"}; + Configurable usePtDiffDcaXYCut{"usePtDiffDcaXYCut", false, "Use pt-differential DCAxy cut for associated tracks"}; + Configurable dcaXYTrkNSigmaMax{"dcaXYTrkNSigmaMax", 7, "Cut on number of sigma deviations from expected DCA in the transverse direction"}; + Configurable dcaXYPtPrimTrkFunc{"dcaXYPtPrimTrkFunc", "(0.0026+0.005/(x^1.01))", "Functional form of pt-dependent DCAxy cut"}; + + // Configurables for tree filling + Configurable fillRecoTrees{"fillRecoTrees", false, "Fill trees of reconstructed charmed hadrons and associated tracks"}; + Configurable downSampleTracks{"downSampleTracks", 1., "Fraction of tracks to keep"}; + Configurable ptMaxDownSampleTracks{"ptMaxDownSampleTracks", 10., "Maximum pt for the application of the downsampling factor"}; + Configurable downSampleCharm{"downSampleCharm", 1., "Fraction of pairs to keep"}; + Configurable ptMaxDownSampleCharm{"ptMaxDownSampleCharm", 10., "Maximum pt for the application of the downsampling factor"}; + + HfcFlowHelper hfcFlowHelper; // flow utilities + o2::framework::Service ccdb{}; + SliceCache cache; + + std::unique_ptr funcDcaXYPtCutPrimTrk{nullptr}; + + using RecoCollisions = soa::Join; + using McCandsDs = soa::Filtered>; + using McCandsDplus = soa::Filtered>; + using McCandsD0 = soa::Filtered>; + using McCandsLc = soa::Filtered>; + using McTracks = soa::Filtered>; + + Filter filterSelectDsCands = aod::hf_sel_candidate_ds::isSelDsToKKPi >= selectionFlag || aod::hf_sel_candidate_ds::isSelDsToPiKK >= selectionFlag; + Filter filterSelectDplusCands = aod::hf_sel_candidate_dplus::isSelDplusToPiKPi >= selectionFlag; + Filter filterSelectD0Cands = aod::hf_sel_candidate_d0::isSelD0 >= selectionFlag || aod::hf_sel_candidate_d0::isSelD0bar >= selectionFlag; + Filter filterSelectLcCands = aod::hf_sel_candidate_lc::isSelLcToPKPi >= selectionFlag || aod::hf_sel_candidate_lc::isSelLcToPiKP >= selectionFlag; + Filter filterSelectTrk = (nabs(aod::track::eta) < etaTrkMax) && (aod::track::pt > ptTrkMin) && (aod::track::pt < ptTrkMax) && (nabs(aod::track::dcaXY) < dcaXYTrkMax) && (nabs(aod::track::dcaZ) < dcaZTrkMax); + + Preslice candsDsPerColl = aod::hf_cand::collisionId; + Preslice candsDplusPerColl = aod::hf_cand::collisionId; + Preslice candsD0PerColl = aod::hf_cand::collisionId; + Preslice candsLcPerColl = aod::hf_cand::collisionId; + Preslice trackIndicesPerColl = aod::track::collisionId; + + PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; + + ConfigurableAxis binsPtCharm{"binsPtCharm", {VARIABLE_WIDTH, 0., 2., 4., 6., 8., 24.}, ""}; + ConfigurableAxis binsPtAssoc{"binsPtAssoc", {20, 0, 10}, ""}; + ConfigurableAxis binsZVtx{"binsZVtx", {120, -11, 11}, ""}; + ConfigurableAxis binsPvContrib{"binsPvContrib", {5001, -0.5, 5000.5}, ""}; + ConfigurableAxis binsEta{"binsEta", {50, -2., 2.}, ""}; + ConfigurableAxis binsPhi{"binsPhi", {64, 0, 6.28}, ""}; + ConfigurableAxis binsMlOne{"binsMlOne", {100, 0., 1.}, ""}; + ConfigurableAxis binsMlTwo{"binsMlTwo", {100, 0., 1.}, ""}; + + HistogramRegistry registry{"registry", {}}; + + void init(InitContext&) + { + + ccdb->setURL(ccdbUrl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + + std::vector classesMl; + for (const auto scoreClass : classMl.value) { + if (scoreClass < 0) { + LOG(fatal) << "ML class indexes must be non-negative"; + return; + } + classesMl.push_back(static_cast(scoreClass)); + } + + if (classesMl.size() > 2) { + LOG(fatal) << "At most two ML classes can be selected"; + return; + } + + hfcFlowHelper.init(ccdb, registry, classesMl, yCandMax.value); + + // Setup pt-dependent DCAxy cut function + if (usePtDiffDcaXYCut) { + funcDcaXYPtCutPrimTrk = std::make_unique("funcDcaXYPtCutPrimTrk", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data()), 0.001, 100.); + funcDcaXYPtCutPrimTrk->SetParameter(0, dcaXYTrkNSigmaMax); + LOGF(info, "DCAxy pt-dependence function: %s", Form("[0]*%s", dcaXYPtPrimTrkFunc.value.data())); + } + + // Add efficiency sparses + // Axes: mass, pt, pv contributors, pseudo-rapidity, ml_score_0, ml_score_1 + const AxisSpec axisPtTrig = {binsPtCharm, "#it{p}_{T} Trig (GeV/#it{c})"}; + const AxisSpec axisPvContrib = {binsPvContrib, "PV contributors"}; + const AxisSpec axisEtaTrig = {binsEta, "#it{#eta}"}; + const AxisSpec axisPhiTrig = {binsPhi, "#it{#varphi}"}; + const AxisSpec axisMlScore0 = {binsMlOne, "ML score 0"}; + const AxisSpec axisMlScore1 = {binsMlTwo, "ML score 1"}; + std::vector axesCharmReco = {axisPtTrig, axisPvContrib, axisEtaTrig, axisPhiTrig, axisMlScore0, axisMlScore1}; + + registry.add("hCharmPromptGen", "Generated prompt charm-hadron candidates;;", HistType::kTHnSparseF, {{axisPtTrig}, {axisEtaTrig}, {axisPhiTrig}}); + registry.add("hCharmFdGen", "Generated non-prompt charm-hadron candidates;;", HistType::kTHnSparseF, {{axisPtTrig}, {axisEtaTrig}, {axisPhiTrig}}); + registry.add("hCharmPromptReco", "Reconstructed prompt charm-hadron candidates;;", HistType::kTHnSparseF, axesCharmReco); + registry.add("hCharmFdReco", "Reconstructed non-prompt charm-hadron candidates;;", HistType::kTHnSparseF, axesCharmReco); + + // Axes: pt, eta, phi, z_vtx, specie + const AxisSpec axisPtAssoc = {binsPtAssoc, "#it{p}_{T} Assoc (GeV/#it{c})"}; + const AxisSpec axisEtaAssoc = {binsEta, "#it{#eta}"}; + const AxisSpec axisPhiAssoc = {binsPhi, "#it{#varphi}"}; + const AxisSpec axisZVtx = {binsZVtx, "#it{z}_{vtx}"}; + const AxisSpec axisAssTrkSpecie = {static_cast(NAssTracksTypes), -0.5, static_cast(NAssTracksTypes) - 0.5, "Particle species"}; + const AxisSpec axisTrackOrigin = {2, -0.5, 1.5, "Particle origin"}; + std::vector axesAssocTracks = {axisPtAssoc, axisEtaAssoc, axisPhiAssoc, axisZVtx, axisAssTrkSpecie, axisTrackOrigin}; + registry.add("hAssTrackReco", "Reconstructed associated tracks;;", HistType::kTHnSparseF, axesAssocTracks); + registry.add("hAssTrackGen", "Generated associated tracks;;", HistType::kTHnSparseF, axesAssocTracks); + } + + /// Match charm-hadron based on its reconstruction flag + /// \param assTrk is the associated track + /// \return the associated track specie index + template + bool isMcMatched(TCharmHad const& cand) + { + int flagMc{-1}; + if constexpr (IsReco) { + flagMc = std::abs(cand.flagMcMatchRec()); + } else { + flagMc = std::abs(cand.flagMcMatchGen()); + } + + if constexpr (CandType == TriggerChannels::D0ToPiK) { + return flagMc == static_cast(hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK); + } else if constexpr (CandType == TriggerChannels::DplusToPiKPi) { + return flagMc == static_cast(hf_decay::hf_cand_3prong::DecayChannelMain::DplusToPiKPi); + } else if constexpr (CandType == TriggerChannels::LcToPKPi) { + return flagMc == static_cast(hf_decay::hf_cand_3prong::DecayChannelMain::LcToPKPi); + } else { // Ds + int8_t reso{0}; + if constexpr (IsReco) { + reso = std::abs(cand.flagMcDecayChanRec()); + } else { + reso = std::abs(cand.flagMcDecayChanGen()); + } + return flagMc == static_cast(hf_decay::hf_cand_3prong::DecayChannelMain::DsToPiKK) && + reso == static_cast(hf_decay::hf_cand_3prong::DecayChannelResonant::DsToPhiPi); + } + } + + /// Fill info for charm-hadron efficiencies + /// \param trigCands are the trigger candidates + /// \param assTrks are the associated tracks + /// \param cent is the collision centrality + /// \param mult is the collision multiplicity + template + void fillCharmEfficiencies(TCharmCands const& charmCands, int numContrib) + { + for (const auto& cand : charmCands) { + if (!hfcFlowHelper.cutCandRapidity(cand)) { + continue; + } + + if (!isMcMatched(cand)) { + continue; + } + + // Fill the sparse + std::array outputMl = hfcFlowHelper.getCandMlScores(cand); + if (cand.originMcRec() == RecoDecay::OriginType::Prompt) { + registry.fill(HIST("hCharmPromptReco"), cand.pt(), numContrib, cand.eta(), cand.phi(), outputMl[0], outputMl[1]); + } + if (cand.originMcRec() == RecoDecay::OriginType::NonPrompt) { + registry.fill(HIST("hCharmFdReco"), cand.pt(), numContrib, cand.eta(), cand.phi(), outputMl[0], outputMl[1]); + } + + if (fillRecoTrees) { + // Downsample charm candidates based on pt + if (downSampleCharm < 1.) { + float const pseudoRndm = cand.pt() * 1000. - static_cast(cand.pt() * 1000); + if (cand.pt() < ptMaxDownSampleCharm && pseudoRndm >= downSampleCharm) { + continue; + } + } + rowTrigMc(cand.pt(), numContrib, cand.eta(), cand.phi(), cand.originMcRec(), outputMl[0], outputMl[1]); + } + } + } + + // Dplus efficiency + void processEffDplus(RecoCollisions const& recoCollisions, + McCandsDplus const& candsDplus, + soa::Join const& genParticles, + aod::McCollisions const&) + { + for (const auto& coll : recoCollisions) { + auto candsInColl = candsDplus.sliceBy(candsDplusPerColl, coll.globalIndex()); + if (candsInColl.size() < 1) { + continue; + } + float cent{-1.}, mult{-1.}; + if (!hfcFlowHelper.checkCollision(coll, cent, mult)) { + continue; + } + fillCharmEfficiencies(candsInColl, coll.numContrib()); + } + + for (const auto& genCand : genParticles) { + if (!genCand.has_mcCollision()) { + continue; + } + if (std::abs(genCand.y()) > yCandMax) { + continue; + } + if (!isMcMatched(genCand)) { + continue; + } + + if (genCand.originMcGen() == RecoDecay::OriginType::Prompt) { + registry.fill(HIST("hCharmPromptGen"), genCand.pt(), genCand.eta(), genCand.phi()); + } + if (genCand.originMcGen() == RecoDecay::OriginType::NonPrompt) { + registry.fill(HIST("hCharmFdGen"), genCand.pt(), genCand.eta(), genCand.phi()); + } + } + } + PROCESS_SWITCH(HfTaskCorrelationsEfficiencies, processEffDplus, "Process efficiencies for Dplus candidates", false); + + /// Match associated tracks based on their PDG code + /// \param assTrk is the associated track + /// \return the associated track specie index + template + int getAssocTrkSpecie(TAssocTrk const& assTrk) + { + switch (std::abs(assTrk.pdgCode())) { + case kPiPlus: + return static_cast(AssTracksTypes::Pi); + case kKPlus: + return static_cast(AssTracksTypes::Ka); + case kProton: + return static_cast(AssTracksTypes::Pr); + case kElectron: + return static_cast(AssTracksTypes::El); + case kMuonMinus: + return static_cast(AssTracksTypes::Mu); + default: + return -1; + } + } + + // Associated Tracks + void processEffAssTracks(RecoCollisions const& recoCollisions, + McTracks const& recoTracks, + aod::McParticles const& genParticles, + aod::McCollisions const&) + { + for (const auto& coll : recoCollisions) { + float cent{-1.}, mult{-1.}; + if (!hfcFlowHelper.checkCollision(coll, cent, mult)) { + continue; + } + if (!coll.has_mcCollision()) { + continue; + } + auto mcCollision = coll.template mcCollision_as(); + + auto tracksInColl = recoTracks.sliceBy(trackIndicesPerColl, coll.globalIndex()); + for (const auto& assTrk : tracksInColl) { + if (!assTrk.has_mcParticle()) { + continue; + } + auto mcParticle = assTrk.template mcParticle_as(); + int matchedSpecie = getAssocTrkSpecie(mcParticle); + if (matchedSpecie < 0) { + continue; + } + + // DCA selection + double assTrkPt = assTrk.pt(); + if (usePtDiffDcaXYCut) { + float const dcaXYTrkCut = funcDcaXYPtCutPrimTrk->Eval(assTrkPt); + if (std::fabs(assTrk.dcaXY()) > dcaXYTrkCut) { + continue; + } + } + int assTrkOrigin = mcParticle.isPhysicalPrimary() ? 0 : 1; // 0 for primary, 1 for secondary + + std::vector values = {assTrk.pt(), assTrk.eta(), assTrk.phi(), mcCollision.posZ(), static_cast(matchedSpecie), static_cast(assTrkOrigin)}; + + // Fill histograms for associated tracks + registry.get(HIST("hAssTrackReco"))->Fill(values.data()); + + if (fillRecoTrees) { + // Downsample tracks based on pt + if (downSampleTracks < 1.) { + float const pseudoRndm = assTrk.pt() * 1000. - static_cast(assTrk.pt() * 1000); + if (assTrk.pt() < ptMaxDownSampleTracks && pseudoRndm >= downSampleTracks) { + continue; + } + } + rowAssocMc(assTrk.pt(), assTrk.eta(), assTrk.phi(), mcCollision.posZ(), matchedSpecie, assTrkOrigin); + rowAssocTrkSels(assTrk.tpcNClsCrossedRows(), assTrk.itsClusterMap(), assTrk.itsNCls(), assTrk.dcaXY(), assTrk.dcaZ()); + } + } + } + + // Fill histograms for generated associated tracks + for (const auto& genTrk : genParticles) { + int matchedSpecie = getAssocTrkSpecie(genTrk); + if (matchedSpecie < 0) { + continue; + } + if (!genTrk.has_mcCollision()) { + continue; + } + + // Pt and eta selection + if (genTrk.pt() < ptTrkMin || genTrk.pt() > ptTrkMax || std::fabs(genTrk.eta()) > etaTrkMax) { + continue; + } + + int genTrkOrigin = genTrk.isPhysicalPrimary() ? 0 : 1; // 0 for primary, 1 for secondary + auto genCollision = genTrk.template mcCollision_as(); + registry.fill(HIST("hAssTrackGen"), genTrk.pt(), genTrk.eta(), genTrk.phi(), genCollision.posZ(), matchedSpecie, genTrkOrigin); + } + } + PROCESS_SWITCH(HfTaskCorrelationsEfficiencies, processEffAssTracks, "Process efficiencies for associated tracks", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGHF/HFC/Utils/utilsFlow.h b/PWGHF/HFC/Utils/utilsFlow.h new file mode 100644 index 00000000000..52b4885bfe4 --- /dev/null +++ b/PWGHF/HFC/Utils/utilsFlow.h @@ -0,0 +1,227 @@ +// 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. + +/// \file utilsFlow.h +/// \brief Utilities for HFC analyses +/// \author Marcello Di Costanzo + +#ifndef PWGHF_HFC_UTILS_UTILSFLOW_H_ +#define PWGHF_HFC_UTILS_UTILSFLOW_H_ + +#include "PWGHF/Core/CentralityEstimation.h" +#include "PWGHF/Core/DecayChannels.h" +#include "PWGHF/Core/HfHelper.h" +#include "PWGHF/Utils/utilsEvSelHf.h" + +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace o2::hfc_utils_flow +{ + +enum TriggerChannels { + DplusToPiKPi = 0, + DsToKKPi, + DsToPiKK, + D0ToPiK, + D0ToKPi, + LcToPKPi, + LcToPiKP, + Hadron +}; + +enum AssTracksTypes { + Pi = 0, + Ka, + Pr, + El, + Mu, + NAssTracksTypes +}; + +class HfcFlowHelper +{ + + public: + HfcFlowHelper() = default; + + /// Initialize the HFC flow helper + /// \param ccdbService is the CCDB service + /// \param histRegistry is the histogram registry + /// \param mlClasses are the machine learning classes + /// \param maxRapidity is the maximum rapidity + void init(o2::framework::Service& ccdbService, + o2::framework::HistogramRegistry& histRegistry, + std::vector mlClasses, + float maxRapidity) + { + + LOG(info) << "Initializing HFC flow helper!"; + this->hfEvSel.addHistograms(histRegistry); + this->ccdb = &ccdbService; + this->registry = &histRegistry; + this->classMl = std::move(mlClasses); + this->yCandMax = maxRapidity; + } + + /// Get charm hadron candidate mass + /// \param candidate is the charm hadron candidate + template + double getCandMass(const TCand& candidate) + { + if constexpr (CandType == TriggerChannels::DsToKKPi) { + return HfHelper::invMassDsToKKPi(candidate); + } + if constexpr (CandType == TriggerChannels::DsToPiKK) { + return HfHelper::invMassDsToPiKK(candidate); + } + if constexpr (CandType == TriggerChannels::DplusToPiKPi) { + return HfHelper::invMassDplusToPiKPi(candidate); + } + if constexpr (CandType == TriggerChannels::D0ToPiK) { + return HfHelper::invMassD0ToPiK(candidate); + } + if constexpr (CandType == TriggerChannels::D0ToKPi) { + return HfHelper::invMassD0barToKPi(candidate); + } + if constexpr (CandType == TriggerChannels::LcToPKPi) { + return HfHelper::invMassLcToPKPi(candidate); + } + if constexpr (CandType == TriggerChannels::LcToPiKP) { + return HfHelper::invMassLcToPiKP(candidate); + } + return -1.; + } + + /// Get charm hadron bdt scores + /// \param candidate is the charm hadron candidate + template + std::array getCandMlScores(const TCand& candidate) + { + std::array outputMl{-999.f, -999.f}; + if constexpr (CandType == TriggerChannels::DsToKKPi) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbDsToKKPi()[classMl.at(iclass)]; + } + } + if constexpr (CandType == TriggerChannels::DsToPiKK) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbDsToPiKK()[classMl.at(iclass)]; + } + } + if constexpr (CandType == TriggerChannels::DplusToPiKPi) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbDplusToPiKPi()[classMl.at(iclass)]; + } + } + if constexpr (CandType == TriggerChannels::D0ToPiK) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbD0()[classMl.at(iclass)]; + } + } + if constexpr (CandType == TriggerChannels::D0ToKPi) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbD0bar()[classMl.at(iclass)]; + } + } + if constexpr (CandType == TriggerChannels::LcToPKPi) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbLcToPKPi()[classMl.at(iclass)]; + } + } + if constexpr (CandType == TriggerChannels::LcToPiKP) { + for (unsigned int iclass = 0; iclass < classMl.size(); iclass++) { + outputMl[iclass] = candidate.mlProbLcToPiKP()[classMl.at(iclass)]; + } + } + return outputMl; + } + + /// Cut on rapidity of the candidate + /// \param candidate is the charm hadron candidate + template + bool cutCandRapidity(const TCand& candidate) + { + double y = 0.0; + if constexpr (CandType == TriggerChannels::DsToKKPi || CandType == TriggerChannels::DsToPiKK) { + y = HfHelper::yDs(candidate); + } else if constexpr (CandType == TriggerChannels::DplusToPiKPi) { + y = HfHelper::yDplus(candidate); + } else if constexpr (CandType == TriggerChannels::D0ToPiK || CandType == TriggerChannels::D0ToKPi) { + y = HfHelper::yD0(candidate); + } else if constexpr (CandType == TriggerChannels::LcToPKPi || CandType == TriggerChannels::LcToPiKP) { + y = HfHelper::yLc(candidate); + } else { + return true; + } + return std::fabs(y) < yCandMax; + } + + /// Check event selections for collision and fill the collision table + /// \param collision is the collision + template + bool checkCollision(Coll const& collision, float& cent, float& mult) + { + o2::hf_evsel::HfCollisionRejectionMask collRejMask{}; + if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FT0A) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, *ccdb, *registry); + mult = collision.multFT0A(); + } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FT0C) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, *ccdb, *registry); + mult = collision.multFT0C(); + } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FT0M) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, *ccdb, *registry); + mult = collision.multFT0M(); + } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::FV0A) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, *ccdb, *registry); + mult = collision.multFV0A(); + } else if constexpr (CentEst == o2::hf_centrality::CentralityEstimator::None) { + collRejMask = hfEvSel.getHfCollisionRejectionMask(collision, cent, *ccdb, *registry); + mult = -1.f; + } else { + LOG(fatal) << "Centrality estimator not recognized for collision selection"; + std::abort(); + } + hfEvSel.fillHistograms(collision, collRejMask, cent); + return collRejMask == 0; + } + + private: + o2::hf_evsel::HfEventSelection hfEvSel; + o2::framework::Service* ccdb{nullptr}; + o2::framework::HistogramRegistry* registry{nullptr}; + std::vector classMl{}; + float yCandMax{-999.f}; +}; + +} // namespace o2::hfc_utils_flow +#endif // PWGHF_HFC_UTILS_UTILSFLOW_H_