diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h index b866649ab04..b4a1ff97539 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseEfficiencyCorrection.h @@ -16,15 +16,22 @@ #ifndef PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_ #define PWGCF_FEMTOUNIVERSE_CORE_FEMTOUNIVERSEEFFICIENCYCORRECTION_H_ +#include +#include +#include +#include +#include + +#include +#include +#include + #include #include #include +#include -#include "Framework/Configurable.h" -#include "CCDB/BasicCCDBManager.h" -#include "TH1.h" -#include "TH2.h" -#include "TH3.h" +#include "PWGCF/FemtoUniverse/DataModel/FemtoDerived.h" namespace o2::analysis::femto_universe::efficiency_correction { @@ -34,30 +41,17 @@ enum ParticleNo : size_t { }; template -concept isOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO; - -template -consteval auto getHistDim() -> int -{ - if (std::is_same_v) - return 1; - else if (std::is_same_v) - return 2; - else if (std::is_same_v) - return 3; - else - return -1; -} +concept IsOneOrTwo = T == ParticleNo::ONE || T == ParticleNo::TWO; struct EffCorConfigurableGroup : framework::ConfigurableGroup { framework::Configurable confEffCorApply{"confEffCorApply", false, "[Efficiency Correction] Should apply efficiency corrections"}; + framework::Configurable confEffCorFillHist{"confEffCorFillHist", false, "[Efficiency Correction] Should fill histograms for efficiency corrections"}; framework::Configurable confEffCorCCDBUrl{"confEffCorCCDBUrl", "http://alice-ccdb.cern.ch", "[Efficiency Correction] CCDB URL to use"}; framework::Configurable confEffCorCCDBPath{"confEffCorCCDBPath", "", "[Efficiency Correction] CCDB path to histograms"}; framework::Configurable> confEffCorCCDBTimestamps{"confEffCorCCDBTimestamps", {}, "[Efficiency Correction] Timestamps of histograms in CCDB (0 can be used as a placeholder, e.g. when running subwagons)"}; + framework::Configurable confEffCorVariables{"confEffCorVariables", "pt", "[Efficiency Correction] Variables for efficiency correction histogram dimensions (available: 'pt'; 'pt,eta'; 'pt,mult'; 'pt,eta,mult')"}; }; -template - requires std::is_base_of_v class EfficiencyCorrection { public: @@ -65,8 +59,23 @@ class EfficiencyCorrection { } - auto init() -> void + auto init(framework::HistogramRegistry* registry, std::vector axisSpecs) -> void { + shouldFillHistograms = config->confEffCorFillHist; + + histRegistry = registry; + if (shouldFillHistograms) { + for (const auto& suffix : histSuffix) { + auto path = std::format("{}/{}", histDirectory, suffix); + registry->add((path + "/hMCTruth").c_str(), "MCTruth; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hPrimary").c_str(), "Primary; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hSecondary").c_str(), "Secondary; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hMaterial").c_str(), "Material; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hFake").c_str(), "Fake; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + registry->add((path + "/hOther").c_str(), "Other; #it{p}_{T} (GeV/#it{c}); #it{#eta}; Mult", framework::kTH3F, axisSpecs); + } + } + ccdb.setURL(config->confEffCorCCDBUrl); ccdb.setLocalObjectValidityChecking(); ccdb.setFatalWhenNull(false); @@ -86,20 +95,123 @@ class EfficiencyCorrection continue; } - hLoaded[idx] = timestamp > 0 ? loadHistFromCCDB(timestamp) : nullptr; + if (timestamp > 0) { + switch (getDimensionFromVariables()) { + case 1: + hLoaded[idx] = loadHistFromCCDB(timestamp); + break; + case 2: + hLoaded[idx] = loadHistFromCCDB(timestamp); + break; + case 3: + hLoaded[idx] = loadHistFromCCDB(timestamp); + break; + default: + LOGF(fatal, notify("Unknown configuration for efficiency variables")); + break; + } + } else { + hLoaded[idx] = nullptr; + } } } } - template - requires(sizeof...(BinVars) == getHistDim()) - auto getWeight(ParticleNo partNo, const BinVars&... binVars) const -> float + template + requires IsOneOrTwo + void fillTruthHist(auto particle) + { + if (!shouldFillHistograms) { + return; + } + + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMCTruth"), + particle.pt(), + particle.eta(), + particle.fdCollision().multV0M()); + } + + template + requires IsOneOrTwo + void fillRecoHist(auto particle, int particlePDG) + { + if (!shouldFillHistograms) { + return; + } + + if (!particle.has_fdMCParticle()) { + return; + } + + auto mcParticle = particle.fdMCParticle(); + + if (mcParticle.pdgMCTruth() == particlePDG) { + switch (mcParticle.partOriginMCTruth()) { + case (o2::aod::femtouniverse_mc_particle::kPrimary): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hPrimary"), + mcParticle.pt(), + mcParticle.eta(), + particle.fdCollision().multV0M()); + break; + + case (o2::aod::femtouniverse_mc_particle::kDaughter): + case (o2::aod::femtouniverse_mc_particle::kDaughterLambda): + case (o2::aod::femtouniverse_mc_particle::kDaughterSigmaplus): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hSecondary"), + mcParticle.pt(), + mcParticle.eta(), + particle.fdCollision().multV0M()); + break; + + case (o2::aod::femtouniverse_mc_particle::kMaterial): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hMaterial"), + mcParticle.pt(), + mcParticle.eta(), + particle.fdCollision().multV0M()); + break; + + case (o2::aod::femtouniverse_mc_particle::kFake): + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hFake"), + mcParticle.pt(), + mcParticle.eta(), + particle.fdCollision().multV0M()); + break; + + default: + histRegistry->fill(HIST(histDirectory) + HIST("/") + HIST(histSuffix[N - 1]) + HIST("/hOther"), + mcParticle.pt(), + mcParticle.eta(), + particle.fdCollision().multV0M()); + break; + } + } + } + + auto getWeight(ParticleNo partNo, auto particle) -> float { auto weight = 1.0f; auto hWeights = hLoaded[partNo - 1]; if (shouldApplyCorrection && hWeights) { - auto bin = hWeights->FindBin(binVars...); + auto dim = static_cast(hWeights->GetDimension()); + if (dim != getDimensionFromVariables()) { + LOGF(fatal, notify("Histogram \"%s\" has wrong dimension %d != %d"), config->confEffCorCCDBPath.value, dim, config->confEffCorVariables.value.size()); + return weight; + } + + auto bin = -1; + if (config->confEffCorVariables.value == "pt") { + bin = hLoaded[partNo - 1]->FindBin(particle.pt()); + } else if (config->confEffCorVariables.value == "pt,eta") { + bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta()); + } else if (config->confEffCorVariables.value == "pt,mult") { + bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.fdCollision().multV0M()); + } else if (config->confEffCorVariables.value == "pt,eta,mult") { + bin = hLoaded[partNo - 1]->FindBin(particle.pt(), particle.eta(), particle.fdCollision().multV0M()); + } else { + LOGF(fatal, notify("Unknown configuration for efficiency variables")); + return weight; + } weight = hWeights->GetBinContent(bin); } @@ -112,7 +224,7 @@ class EfficiencyCorrection return fmt::format("[EFFICIENCY CORRECTION] {}", msg); } - static auto isHistEmpty(HistType* hist) -> bool + static auto isHistEmpty(TH1* hist) -> bool { if (!hist) { return true; @@ -125,9 +237,10 @@ class EfficiencyCorrection return true; } - auto loadHistFromCCDB(const int64_t timestamp) const -> HistType* + template + auto loadHistFromCCDB(const int64_t timestamp) const -> H* { - auto hWeights = ccdb.getForTimeStamp(config->confEffCorCCDBPath, timestamp); + auto hWeights = ccdb.getForTimeStamp(config->confEffCorCCDBPath, timestamp); if (!hWeights || hWeights->IsZombie()) { LOGF(error, notify("Could not load histogram \"%s/%ld\""), config->confEffCorCCDBPath.value, timestamp); return nullptr; @@ -141,12 +254,23 @@ class EfficiencyCorrection return hWeights; } + auto getDimensionFromVariables() -> size_t + { + auto parts = std::views::split(config->confEffCorVariables.value, ','); + return std::ranges::distance(parts); + } + EffCorConfigurableGroup* config{}; - bool shouldApplyCorrection = false; + bool shouldApplyCorrection{false}; + bool shouldFillHistograms{false}; o2::ccdb::BasicCCDBManager& ccdb{o2::ccdb::BasicCCDBManager::instance()}; - std::array hLoaded{}; + std::array hLoaded{}; + + framework::HistogramRegistry* histRegistry{}; + static constexpr std::string_view histDirectory{"EfficiencyCorrection"}; + static constexpr std::string_view histSuffix[2]{"one", "two"}; }; } // namespace o2::analysis::femto_universe::efficiency_correction diff --git a/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx b/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx new file mode 100644 index 00000000000..e99331231e2 --- /dev/null +++ b/PWGCF/FemtoUniverse/Macros/calculateEfficiencyCorrection.cxx @@ -0,0 +1,240 @@ +// 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 calculateEfficiencyCorrection.cxx +/// \brief Macro for calculating efficiency corrections based on 3D histograms +/// \author Dawid Karpiński, WUT Warsaw, dawid.karpinski@cern.ch + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include // NOLINT +#include +#include +#include + +namespace fs = std::filesystem; + +auto* getHistogram(TFile* file, const std::string& name) +{ + return dynamic_cast(file->Get(name.c_str())); +} + +auto* projectHistogram(TH3* hist, const std::string& projection) +{ + return hist->Project3D(projection.c_str()); +} + +template +H* cloneHistogram(H* hist, const std::string& name) +{ + return dynamic_cast(hist->Clone(name.c_str())); +} + +auto forEachBin(TH1* hist, auto func) -> void +{ + if (hist->GetDimension() == 1) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + func(x, 0, 0); + } + } else if (hist->GetDimension() == 2) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + for (auto y{1}; y <= hist->GetNbinsY(); ++y) { + func(x, y, 0); + } + } + } else if (hist->GetDimension() == 3) { + for (auto x{1}; x <= hist->GetNbinsX(); ++x) { + for (auto y{1}; y <= hist->GetNbinsY(); ++y) { + for (auto z{1}; z <= hist->GetNbinsZ(); ++z) { + func(x, y, z); + } + } + } + } else { + assert(false && "should not happen"); + } +} + +auto setAxisTitles(TH1* hist, const std::string& projection) -> void +{ + auto* xAxis{hist->GetXaxis()}; + auto* yAxis{hist->GetYaxis()}; + auto* zAxis{hist->GetZaxis()}; + + xAxis->SetTitle("#it{p}_{T} (GeV/#it{c})"); + + if (hist->GetDimension() == 2) { + if (projection == "yx") { + yAxis->SetTitle("#it{#eta}"); + } else if (projection == "zx") { + yAxis->SetTitle("mult"); + } + } else if (hist->GetDimension() == 3) { + yAxis->SetTitle("#it{#eta}"); + zAxis->SetTitle("mult"); + } +} + +auto calculateEfficiencyCorrection(const fs::path& resultsPath, const fs::path& histPath, const std::string& projection) -> void +{ + assert(!resultsPath.empty() && !histPath.empty()); + if (projection != "" && projection != "x" && projection != "yx" && projection != "zx") { + std::cerr << "Error: projection must be one of: x, yx, zx\n"; + std::exit(1); + return; + } + + auto isAlien{false}; + if (resultsPath.string().starts_with("alien://")) { + TGrid::Connect("alien://"); + isAlien = true; + } + + auto* resultFile{TFile::Open(resultsPath.c_str())}; + assert(resultFile != nullptr && !resultFile->IsZombie()); + + using namespace std::chrono; + auto now{duration_cast(system_clock::now().time_since_epoch()).count()}; + auto outputPath{isAlien ? std::filesystem::current_path() : resultsPath.parent_path()}; + outputPath /= std::format("{}-effcor-{}.root", resultsPath.stem().string(), now); + + auto* outputFile{TFile::Open(outputPath.c_str(), "RECREATE")}; + assert(outputFile != nullptr && !outputFile->IsZombie()); + + auto* histTruthBase{getHistogram(resultFile, histPath / "hMCTruth")}; + auto* histPrimaryBase{getHistogram(resultFile, histPath / "hPrimary")}; + auto* histSecondaryBase{getHistogram(resultFile, histPath / "hSecondary")}; + + assert(histTruthBase); + assert(histPrimaryBase); + assert(histSecondaryBase); + + TH1* histTruth{histTruthBase}; + TH1* histPrimary{histPrimaryBase}; + TH1* histSecondary{histSecondaryBase}; + + if (projection != "") { + histPrimary = projectHistogram(histPrimaryBase, projection); + histSecondary = projectHistogram(histSecondaryBase, projection); + histTruth = projectHistogram(histTruthBase, projection); + } + + auto* histTotal{cloneHistogram(histPrimary, "hTotal")}; + histTotal->Add(histSecondary); + + auto* histEfficiency{cloneHistogram(histPrimary, "hEfficiency")}; + histEfficiency->Reset(); + setAxisTitles(histEfficiency, projection); + + auto* histWeights{cloneHistogram(histPrimary, "hWeights")}; + histWeights->Reset(); + setAxisTitles(histWeights, projection); + + forEachBin(histPrimary, [&](int x, int y, int z) { + auto primVal{histPrimary->GetBinContent(x, y, z)}; + auto primErr{histPrimary->GetBinError(x, y, z)}; + + auto secVal{histSecondary->GetBinContent(x, y, z)}; + auto secErr{histSecondary->GetBinError(x, y, z)}; + + auto truthVal{histTruth->GetBinContent(x, y, z)}; + auto truthErr{histTruth->GetBinError(x, y, z)}; + + auto effVal{0.}; + auto effErr{0.}; + if (truthVal > 0) { + effVal = primVal / truthVal; + effErr = std::sqrt(std::pow(primErr / truthVal, 2) + std::pow((primVal * truthErr / std::pow(truthVal, 2)), 2)); + } + + histEfficiency->SetBinContent(x, y, z, effVal); + histEfficiency->SetBinError(x, y, z, effErr); + + auto totalVal{primVal + secVal}; + auto totalErr{std::hypot(primErr, secErr)}; + + auto contVal{0.}; + auto contErr{0.}; + if (totalVal > 0) { + contVal = secVal / totalVal; + contErr = std::sqrt(std::pow(secErr / totalVal, 2) + std::pow((secVal * totalErr / std::pow(totalVal, 2)), 2)); + } + + auto weightVal{0.}; + auto weightErr{0.}; + if (effVal > 0) { + weightVal = (1 - contVal) / effVal; + weightErr = std::sqrt(std::pow(contErr / effVal, 2) + std::pow((1 - contVal) * effErr / std::pow(effVal, 2), 2)); + } + + histWeights->SetBinContent(x, y, z, weightVal); + histWeights->SetBinError(x, y, z, weightErr); + }); + + outputFile->WriteTObject(histEfficiency); + outputFile->WriteTObject(histWeights); + + outputFile->Close(); + resultFile->Close(); +} + +auto printUsage(const char* name) -> void +{ + std::cerr << "Usage: " << name << "\n" + << " -f \n" + << " -d \n" + << " -p [optional, default: no projection, 3D histogram]\n" + << " Available projections:\n" + << " x - projection onto pT axis (1D histogram)\n" + << " yx - projection onto pT, eta (2D histogram)\n" + << " zx - projection onto pT, mult (2D histogram)\n"; +} + +int main(int argc, char** argv) +{ + std::string results{""}, hist{""}, proj{""}; + auto flag{0}; + + while ((flag = getopt(argc, argv, "f:d:p:")) != -1) { + switch (flag) { + case 'f': + results = optarg; + break; + case 'd': + hist = optarg; + break; + case 'p': + proj = optarg; + break; + default: + printUsage(argv[0]); + return 1; + } + } + + if (results.empty() || hist.empty()) { + printUsage(argv[0]); + return 1; + } + + calculateEfficiencyCorrection(results, hist, proj); + return 0; +} diff --git a/PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_calculator.py b/PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_calculator.py similarity index 100% rename from PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_calculator.py rename to PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_calculator.py diff --git a/PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_phi_calculator.py b/PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_phi_calculator.py similarity index 100% rename from PWGCF/FemtoUniverse/Scripts/femto_universe_efficiency_phi_calculator.py rename to PWGCF/FemtoUniverse/Macros/femto_universe_efficiency_phi_calculator.py diff --git a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx index ce06e110b2c..55a19bf897c 100644 --- a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx +++ b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx @@ -708,10 +708,11 @@ struct FemtoUniverseProducerTask { particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kPrimary; } else if (!motherparticlesMC.empty()) { auto motherparticleMC = motherparticlesMC.front(); - if (motherparticleMC.producedByGenerator()) + if (motherparticleMC.producedByGenerator()) { particleOrigin = checkDaughterType(fdparttype, motherparticleMC.pdgCode()); - } else { - particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kMaterial; + } else { + particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kMaterial; + } } } else { particleOrigin = aod::femtouniverse_mc_particle::ParticleOriginMCTruth::kFake; @@ -902,8 +903,8 @@ struct FemtoUniverseProducerTask { { for (const auto& c : col) { const auto vtxZ = c.posZ(); - float mult = 0; - int multNtr = 0; + float mult = confIsRun3 ? c.multFV0M() : 0.5 * (c.multFV0M()); + int multNtr = confIsRun3 ? c.multNTracksPV() : c.multTracklets(); if (std::abs(vtxZ) > confEvtZvtx) { continue; diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx index 27e865ce552..b85e4832089 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackExtended.cxx @@ -105,7 +105,6 @@ struct FemtoUniversePairTaskTrackTrackExtended { /// Histogramming for particle 1 FemtoUniverseParticleHisto trackHistoPartOne; - FemtoUniverseParticleHisto hMCTruth1; /// Particle 2 Configurable confIsSame{"confIsSame", false, "Pairs of the same particle"}; @@ -129,7 +128,6 @@ struct FemtoUniversePairTaskTrackTrackExtended { /// Histogramming for particle 2 FemtoUniverseParticleHisto trackHistoPartTwo; - FemtoUniverseParticleHisto hMCTruth2; /// Histogramming for Event FemtoUniverseEventHisto eventHisto; @@ -176,9 +174,10 @@ struct FemtoUniversePairTaskTrackTrackExtended { HistogramRegistry qaRegistry{"TrackQA", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry resultRegistry{"Correlations", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry mixQaRegistry{"mixQaRegistry", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry effCorrRegistry{"EfficiencyCorrection", {}, OutputObjHandlingPolicy::AnalysisObject}; EffCorConfigurableGroup effCorConfGroup; - EfficiencyCorrection effCorrection{&effCorConfGroup}; + EfficiencyCorrection effCorrection{&effCorConfGroup}; /// @brief Counter for particle swapping int fNeventsProcessed = 0; @@ -322,18 +321,10 @@ struct FemtoUniversePairTaskTrackTrackExtended { void init(InitContext&) { - if (twotracksconfigs.confIsMC) { - hMCTruth1.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarPDGBins, false, trackonefilter.confPDGCodePartOne, false); - if (!confIsSame) { - hMCTruth2.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarPDGBins, false, tracktwofilter.confPDGCodePartTwo, false); - } - } - effCorrection.init(); - eventHisto.init(&qaRegistry); - trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true); // last true = isDebug + trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true, std::nullopt); // last true = isDebug if (!confIsSame) { - trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true); // last true = isDebug + trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true, std::nullopt); // last true = isDebug } mixQaRegistry.add("MixingQA/hSECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); @@ -351,6 +342,14 @@ struct FemtoUniversePairTaskTrackTrackExtended { vPIDPartOne = trackonefilter.confPIDPartOne.value; vPIDPartTwo = tracktwofilter.confPIDPartTwo.value; kNsigma = twotracksconfigs.confTrkPIDnSigmaMax.value; + + effCorrection.init( + &effCorrRegistry, + { + static_cast(confTempFitVarpTBins), + {confEtaBins, -2, 2}, + confMultBins, + }); } template @@ -361,23 +360,12 @@ struct FemtoUniversePairTaskTrackTrackExtended { } template - requires isOneOrTwo - auto doMCTruth(FemtoUniverseParticleHisto hist, auto parts) -> void + requires IsOneOrTwo + auto doMCTruth(auto parts, int partPDG, int partCharge) -> void { - auto expectedPDG = 0; - auto expectedCharge = 0.0l; - - if constexpr (N == ParticleNo::ONE) { - expectedPDG = trackonefilter.confPDGCodePartOne; - expectedCharge = trackonefilter.confChargePart1; - } else if constexpr (N == ParticleNo::TWO) { - expectedPDG = tracktwofilter.confPDGCodePartTwo; - expectedCharge = tracktwofilter.confChargePart2; - } - for (const auto& particle : parts) { auto pdgCode = static_cast(particle.pidCut()); - if (pdgCode != expectedPDG) { + if (pdgCode != partPDG) { continue; } @@ -386,9 +374,11 @@ struct FemtoUniversePairTaskTrackTrackExtended { continue; } - if (sign(pdgParticle->Charge()) == sign(expectedCharge)) { - hist.template fillQA(particle); + if (sign(pdgParticle->Charge()) != sign(partCharge)) { + continue; } + + effCorrection.fillTruthHist(particle); } } @@ -432,6 +422,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { } trackHistoPartOne.fillQA(part); + if constexpr (isMC) { + effCorrection.fillRecoHist(part, trackonefilter.confPDGCodePartOne); + } } if (!confIsSame) { @@ -456,6 +449,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { } } trackHistoPartTwo.fillQA(part); + if constexpr (isMC) { + effCorrection.fillRecoHist(part, tracktwofilter.confPDGCodePartTwo); + } } /// Now build the combinations for non-identical particle pairs @@ -505,9 +501,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { continue; } - float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt()); + float weight = effCorrection.getWeight(ParticleNo::ONE, p1); if (!confIsSame) { - weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt()); + weight *= effCorrection.getWeight(ParticleNo::TWO, p2); } if (swpart) @@ -566,9 +562,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { continue; } - float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt()); + float weight = effCorrection.getWeight(ParticleNo::ONE, p1); if (!confIsSame) { - weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt()); + weight *= effCorrection.getWeight(ParticleNo::TWO, p2); } sameEventCont.setPair(p1, p2, multCol, twotracksconfigs.confUse3D, weight); @@ -602,11 +598,11 @@ struct FemtoUniversePairTaskTrackTrackExtended { fillCollision(col); auto groupMCTruth1 = partsOneMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); - doMCTruth<1>(hMCTruth1, groupMCTruth1); + doMCTruth(groupMCTruth1, trackonefilter.confPDGCodePartOne, trackonefilter.confChargePart1); if (!confIsSame) { auto groupMCTruth2 = partsTwoMCTruth->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); - doMCTruth<2>(hMCTruth2, groupMCTruth2); + doMCTruth(groupMCTruth2, tracktwofilter.confPDGCodePartTwo, tracktwofilter.confChargePart2); } auto groupMCReco1 = partsOneMCReco->sliceByCached(aod::femtouniverseparticle::fdCollisionId, col.globalIndex(), cache); @@ -674,9 +670,9 @@ struct FemtoUniversePairTaskTrackTrackExtended { } } - float weight = effCorrection.getWeight(ParticleNo::ONE, p1.pt()); + float weight = effCorrection.getWeight(ParticleNo::ONE, p1); if (!confIsSame) { - weight *= effCorrection.getWeight(ParticleNo::TWO, p2.pt()); + weight *= effCorrection.getWeight(ParticleNo::TWO, p2); } if (swpart)