diff --git a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx index e848d353761..9c7f14a842e 100644 --- a/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx +++ b/PWGCF/GenericFramework/Tasks/flowGenericFramework.cxx @@ -137,11 +137,13 @@ struct FlowGenericFramework { O2_DEFINE_CONFIGURABLE(cfgUsePID, bool, true, "Enable PID information") O2_DEFINE_CONFIGURABLE(cfgUseGapMethod, bool, false, "Use gap method in vn-pt calculations") O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object") + O2_DEFINE_CONFIGURABLE(cfgUsePIDEfficiencies, bool, false, "Use species dependent efficiencies") O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object") O2_DEFINE_CONFIGURABLE(cfgPtmin, float, 0.2, "minimum pt (GeV/c)"); O2_DEFINE_CONFIGURABLE(cfgPtmax, float, 10, "maximum pt (GeV/c)"); O2_DEFINE_CONFIGURABLE(cfgEta, float, 0.8, "eta cut"); O2_DEFINE_CONFIGURABLE(cfgEtaPtPt, float, 0.4, "eta cut for pt-pt correlations"); + O2_DEFINE_CONFIGURABLE(cfgEtaNch, float, 0.4, "eta cut for nch selection"); O2_DEFINE_CONFIGURABLE(cfgUsePIDTotal, bool, false, "use fraction of PID total"); O2_DEFINE_CONFIGURABLE(cfgVtxZ, float, 10, "vertex cut (cm)"); struct : ConfigurableGroup { @@ -223,6 +225,7 @@ struct FlowGenericFramework { struct Config { TH1D* mEfficiency = nullptr; + std::vector mPIDEfficiencies; std::vector mAcceptance; bool correctionsLoaded = false; } cfg; @@ -845,11 +848,21 @@ struct FlowGenericFramework { } } if (!cfgEfficiency.value.empty()) { - cfg.mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); - if (cfg.mEfficiency == nullptr) { - LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str()); + if (!cfgUsePIDEfficiencies) { + cfg.mEfficiency = ccdb->getForTimeStamp(cfgEfficiency, timestamp); + if (cfg.mEfficiency == nullptr) { + LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str()); + } + LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency); + } else { + std::vector species = {"ch", "pi", "ka", "pr"}; + for (const auto& sp : species) { + cfg.mPIDEfficiencies.push_back(ccdb->getForTimeStamp(cfgEfficiency.value + "/" + sp, timestamp)); + if (cfg.mPIDEfficiencies.back() == nullptr) + LOGF(fatal, "Could not load PID efficiency histograms from %s", cfgEfficiency.value + "/" + sp); + LOGF(info, "Loaded PID efficiency histogram from %s", cfgEfficiency.value + "/" + sp); + } } - LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency); } cfg.correctionsLoaded = true; } @@ -864,11 +877,15 @@ struct FlowGenericFramework { } template - double getEfficiency(TTrack track) + double getEfficiency(TTrack track, int pidIndex = 0) { //-1 ref, 0 ch, 1 pi, 2 ka, 3 pr double eff = 1.; - if (cfg.mEfficiency) - eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt())); + if (!cfgUsePIDEfficiencies) { + if (cfg.mEfficiency) + eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt())); + } else { + eff = cfg.mPIDEfficiencies[pidIndex]->GetBinContent(cfg.mPIDEfficiencies[pidIndex]->FindBin(track.pt())); + } if (eff == 0) return -1.; else @@ -1231,7 +1248,7 @@ struct FlowGenericFramework { float total = chtotal; if (cfgUsePIDTotal) - (pidcounter) ? acceptedtracks.pidtotal[pidcounter] : chtotal; + total = (pidcounter) ? acceptedtracks.pidtotal[pidcounter] : chtotal; if (total == 0.) { ++pidcounter; @@ -1336,11 +1353,11 @@ struct FlowGenericFramework { int ptBinIndex = fPtAxis->FindBin(v0.pt()) - 1; if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { if (v0.mK0Short() < cfgPIDCuts.cfgK0SignalMin) - npt_resonances[0][ptBinIndex] += getEfficiency(v0); + npt_resonances[0][ptBinIndex]++; else if (v0.mK0Short() > cfgPIDCuts.cfgK0SignalMax) - npt_resonances[2][ptBinIndex] += getEfficiency(v0); + npt_resonances[2][ptBinIndex]++; else - npt_resonances[1][ptBinIndex] += getEfficiency(v0); + npt_resonances[1][ptBinIndex]++; registry.fill(HIST("K0/hK0s"), 1); } } @@ -1351,11 +1368,11 @@ struct FlowGenericFramework { int ptBinIndex = fPtAxis->FindBin(v0.pt()) - 1; if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { if (v0.mLambda() < cfgPIDCuts.cfgLambdaSignalMin) - npt_resonances[3][ptBinIndex] += getEfficiency(v0); + npt_resonances[3][ptBinIndex]++; else if (v0.mLambda() > cfgPIDCuts.cfgLambdaSignalMax) - npt_resonances[5][ptBinIndex] += getEfficiency(v0); + npt_resonances[5][ptBinIndex]++; else - npt_resonances[4][ptBinIndex] += getEfficiency(v0); + npt_resonances[4][ptBinIndex]++; registry.fill(HIST("Lambda/hLambdas"), 1); } } @@ -1368,7 +1385,7 @@ struct FlowGenericFramework { for (const auto& vec : fractions_resonances) { float total = chtotal; if (cfgUsePIDTotal) - (pidcounter) ? std::accumulate(vec.begin(), vec.end(), 0) : chtotal; + total = (pidcounter) ? std::accumulate(vec.begin(), vec.end(), 0) : chtotal; if (total == 0.) { ++pidcounter; @@ -1441,8 +1458,12 @@ struct FlowGenericFramework { if (!nchSelected(track)) return; - acceptedTracks.total += getEfficiency(track); - ++acceptedTracks.total_uncorr; + double weffCh = getEfficiency(track, 0); + if (std::abs(track.eta()) < cfgEtaNch) { + if (weffCh > 0) + acceptedTracks.total += weffCh; + ++acceptedTracks.total_uncorr; + } if (!trackSelected(track, field)) return; @@ -1457,18 +1478,23 @@ struct FlowGenericFramework { pidIndex = 3; } - if (pidIndex) - acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track); - int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; - - if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { - acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 1) - acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 2) - acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 3) - acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + if (std::abs(track.eta()) < cfgEtaNch) { + double weff = getEfficiency(track, pidIndex); + + if (pidIndex && weff > 0) + acceptedTracks.pidtotal[pidIndex - 1] += weff; + + int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + if (weffCh > 0) + acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? weffCh : 1.0; + if (pidIndex == 1 && weff > 0) + acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0; + if (pidIndex == 2 && weff > 0) + acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0; + if (pidIndex == 3 && weff > 0) + acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0; + } } if (cfgFillWeights) { @@ -1504,21 +1530,24 @@ struct FlowGenericFramework { if (std::abs(track.pdgCode()) == kProton) pidIndex = 3; } - ++acceptedTracks.total; - ++acceptedTracks.total_uncorr; - if (pidIndex) - acceptedTracks.pidtotal[pidIndex - 1] += 1; - int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; - - if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { - acceptedTracks.nch[ptBinIndex] += 1.0; - if (pidIndex == 1) - acceptedTracks.npi[ptBinIndex] += 1.0; - if (pidIndex == 2) - acceptedTracks.nka[ptBinIndex] += 1.0; - if (pidIndex == 3) - acceptedTracks.npr[ptBinIndex] += 1.0; + if (std::abs(track.eta()) < cfgEtaNch) { + ++acceptedTracks.total; + ++acceptedTracks.total_uncorr; + + if (pidIndex) + acceptedTracks.pidtotal[pidIndex - 1] += 1; + int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; + + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + acceptedTracks.nch[ptBinIndex] += 1.0; + if (pidIndex == 1) + acceptedTracks.npi[ptBinIndex] += 1.0; + if (pidIndex == 2) + acceptedTracks.nka[ptBinIndex] += 1.0; + if (pidIndex == 3) + acceptedTracks.npr[ptBinIndex] += 1.0; + } } fillPtSums(track, vtxz, pidIndex); @@ -1532,8 +1561,13 @@ struct FlowGenericFramework { // Select tracks with nominal cuts always if (!nchSelected(track)) return; - acceptedTracks.total += getEfficiency(track); - ++acceptedTracks.total_uncorr; + + double weffCh = getEfficiency(track, 0); + if (std::abs(track.eta()) < cfgEtaNch) { + if (weffCh > 0) + acceptedTracks.total += weffCh; + ++acceptedTracks.total_uncorr; + } if (!trackSelected(track, field)) return; @@ -1541,23 +1575,25 @@ struct FlowGenericFramework { // if (cfgUsePID) Need PID for v02 int pidIndex = getNsigmaPID(track); - if (pidIndex) - acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track); + if (std::abs(track.eta()) < cfgEtaNch) { + double weff = getEfficiency(track, pidIndex); + if (pidIndex && weff > 0) + acceptedTracks.pidtotal[pidIndex - 1] += weff; - if (pidIndex) - acceptedTracks.pidtotal[pidIndex - 1] += getEfficiency(track); - int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; + int ptBinIndex = fPtAxis->FindBin(track.pt()) - 1; - if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { - acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - if (pidIndex == 1) { - acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - } - if (pidIndex == 2) { - acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; - } - if (pidIndex == 3) { - acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? getEfficiency(track) : 1.0; + if (!(ptBinIndex < 0 || ptBinIndex >= static_cast(o2::analysis::gfw::ptbinning.size()))) { + if (weffCh > 0) + acceptedTracks.nch[ptBinIndex] += (cfgUseNchCorrection) ? weffCh : 1.0; + if (pidIndex == 1 && weff > 0) { + acceptedTracks.npi[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0; + } + if (pidIndex == 2 && weff > 0) { + acceptedTracks.nka[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0; + } + if (pidIndex == 3 && weff > 0) { + acceptedTracks.npr[ptBinIndex] += (cfgUseNchCorrection) ? weff : 1.0; + } } } @@ -1779,27 +1815,31 @@ struct FlowGenericFramework { bool withinPtNch = (track.pt() > ptmins[0] && track.pt() < ptmaxs[0]); if (!withinPtPOI && !withinPtRef) return; + double weff = (dt == kGen) ? 1. : getEfficiency(track, pid_index); + double weffInclusive = (dt == kGen) ? 1. : getEfficiency(track, 0); + if (weff < 0) + return; double waccRef = (dt == kGen) ? 1. : getAcceptance(track, vtxz, 0); double waccPOI = (dt == kGen) ? 1. : withinPtPOI ? getAcceptance(track, vtxz, pid_index + 1) : getAcceptance(track, vtxz, 0); // if (withinPtRef && withinPtPOI && pid_index) waccRef = waccPOI; // if particle is both (then it's overlap), override ref with POI if (withinPtRef) - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccRef, 1); + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccRef * weffInclusive, 1); if (withinPtPOI && pid_index) - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, (1 << (pid_index + 1))); + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff, (1 << (pid_index + 1))); if (withinPtNch) - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, 2); + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff * weffInclusive, 2); if (withinPtPOI && withinPtRef && pid_index) - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, (1 << (pid_index + 5))); + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff, (1 << (pid_index + 5))); if (withinPtNch && withinPtRef) - fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI, 32); + fGFW->Fill(track.eta(), fPtAxis->FindBin(track.pt()) - 1, track.phi(), waccPOI * weff * weffInclusive, 32); } else { // Analysing only integrated flow bool withinPtRef = (track.pt() > o2::analysis::gfw::ptreflow && track.pt() < o2::analysis::gfw::ptrefup); bool withinPtPOI = (track.pt() > o2::analysis::gfw::ptpoilow && track.pt() < o2::analysis::gfw::ptpoiup); if (!withinPtPOI && !withinPtRef) return; - double weff = (dt == kGen) ? 1. : getEfficiency(track); + double weff = (dt == kGen) ? 1. : getEfficiency(track, 0); if (weff < 0) return; if (cfgUseDensityDependentCorrection && withinPtRef && dt != kGen) { @@ -1886,7 +1926,7 @@ struct FlowGenericFramework { } o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ; - o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz; + o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true)) && (aod::track::itsChi2NCl < cfgTrackCuts.cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgTrackCuts.cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgTrackCuts.cfgDCAz; using GFWCollisions = soa::Filtered>; // using GFWTracks = soa::Filtered>;