Skip to content

Commit 99cd23c

Browse files
RrantuRrantu
committed
add UPC infos in taskLc.cxx
1 parent c585c36 commit 99cd23c

File tree

1 file changed

+183
-35
lines changed

1 file changed

+183
-35
lines changed

‎PWGHF/D2H/Tasks/taskLc.cxx‎

Lines changed: 183 additions & 35 deletions
Original file line numberDiff line numberDiff line change
@@ -29,6 +29,7 @@
2929
#include "PWGHF/DataModel/CandidateSelectionTables.h"
3030
#include "PWGHF/DataModel/TrackIndexSkimmingTables.h"
3131
#include "PWGHF/Utils/utilsEvSelHf.h"
32+
#include "PWGHF/Utils/utilsUpcHf.h"
3233
#include "PWGUD/Core/UPCHelpers.h"
3334

3435
#include "Common/Core/RecoDecay.h"
@@ -67,6 +68,7 @@ using namespace o2::framework::expressions;
6768
using namespace o2::hf_centrality;
6869
using namespace o2::hf_occupancy;
6970
using namespace o2::hf_evsel;
71+
using namespace o2::analysis::hf_upc;
7072

7173
/// Λc± → p± K∓ π± analysis task
7274
struct HfTaskLc {
@@ -84,7 +86,8 @@ struct HfTaskLc {
8486
Configurable<std::string> ccdbPathGrp{"ccdbPathGrp", "GLO/GRP/GRP", "Path of the grp file (Run 2)"};
8587
Configurable<std::string> ccdbPathGrpMag{"ccdbPathGrpMag", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object (Run 3)"};
8688

87-
HfEventSelection hfEvSel; // event selection and monitoring
89+
HfEventSelection hfEvSel; // event selection and monitoring
90+
HfUpcGapThresholds upcThresholds; // UPC gap determination thresholds
8891
SliceCache cache;
8992
Service<o2::ccdb::BasicCCDBManager> ccdb;
9093

@@ -121,7 +124,11 @@ struct HfTaskLc {
121124
ConfigurableAxis thnConfigAxisNumPvContr{"thnConfigAxisNumPvContr", {200, -0.5, 199.5}, "Number of PV contributors"};
122125
ConfigurableAxis thnConfigAxisOccupancy{"thnConfigAxisOccupancy", {14, 0, 14000}, "axis for centrality"};
123126
ConfigurableAxis thnConfigAxisProperLifetime{"thnConfigAxisProperLifetime", {200, 0, 2}, "Proper lifetime, ps"};
124-
127+
ConfigurableAxis thnConfigAxisGapType{"thnConfigAxisGapType", {7, -1.5, 5.5}, "axis for UPC gap type (see TrueGap enum in o2::aod::sgselector)"};
128+
ConfigurableAxis thnConfigAxisFT0{"thnConfigAxisFT0", {1001, -1.5, 999.5}, "axis for FT0 amplitude (a.u.)"};
129+
ConfigurableAxis thnConfigAxisFV0A{"thnConfigAxisFV0A", {2001, -1.5, 1999.5}, "axis for FV0-A amplitude (a.u.)"};
130+
ConfigurableAxis thnConfigAxisFDD{"thnConfigAxisFDD", {200, 0., 4000.}, "axis for FDD amplitude (a.u.)"};
131+
ConfigurableAxis thnConfigAxisZN{"thnConfigAxisZN", {510, -1.5, 49.5}, "axis for ZN energy (a.u.)"};
125132
HistogramRegistry registry{"registry", {}};
126133
HistogramRegistry qaRegistry{"QAHistos", {}, OutputObjHandlingPolicy::AnalysisObject};
127134

@@ -145,12 +152,6 @@ struct HfTaskLc {
145152
NonPrompt
146153
};
147154

148-
enum class GapType {
149-
GapA = 0,
150-
GapC = 1,
151-
DoubleGap = 2,
152-
};
153-
154155
void init(InitContext&)
155156
{
156157
const std::array<bool, 14> doprocess{doprocessDataStd, doprocessDataStdWithFT0C, doprocessDataStdWithFT0M, doprocessDataWithMl, doprocessDataWithMlWithFT0C, doprocessDataWithMlWithFT0M, doprocessDataWithMlWithUpc, doprocessMcStd, doprocessMcStdWithFT0C, doprocessMcStdWithFT0M, doprocessMcWithMl, doprocessMcWithMlWithFT0C, doprocessMcWithMlWithFT0M, doprocessDataStdWithUpc};
@@ -265,10 +266,7 @@ struct HfTaskLc {
265266
if (isUpc) {
266267
qaRegistry.add("Data/fitInfo/ampFT0A_vs_ampFT0C", "FT0-A vs FT0-C amplitude;FT0-A amplitude (a.u.);FT0-C amplitude (a.u.)", {HistType::kTH2F, {{1500, 0., 1500}, {1500, 0., 1500}}});
267268
qaRegistry.add("Data/zdc/energyZNA_vs_energyZNC", "ZNA vs ZNC common energy;E_{ZNA}^{common} (a.u.);E_{ZNC}^{common} (a.u.)", {HistType::kTH2F, {{200, 0., 20}, {200, 0., 20}}});
268-
qaRegistry.add("Data/hUpcGapAfterSelection", "UPC gap type after selection;Gap side;Counts", {HistType::kTH1F, {{3, -0.5, 2.5}}});
269-
qaRegistry.get<TH1>(HIST("Data/hUpcGapAfterSelection"))->GetXaxis()->SetBinLabel(static_cast<int>(GapType::GapA) + 1, "A");
270-
qaRegistry.get<TH1>(HIST("Data/hUpcGapAfterSelection"))->GetXaxis()->SetBinLabel(static_cast<int>(GapType::GapC) + 1, "C");
271-
qaRegistry.get<TH1>(HIST("Data/hUpcGapAfterSelection"))->GetXaxis()->SetBinLabel(static_cast<int>(GapType::DoubleGap) + 1, "Double");
269+
qaRegistry.add("Data/hUpcGapAfterSelection", "UPC gap type after selection;Gap side;Counts", {HistType::kTH1F, {{7, -1.5, 5.5}}});
272270
}
273271
if (fillTHn) {
274272
const AxisSpec thnAxisMass{thnConfigAxisMass, "inv. mass (p K #pi) (GeV/#it{c}^{2})"};
@@ -289,6 +287,14 @@ struct HfTaskLc {
289287
const AxisSpec thnAxisTracklets{thnConfigAxisNumPvContr, "Number of PV contributors"};
290288
const AxisSpec thnAxisOccupancy{thnConfigAxisOccupancy, "Occupancy"};
291289
const AxisSpec thnAxisProperLifetime{thnConfigAxisProperLifetime, "T_{proper} (ps)"};
290+
const AxisSpec thnAxisGapType{thnConfigAxisGapType, "Gap type"};
291+
const AxisSpec thnAxisFT0A{thnConfigAxisFT0, "FT0-A amplitude"};
292+
const AxisSpec thnAxisFT0C{thnConfigAxisFT0, "FT0-C amplitude"};
293+
const AxisSpec thnAxisFV0A{thnConfigAxisFV0A, "FV0-A amplitude"};
294+
const AxisSpec thnAxisFDDA{thnConfigAxisFDD, "FDD-A amplitude"};
295+
const AxisSpec thnAxisFDDC{thnConfigAxisFDD, "FDD-C amplitude"};
296+
const AxisSpec thnAxisZNA{thnConfigAxisZN, "ZNA energy"};
297+
const AxisSpec thnAxisZNC{thnConfigAxisZN, "ZNC energy"};
292298

293299
bool const isDataWithMl = doprocessDataWithMl || doprocessDataWithMlWithFT0C || doprocessDataWithMlWithFT0M || doprocessDataWithMlWithUpc;
294300
bool const isMcWithMl = doprocessMcWithMl || doprocessMcWithMlWithFT0C || doprocessMcWithMlWithFT0M;
@@ -327,7 +333,20 @@ struct HfTaskLc {
327333
}
328334
}
329335
}
330-
336+
if (isUpc) {
337+
for (const auto& axes : std::array<std::vector<AxisSpec>*, 3>{&axesWithBdt, &axesStd, &axesGen}) {
338+
if (!axes->empty()) {
339+
axes->push_back(thnAxisGapType);
340+
axes->push_back(thnAxisFT0A);
341+
axes->push_back(thnAxisFT0C);
342+
axes->push_back(thnAxisFV0A);
343+
axes->push_back(thnAxisFDDA);
344+
axes->push_back(thnAxisFDDC);
345+
axes->push_back(thnAxisZNA);
346+
axes->push_back(thnAxisZNC);
347+
}
348+
}
349+
}
331350
if (isDataWithMl) {
332351
registry.add("hnLcVarsWithBdt", "THn for Lambdac candidates with BDT scores for data with ML", HistType::kTHnSparseF, axesWithBdt);
333352
} else if (isMcWithMl) {
@@ -727,36 +746,165 @@ struct HfTaskLc {
727746
/// at least one event selection not satisfied --> reject the candidate
728747
continue;
729748
}
749+
const auto thisCollId = collision.globalIndex();
750+
const auto& groupedLcCandidates = candidates.sliceBy(candLcPerCollision, thisCollId);
751+
const auto numPvContributors = collision.numContrib();
730752
const auto& bc = collision.template bc_as<BCsType>();
753+
754+
// Determine gap type using SGSelector with BC range checking
755+
const auto gapResult = hf_upc::determineGapType(collision, bcs, upcThresholds);
756+
const int gap = gapResult.value;
757+
758+
// Use the BC with FIT activity if available from SGSelector
759+
auto bcForUPC = bc;
760+
if (gapResult.bc) {
761+
bcForUPC = *(gapResult.bc);
762+
}
763+
764+
// Get FIT information from the UPC BC
731765
upchelpers::FITInfo fitInfo{};
732-
udhelpers::getFITinfo(fitInfo, bc, bcs, ft0s, fv0as, fdds);
766+
udhelpers::getFITinfo(fitInfo, bcForUPC, bcs, ft0s, fv0as, fdds);
733767

734-
GapType gap = GapType::DoubleGap;
735-
if (bc.has_zdc()) {
736-
const auto zdc = bc.zdc();
768+
// Get ZDC energies if available (extract once and reuse)
769+
const bool hasZdc = bcForUPC.has_zdc();
770+
float zdcEnergyZNA = -1.f;
771+
float zdcEnergyZNC = -1.f;
772+
773+
if (hasZdc) {
774+
const auto zdc = bcForUPC.zdc();
775+
zdcEnergyZNA = zdc.energyCommonZNA();
776+
zdcEnergyZNC = zdc.energyCommonZNC();
737777
qaRegistry.fill(HIST("Data/fitInfo/ampFT0A_vs_ampFT0C"), fitInfo.ampFT0A, fitInfo.ampFT0C);
738-
qaRegistry.fill(HIST("Data/zdc/energyZNA_vs_energyZNC"), zdc.energyCommonZNA(), zdc.energyCommonZNC());
739-
gap = determineGapType(fitInfo.ampFT0A, fitInfo.ampFT0C, zdc.energyCommonZNA(), zdc.energyCommonZNC());
778+
qaRegistry.fill(HIST("Data/zdc/energyZNA_vs_energyZNC"), zdcEnergyZNA, zdcEnergyZNC);
740779
qaRegistry.fill(HIST("Data/hUpcGapAfterSelection"), static_cast<int>(gap));
741780
}
742-
if (gap == GapType::GapA || gap == GapType::GapC) {
743-
fillHistosData<FillMl>(collision, candidates);
744-
}
745-
}
746-
}
781+
for (const auto& candidate : groupedLcCandidates) {
782+
if (!(candidate.hfflag() & 1 << aod::hf_cand_3prong::DecayType::LcToPKPi)) {
783+
continue;
784+
}
785+
if (yCandRecoMax >= 0. && std::abs(HfHelper::yLc(candidate)) > yCandRecoMax) {
786+
continue;
787+
}
788+
const auto pt = candidate.pt();
789+
const auto ptProng0 = candidate.ptProng0();
790+
const auto ptProng1 = candidate.ptProng1();
791+
const auto ptProng2 = candidate.ptProng2();
792+
const auto decayLength = candidate.decayLength();
793+
const auto decayLengthXY = candidate.decayLengthXY();
794+
const auto chi2PCA = candidate.chi2PCA();
795+
const auto cpa = candidate.cpa();
796+
const auto cpaXY = candidate.cpaXY();
747797

748-
GapType determineGapType(float FT0A, float FT0C, float ZNA, float ZNC)
749-
{
750-
constexpr float FT0AThreshold = 100.0;
751-
constexpr float FT0CThreshold = 50.0;
752-
constexpr float ZDCThreshold = 1.0;
753-
if (FT0A < FT0AThreshold && FT0C > FT0CThreshold && ZNA < ZDCThreshold && ZNC > ZDCThreshold) {
754-
return GapType::GapA;
755-
}
756-
if (FT0A > FT0AThreshold && FT0C < FT0CThreshold && ZNA > ZDCThreshold && ZNC < ZDCThreshold) {
757-
return GapType::GapC;
798+
if (candidate.isSelLcToPKPi() >= selectionFlagLc) {
799+
registry.fill(HIST("Data/hMass"), HfHelper::invMassLcToPKPi(candidate));
800+
registry.fill(HIST("Data/hMassVsPtVsNPvContributors"), HfHelper::invMassLcToPKPi(candidate), pt, numPvContributors);
801+
registry.fill(HIST("Data/hMassVsPt"), HfHelper::invMassLcToPKPi(candidate), pt);
802+
}
803+
if (candidate.isSelLcToPiKP() >= selectionFlagLc) {
804+
registry.fill(HIST("Data/hMass"), HfHelper::invMassLcToPiKP(candidate));
805+
registry.fill(HIST("Data/hMassVsPtVsNPvContributors"), HfHelper::invMassLcToPiKP(candidate), pt, numPvContributors);
806+
registry.fill(HIST("Data/hMassVsPt"), HfHelper::invMassLcToPiKP(candidate), pt);
807+
}
808+
registry.fill(HIST("Data/hPt"), pt);
809+
registry.fill(HIST("Data/hPtProng0"), ptProng0);
810+
registry.fill(HIST("Data/hPtProng1"), ptProng1);
811+
registry.fill(HIST("Data/hPtProng2"), ptProng2);
812+
registry.fill(HIST("Data/hd0Prong0"), candidate.impactParameter0());
813+
registry.fill(HIST("Data/hd0Prong1"), candidate.impactParameter1());
814+
registry.fill(HIST("Data/hd0Prong2"), candidate.impactParameter2());
815+
registry.fill(HIST("Data/hd0VsPtProng0"), candidate.impactParameter0(), pt);
816+
registry.fill(HIST("Data/hd0VsPtProng1"), candidate.impactParameter1(), pt);
817+
registry.fill(HIST("Data/hd0VsPtProng2"), candidate.impactParameter2(), pt);
818+
registry.fill(HIST("Data/hDecLength"), decayLength);
819+
registry.fill(HIST("Data/hDecLengthVsPt"), decayLength, pt);
820+
registry.fill(HIST("Data/hDecLengthxy"), decayLengthXY);
821+
registry.fill(HIST("Data/hDecLengthxyVsPt"), decayLengthXY, pt);
822+
registry.fill(HIST("Data/hCt"), HfHelper::ctLc(candidate));
823+
registry.fill(HIST("Data/hCtVsPt"), HfHelper::ctLc(candidate), pt);
824+
registry.fill(HIST("Data/hCPA"), cpa);
825+
registry.fill(HIST("Data/hCPAVsPt"), cpa, pt);
826+
registry.fill(HIST("Data/hCPAxy"), cpaXY);
827+
registry.fill(HIST("Data/hCPAxyVsPt"), cpaXY, pt);
828+
registry.fill(HIST("Data/hDca2"), chi2PCA);
829+
registry.fill(HIST("Data/hDca2VsPt"), chi2PCA, pt);
830+
registry.fill(HIST("Data/hEta"), candidate.eta());
831+
registry.fill(HIST("Data/hEtaVsPt"), candidate.eta(), pt);
832+
registry.fill(HIST("Data/hPhi"), candidate.phi());
833+
registry.fill(HIST("Data/hPhiVsPt"), candidate.phi(), pt);
834+
registry.fill(HIST("hSelectionStatus"), candidate.isSelLcToPKPi(), pt);
835+
registry.fill(HIST("hSelectionStatus"), candidate.isSelLcToPiKP(), pt);
836+
registry.fill(HIST("Data/hImpParErrProng0VsPt"), candidate.errorImpactParameter0(), pt);
837+
registry.fill(HIST("Data/hImpParErrProng1VsPt"), candidate.errorImpactParameter1(), pt);
838+
registry.fill(HIST("Data/hImpParErrProng2VsPt"), candidate.errorImpactParameter2(), pt);
839+
registry.fill(HIST("Data/hDecLenErrVsPt"), candidate.errorDecayLength(), pt);
840+
841+
if (fillTHn) {
842+
float const cent = evaluateCentralityColl(collision);
843+
float occ{-1.};
844+
if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) {
845+
occ = o2::hf_occupancy::getOccupancyColl(collision, occEstimator);
846+
}
847+
double outputBkg(-1), outputPrompt(-1), outputFD(-1);
848+
const float properLifetime = HfHelper::ctLc(candidate) * CtToProperLifetimePs;
849+
850+
auto fillTHnData = [&](bool isPKPi) {
851+
const auto massLc = isPKPi ? HfHelper::invMassLcToPKPi(candidate) : HfHelper::invMassLcToPiKP(candidate);
852+
853+
if constexpr (FillMl) {
854+
const auto& mlProb = isPKPi ? candidate.mlProbLcToPKPi() : candidate.mlProbLcToPiKP();
855+
if (mlProb.size() == NumberOfMlClasses) {
856+
outputBkg = mlProb[MlClassBackground]; /// bkg score
857+
outputPrompt = mlProb[MlClassPrompt]; /// prompt score
858+
outputFD = mlProb[MlClassNonPrompt]; /// non-prompt score
859+
}
860+
/// Fill the ML outputScores and variables of candidate
861+
std::vector<double> valuesToFill{massLc, pt, cent, outputBkg, outputPrompt, outputFD, static_cast<double>(numPvContributors)};
862+
if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) {
863+
valuesToFill.push_back(occ);
864+
}
865+
if (storeProperLifetime) {
866+
valuesToFill.push_back(properLifetime);
867+
}
868+
869+
valuesToFill.push_back(static_cast<double>(gap));
870+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFT0A));
871+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFT0C));
872+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFV0A));
873+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFDDA));
874+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFDDC));
875+
valuesToFill.push_back(static_cast<double>(zdcEnergyZNA));
876+
valuesToFill.push_back(static_cast<double>(zdcEnergyZNC));
877+
878+
registry.get<THnSparse>(HIST("hnLcVarsWithBdt"))->Fill(valuesToFill.data());
879+
} else {
880+
std::vector<double> valuesToFill{massLc, pt, cent, ptProng0, ptProng1, ptProng2, chi2PCA, decayLength, cpa, static_cast<double>(numPvContributors)};
881+
if (storeOccupancy && occEstimator != o2::hf_occupancy::OccupancyEstimator::None) {
882+
valuesToFill.push_back(occ);
883+
}
884+
if (storeProperLifetime) {
885+
valuesToFill.push_back(properLifetime);
886+
}
887+
valuesToFill.push_back(static_cast<double>(gap));
888+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFT0A));
889+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFT0C));
890+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFV0A));
891+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFDDA));
892+
valuesToFill.push_back(static_cast<double>(fitInfo.ampFDDC));
893+
valuesToFill.push_back(static_cast<double>(zdcEnergyZNA));
894+
valuesToFill.push_back(static_cast<double>(zdcEnergyZNC));
895+
registry.get<THnSparse>(HIST("hnLcVars"))->Fill(valuesToFill.data());
896+
}
897+
};
898+
899+
if (candidate.isSelLcToPKPi() >= selectionFlagLc) {
900+
fillTHnData(true);
901+
}
902+
if (candidate.isSelLcToPiKP() >= selectionFlagLc) {
903+
fillTHnData(false);
904+
}
905+
}
906+
}
758907
}
759-
return GapType::DoubleGap;
760908
}
761909

762910
/// Run the analysis on MC data

0 commit comments

Comments
 (0)