diff --git a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx index 86ce8a1a340..c37bdc6aeb9 100644 --- a/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx +++ b/PWGLF/Tasks/Nuspex/antinucleiInJets.cxx @@ -142,6 +142,22 @@ struct AntinucleiInJets { Configurable cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {0.6539, 1.591, 0.8225, 2.363, 0.09}, "TPC Bethe-Bloch parameterisation for He3"}; + // Configuration parameters for CCDB access and reweighting input files + Configurable urlToCcdb{"urlToCcdb", "http://alice-ccdb.cern.ch", "url of the personal ccdb"}; + Configurable pathToFile{"pathToFile", "Users/a/alcaliva/PrimaryFractionAntip", "path to file with reweighting"}; + Configurable weightsProton{"weightsProton", "", "weightsProton"}; + Configurable weightsLambda{"weightsLambda", "", "weightsLambda"}; + Configurable weightsSigma{"weightsSigma", "", "weightsSigma"}; + Configurable weightsXi{"weightsXi", "", "weightsXi"}; + Configurable weightsOmega{"weightsOmega", "", "weightsOmega"}; + + // Reweighting histograms + TH1F* primaryAntiprotons; + TH1F* primaryAntiLambda; + TH1F* primaryAntiSigma; + TH1F* primaryAntiXi; + TH1F* primaryAntiOmega; + // CCDB manager service for accessing condition data Service ccdb; @@ -179,6 +195,16 @@ struct AntinucleiInJets { itsResponse.setMCDefaultParameters(); } + // Load reweighting histograms from CCDB if antinuclei efficiency processing is enabled + if (doprocessAntinucleiEfficiency) { + ccdb->setURL(urlToCcdb.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + ccdb->setFatalWhenNull(false); + getReweightingHistograms(ccdb, TString(pathToFile), TString(weightsProton), TString(weightsLambda), TString(weightsSigma), TString(weightsXi), TString(weightsOmega)); + } + // Binning double min = 0.0; double max = 6.0; @@ -361,6 +387,17 @@ struct AntinucleiInJets { // nsigmaITS for antiproton candidates registryMC.add("antiproton_nsigma_its_mc", "antiproton_nsigma_its_mc", HistType::kTH2F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}, {400, -20.0, 20.0, "n#sigma_{ITS}"}}); + + // Systematics on the fraction of primary antiprotons + registryMC.add("antip_prim_pythia", "antip_prim_pythia", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("antip_prim_std", "antip_prim_std", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("antip_prim_up", "antip_prim_up", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("antip_prim_low", "antip_prim_low", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + + registryMC.add("antip_sec_pythia", "antip_sec_pythia", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("antip_sec_std", "antip_sec_std", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("antip_sec_up", "antip_sec_up", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); + registryMC.add("antip_sec_low", "antip_sec_low", HistType::kTH1F, {{nbins, min, max, "#it{p}_{T} (GeV/#it{c})"}}); } // Systematic uncertainties (Data) @@ -442,6 +479,27 @@ struct AntinucleiInJets { } } + void getReweightingHistograms(o2::framework::Service const& ccdbObj, TString filepath, TString antip, TString antilambda, TString antisigma, TString antixi, TString antiomega) + { + TList* list = ccdbObj->get(filepath.Data()); + if (!list) { + LOGP(error, "Could not retrieve the list from CCDB"); + return; + } + + primaryAntiprotons = static_cast(list->FindObject(antip)); + primaryAntiLambda = static_cast(list->FindObject(antilambda)); + primaryAntiSigma = static_cast(list->FindObject(antisigma)); + primaryAntiXi = static_cast(list->FindObject(antixi)); + primaryAntiOmega = static_cast(list->FindObject(antiomega)); + + if (!primaryAntiprotons || !primaryAntiSigma || !primaryAntiLambda || !primaryAntiXi || !primaryAntiOmega) { + LOGP(error, "Missing one or more reweighting histograms in CCDB list"); + } + + LOGP(info, "Successfully loaded reweighting histograms from CCDB path"); + } + // Compute two unit vectors perpendicular to p void getPerpendicularAxis(const TVector3& p, TVector3& u, double sign) { @@ -1344,6 +1402,94 @@ struct AntinucleiInJets { continue; const auto particle = track.mcParticle(); + // **************************************************************************************************************** + + // Systematic uncertainty on primary fraction + if (track.sign() < 0 && particle.pdgCode() == PDG_t::kProtonBar) { + + // Primary antiprotons + if (particle.isPhysicalPrimary()) { + + // Initialize weights + double wPrimStd(1.0), wPrimUp(1.0), wPrimLow(1.0); + + // Weight assignment + if (particle.pt() < primaryAntiprotons->GetXaxis()->GetXmax()) { + int ipt = primaryAntiprotons->FindBin(particle.pt()); + wPrimStd = primaryAntiprotons->GetBinContent(ipt); + wPrimUp = wPrimStd + primaryAntiprotons->GetBinError(ipt); + wPrimLow = wPrimStd - primaryAntiprotons->GetBinError(ipt); + } + + // Fill histograms + registryMC.fill(HIST("antip_prim_pythia"), track.pt()); + registryMC.fill(HIST("antip_prim_std"), track.pt(), wPrimStd); + registryMC.fill(HIST("antip_prim_up"), track.pt(), wPrimUp); + registryMC.fill(HIST("antip_prim_low"), track.pt(), wPrimLow); + } + + // Secondary antiprotons from week decays + if (!particle.isPhysicalPrimary() && particle.has_mothers()) { + + // Initialize weights + double wSecStd(1.0), wSecUp(1.0), wSecLow(1.0); + auto mother = mcParticles.iteratorAt(particle.mothersIds()[0]); + + // Antiprotons from sigma + if (std::abs(mother.pdgCode()) == PDG_t::kSigmaBarMinus) { + if (mother.pt() < primaryAntiSigma->GetXaxis()->GetXmax()) { + int ipt = primaryAntiSigma->FindBin(mother.pt()); + wSecStd = primaryAntiSigma->GetBinContent(ipt); + wSecUp = wSecStd + primaryAntiSigma->GetBinError(ipt); + wSecLow = wSecStd - primaryAntiSigma->GetBinError(ipt); + } + } + + // Antiprotons from primary Lambda0 + if (std::abs(mother.pdgCode()) == kLambda0Bar) { + if (mother.isPhysicalPrimary()) { + if (mother.pt() < primaryAntiLambda->GetXaxis()->GetXmax()) { + int ipt = primaryAntiLambda->FindBin(mother.pt()); + wSecStd = primaryAntiLambda->GetBinContent(ipt); + wSecUp = wSecStd + primaryAntiLambda->GetBinError(ipt); + wSecLow = wSecStd - primaryAntiLambda->GetBinError(ipt); + } + } + + // Antiprotons from secondary Lambda0 (Xi -> Lambda0) + if (!mother.isPhysicalPrimary()) { + auto grandmother = mcParticles.iteratorAt(mother.mothersIds()[0]); + if (std::abs(grandmother.pdgCode()) == kXiMinus) { + if (grandmother.pt() < primaryAntiXi->GetXaxis()->GetXmax()) { + int ipt = primaryAntiXi->FindBin(grandmother.pt()); + wSecStd = primaryAntiXi->GetBinContent(ipt); + wSecUp = wSecStd + primaryAntiXi->GetBinError(ipt); + wSecLow = wSecStd - primaryAntiXi->GetBinError(ipt); + } + } + + // Antiprotons from secondary Lambda0 (Omega -> Lambda0) + if (std::abs(grandmother.pdgCode()) == kOmegaMinus) { + if (grandmother.pt() < primaryAntiOmega->GetXaxis()->GetXmax()) { + int ipt = primaryAntiOmega->FindBin(grandmother.pt()); + wSecStd = primaryAntiOmega->GetBinContent(ipt); + wSecUp = wSecStd + primaryAntiOmega->GetBinError(ipt); + wSecLow = wSecStd - primaryAntiOmega->GetBinError(ipt); + } + } + } + } + + // Fill histograms + registryMC.fill(HIST("antip_sec_pythia"), track.pt()); + registryMC.fill(HIST("antip_sec_std"), track.pt(), wSecStd); + registryMC.fill(HIST("antip_sec_up"), track.pt(), wSecUp); + registryMC.fill(HIST("antip_sec_low"), track.pt(), wSecLow); + } + } + + // **************************************************************************************************************** + // Select only physical primary particles if (!particle.isPhysicalPrimary()) continue;