Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
81 changes: 58 additions & 23 deletions PWGHF/D2H/TableProducer/dataCreatorHiddenCharmReduced.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -76,9 +76,10 @@ struct HfDataCreatorHiddenCharmReduced {
Produces<aod::HcSelTracks> hfTrackLite;

struct : ConfigurableGroup {
// track quality
Configurable<bool> fillHistograms{"fillHistograms", true, "Fill proton QA histograms"};
Configurable<bool> fillSparses{"fillSparses", false, "Flag to enable sparse filling"};
Configurable<bool> selectProtons{"selectProtons", true, "Select protons"};
// track quality
Configurable<int> itsNClsMin{"itsNClsMin", 5, "Minimum number of ITS clusters"};
Configurable<int> tpcNClsFoundMin{"tpcNClsFoundMin", 50, "Minimum number of found TPC clusters"};
Configurable<int> tpcNClsCrossedRowsMin{"tpcNClsCrossedRowsMin", 80, "Minimum number of crossed TPC rows"};
Expand All @@ -90,6 +91,7 @@ struct HfDataCreatorHiddenCharmReduced {
Configurable<std::vector<double>> binsPtTrack{"binsPtTrack", std::vector<double>{hf_cuts_single_track::vecBinsPtTrack}, "Track pT bin limits for DCA cuts"};
Configurable<LabeledArray<double>> cutsTrack{"cutsTrack", {hf_cuts_single_track::CutsTrack[0], hf_cuts_single_track::NBinsPtTrack, hf_cuts_single_track::NCutVarsTrack, hf_cuts_single_track::labelsPtTrack, hf_cuts_single_track::labelsCutVarTrack}, "Single-track DCA selections per pT bin"};
// DCA
Configurable<bool> applyPtDependentDcaTrackSelection{"applyPtDependentDcaTrackSelection", false, "Flag to enable pT-dependent DCAxy and DCAz selection on single tracks"};
Configurable<std::vector<float>> paramsDCAxyPtDep{"paramsDCAxyPtDep", std::vector<float>{0.0010, 0.0080, 0.73}, "Parameters for pT-dependent DCAxy cut: p0, p1, p2 for cut = p0 + p1/pt^p2"};
Configurable<std::vector<float>> paramsDCAzPtDep{"paramsDCAzPtDep", std::vector<float>{-0.0044, 0.0152, 0.47}, "Parameters for pT-dependent DCAz cut: p0, p1, p2 for cut = p0 + p1/pt^p2"};
// PID
Expand Down Expand Up @@ -147,17 +149,30 @@ struct HfDataCreatorHiddenCharmReduced {
const AxisSpec axisPt{360, 0., 36., "#it{p}_{T}^{proton} (GeV/#it{c})"};
const AxisSpec axisEta{100, -1., 1., "#eta"};
const AxisSpec axisDca{400, -2., 2., "DCA_{xy} to primary vertex (cm)"};
const AxisSpec axisNSigma{100, -5., 5., "n#sigma"};
const AxisSpec axisNSigmaTPC{100, -5., 5., "n#sigma_{TPC}"};
const AxisSpec axisNSigmaTOF{100, -5., 5., "n#sigma_{TOF}"};
const AxisSpec axisMass{100, 2.85, 3.25, "p#overline{p} invariant mass (GeV/#it{c}^{2})"};

registry.add("hPzVtx", "Z position of primary vertex for selected tracks;z_{vtx} (cm);entries", {HistType::kTH1D, {AxisSpec{200, -20., 20., "z_{vtx} (cm)"}}});
registry.add("hPtNoCuts", "All associated tracks;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1D, {axisPt}});
registry.add("hPtCutsProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});entries", {HistType::kTH1D, {axisPt}});
registry.add("hEtaCutsProton", "Selected proton tracks;#eta;entries", {HistType::kTH1D, {axisEta}});
registry.add("hDCAToPrimXYVsPtCutsProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});DCA_{xy} to primary vertex (cm)", {HistType::kTH2D, {axisPt, axisDca}});
registry.add("hNSigmaTPCProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});n#sigma_{TPC}", {HistType::kTH2D, {axisPt, axisNSigma}});
registry.add("hNSigmaTOFProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});n#sigma_{TOF}", {HistType::kTH2D, {axisPt, axisNSigma}});
registry.add("hInvMass", "Invariant mass of selected proton with all other tracks in the event;#it{p}_{T}^{proton} (GeV/#it{c});invariant mass with other tracks (GeV/#it{c}^{2})", {HistType::kTH2D, {axisPt, AxisSpec{100, 2.85, 3.25, "invariant mass with other tracks (GeV/#it{c}^{2})"}}});
registry.add("hNSigmaTPCProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});n#sigma_{TPC}^{p}", {HistType::kTH2D, {axisPt, axisNSigmaTPC}});
registry.add("hNSigmaTOFProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});n#sigma_{TOF}^{p}", {HistType::kTH2D, {axisPt, axisNSigmaTOF}});
registry.add("hNSigmaTPCKaon", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); n#sigma_{TPC}^{K}", {HistType::kTH2D, {axisPt, axisNSigmaTPC}});
registry.add("hNSigmaTOFKaon", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); n#sigma_{TOF}^{K}", {HistType::kTH2D, {axisPt, axisNSigmaTOF}});
registry.add("hNSigmaTPCPion", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); n#sigma_{TPC}^{#pi}", {HistType::kTH2D, {axisPt, axisNSigmaTPC}});
registry.add("hNSigmaTOFPion", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); n#sigma_{TOF}^{#pi}", {HistType::kTH2D, {axisPt, axisNSigmaTOF}});
registry.add("hTpcNClsFound", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); TpcNClsFound", {HistType::kTH2D, {axisPt, {200, 0., 200, "TpcNClsFound"}}});
registry.add("hTpcNClsCrossedRows", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); TpcNClsCrossedRows", {HistType::kTH2D, {axisPt, {200, 0., 200, "TpcNClsCrossedRows"}}});
registry.add("hTpcChi2NCl", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); TpcChi2NCl", {HistType::kTH2D, {axisPt, {10, 0., 10, "TpcChi2NCl"}}});
registry.add("hItsNCls", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});ItsNCls", {HistType::kTH2D, {axisPt, {7, 0, 7, "ItsNCls"}}});
registry.add("hItsChi2NCl", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c}); ItsChi2NCl", {HistType::kTH2D, {axisPt, {10, 0., 10, "ItsChi2NCl"}}});
registry.add("hInvMass", "Invariant mass of selected proton with all other tracks in the event;#it{p}_{T}^{proton} (GeV/#it{c});", {HistType::kTH2D, {axisPt, axisMass}});
registry.add("hInvMassSparse", "Invariant mass of selected proton with all other tracks in the event;#it{p}_{T}^{proton} (GeV/#it{c});", {HistType::kTH2D, {axisPt, axisMass}});
registry.add("hDeDxTPCProton", "Selected proton tracks;#it{p}_{T}^{track} (GeV/#it{c});TPC dE/dx (a.u.)", {HistType::kTH2D, {axisPt, AxisSpec{100, 0., 200., "TPC dE/dx (a.u.)"}}});
registry.add("hMassPtCutVars", "Charmonia candidates;#it{M} (p#overline{p}) (GeV/#it{c}^{2});#it{p}_{T}(cc) (GeV/#it{c}); n#sigma_{TPC}^{p1}; n#sigma_{TOF}^{p1}; dca_{xy}^{p1}; n#sigma_{TPC}^{p2}; n#sigma_{TOF}^{p2}; dca_{xy}^{p2}", {HistType::kTHnSparseF, {axisMass, axisPt, axisNSigmaTPC, axisNSigmaTOF, axisDca, axisNSigmaTPC, axisNSigmaTOF, axisDca}});
}

// init HF event selection helper
Expand Down Expand Up @@ -266,8 +281,10 @@ struct HfDataCreatorHiddenCharmReduced {
if (!isSelectedTrackDca(config.binsPtTrack, config.cutsTrack, track.pt(), track.dcaXY(), track.dcaZ())) {
return false;
}
if (dcaSigma(track.pt(), config.paramsDCAxyPtDep.value[0], config.paramsDCAxyPtDep.value[1], config.paramsDCAxyPtDep.value[2]) > std::abs(dcaXY) || dcaSigma(track.pt(), config.paramsDCAzPtDep.value[0], config.paramsDCAzPtDep.value[1], config.paramsDCAzPtDep.value[2]) > std::abs(dcaZ)) {
return false;
if (config.applyPtDependentDcaTrackSelection.value) {
if (std::abs(dcaXY) > dcaSigma(track.pt(), config.paramsDCAxyPtDep.value[0], config.paramsDCAxyPtDep.value[1], config.paramsDCAxyPtDep.value[2]) || std::abs(dcaZ) > dcaSigma(track.pt(), config.paramsDCAzPtDep.value[0], config.paramsDCAzPtDep.value[1], config.paramsDCAzPtDep.value[2])) {
return false;
}
}
return isSelectedPid(track);
}
Expand Down Expand Up @@ -304,37 +321,55 @@ struct HfDataCreatorHiddenCharmReduced {
if (!isSelectedTrack(trk)) {
continue;
}
std::array pVecProton{trk.pVector()};
hfTrackLite(trk.globalIndex(), collision.globalIndex(), pVecProton[0], pVecProton[1], pVecProton[2], trk.sign(), static_cast<uint8_t>(TrackType::Proton));
selectedTrackIds.push_back(trk.globalIndex());
selectedTrackIds.push_back(trkId.trackId());

if (config.fillHistograms) {
registry.fill(HIST("hPtCutsProton"), trk.pt());
registry.fill(HIST("hEtaCutsProton"), trk.eta());
registry.fill(HIST("hDCAToPrimXYVsPtCutsProton"), trk.pt(), trk.dcaXY());
registry.fill(HIST("hNSigmaTPCProton"), trk.pt(), trk.tpcNSigmaPr());
registry.fill(HIST("hDeDxTPCProton"), trk.pt(), trk.tpcSignal());
registry.fill(HIST("hNSigmaTPCPion"), trk.pt(), trk.tpcNSigmaPi());
registry.fill(HIST("hNSigmaTPCKaon"), trk.pt(), trk.tpcNSigmaKa());
if (trk.hasTOF()) {
registry.fill(HIST("hNSigmaTOFProton"), trk.pt(), trk.tofNSigmaPr());
registry.fill(HIST("hNSigmaTOFPion"), trk.pt(), trk.tofNSigmaPi());
registry.fill(HIST("hNSigmaTOFKaon"), trk.pt(), trk.tofNSigmaKa());
}
registry.fill(HIST("hTpcNClsFound"), trk.pt(), trk.tpcNClsFound());
registry.fill(HIST("hTpcNClsCrossedRows"), trk.pt(), trk.tpcNClsCrossedRows());
registry.fill(HIST("hTpcChi2NCl"), trk.pt(), trk.tpcChi2NCl());
registry.fill(HIST("hItsNCls"), trk.pt(), trk.itsNCls());
registry.fill(HIST("hItsChi2NCl"), trk.pt(), trk.itsChi2NCl());
}
}
// skip event if not at least a pair of tracks skimming selection
if (selectedTrackIds.size() < NDaughtersCharmMeson) {
continue;
}
hfReducedCollision(collision.posX(), collision.posY(), collision.posZ(), collision.numContrib(), hfRejMap, bz);
for (size_t i = 0; i < selectedTrackIds.size(); ++i) {
auto t1 = tracks.rawIteratorAt(selectedTrackIds[i]);
std::array pVec1{t1.pVector()};
for (size_t j = i + 1; j < selectedTrackIds.size(); ++j) {
auto t2 = tracks.rawIteratorAt(selectedTrackIds[j]);
if (t1.sign() * t2.sign() > 0) {
continue;
}
std::array pVec2{t2.pVector()};
float invMass = RecoDecay::m(std::array{pVec1, pVec2}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassProton});
float ptEtac = RecoDecay::pt(RecoDecay::sumOfVec(pVec1, pVec2));
if (config.fillHistograms) {

auto reducedEventId = hfReducedCollision.lastIndex();
for (const auto& trkId : selectedTrackIds) {
auto track = tracks.rawIteratorAt(trkId);
std::array pVec{track.pVector()};
hfTrackLite(track.globalIndex(), reducedEventId, pVec[0], pVec[1], pVec[2], track.sign(), static_cast<uint8_t>(TrackType::Proton));
}

if (config.fillHistograms.value) {
for (size_t i = 0; i < selectedTrackIds.size(); ++i) {
auto t1 = tracks.rawIteratorAt(selectedTrackIds[i]);
std::array pVec1{t1.pVector()};
for (size_t j = i + 1; j < selectedTrackIds.size(); ++j) {
auto t2 = tracks.rawIteratorAt(selectedTrackIds[j]);
if (t1.sign() * t2.sign() > 0)
continue;
std::array pVec2{t2.pVector()};
float invMass = RecoDecay::m(std::array{pVec1, pVec2}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassProton});
float ptEtac = RecoDecay::pt(RecoDecay::sumOfVec(pVec1, pVec2));
registry.fill(HIST("hInvMass"), ptEtac, invMass);
if (config.fillSparses.value) {
registry.fill(HIST("hMassPtCutVars"), invMass, ptEtac, t1.tpcNSigmaPr(), t1.tofNSigmaPr(), t1.dcaXY(), t2.tpcNSigmaPr(), t2.tofNSigmaPr(), t2.dcaXY());
}
}
}
}
Expand Down
38 changes: 22 additions & 16 deletions PWGHF/D2H/Tasks/taskHiddenCharm.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -52,8 +52,8 @@ struct HfTaskHiddenCharm {
Configurable<int> centEstimator{"centEstimator", 0, "Centrality estimation (None: 0, FT0A: 1, FT0C: 2, FT0M: 3, FV0A: 4, NTracksPV: 5, FT0CVariant2: 6)"};
Configurable<float> centralityMin{"centralityMin", 0.f, "Minimum accepted centrality"};
Configurable<float> centralityMax{"centralityMax", 100.f, "Maximum accepted centrality"};
Configurable<bool> fillOnlyUnlikeSign{"fillOnlyUnlikeSign", true, "Fill only unlike-sign proton pairs"};
Configurable<bool> fillOnlyLikeSign{"fillOnlyLikeSign", true, "Fill only like-sign proton pairs"};
Configurable<bool> fillUnlikeSign{"fillUnlikeSign", true, "Fill unlike-sign proton pairs"};
Configurable<bool> fillLikeSign{"fillLikeSign", true, "Fill like-sign proton pairs"};

SliceCache cache;

Expand Down Expand Up @@ -94,34 +94,40 @@ struct HfTaskHiddenCharm {
void fillEtac(TCollisions const& collision,
TProtonIds const& protonIds)
{
float cent{-1.f};
float cent{0.f};
if constexpr (CentEstimator != o2::hf_centrality::CentralityEstimator::None) {
cent = o2::hf_centrality::getCentralityColl(collision, centEstimator);
if (cent < centralityMin || cent >= centralityMax) {
return; // skip events outside the centrality range
}
}

for (const auto& proton1 : protonIds) {
for (const auto& proton2 : protonIds) {
if (proton1.trackId() >= proton2.trackId()) {
continue; // avoid double counting and self-pairs
}
for (auto it1 = protonIds.begin(); it1 != protonIds.end(); ++it1) {
for (auto it2 = it1 + 1; it2 != protonIds.end(); ++it2) {

const auto& proton1 = *it1;
const auto& proton2 = *it2;
const int sign = (proton1.sign() * proton2.sign() > 0) ? 1 : -1;
if ((sign == 1 && !fillOnlyLikeSign) || (sign == -1 && !fillOnlyUnlikeSign)) {

if (sign == 1 && !fillLikeSign)
continue;
}
if (sign == -1 && !fillUnlikeSign)
continue;

std::array<float, 3> pVec1{proton1.px(), proton1.py(), proton1.pz()};
std::array<float, 3> pVec2{proton2.px(), proton2.py(), proton2.pz()};

float invMass = RecoDecay::m(std::array{pVec1, pVec2}, std::array{o2::constants::physics::MassProton, o2::constants::physics::MassProton});
float ptEtac = RecoDecay::pt(RecoDecay::sumOfVec(pVec1, pVec2));
registry.fill(HIST("hSparseHiddenCharm"), invMass, ptEtac, sign, cent);
float pt = RecoDecay::pt(RecoDecay::sumOfVec(pVec1, pVec2));

registry.fill(HIST("hSparseHiddenCharm"), invMass, pt, sign, cent);
if (sign == 1) {
registry.fill(HIST("hPtVsInvMassLikeSign"), invMass, ptEtac);
} else if (sign == -1) {
registry.fill(HIST("hPtVsInvMassUnlikeSign"), invMass, ptEtac);
registry.fill(HIST("hPtVsInvMassAllSign"), invMass, ptEtac);
registry.fill(HIST("hPtVsInvMassLikeSign"), invMass, pt);
} else {
registry.fill(HIST("hPtVsInvMassUnlikeSign"), invMass, pt);
}

registry.fill(HIST("hPtVsInvMassAllSign"), invMass, pt);
}
}
}
Expand Down
Loading