Skip to content

Commit dd3985c

Browse files
zhangbiao-phyzhangbiao-phy
authored
[PWGHF] correct daughter-track removal from TPC Q-vector in SP method (#14071)
1 parent 1c2381a commit dd3985c

File tree

1 file changed

+123
-96
lines changed

1 file changed

+123
-96
lines changed

‎PWGHF/D2H/Tasks/taskFlowCharmHadrons.cxx‎

Lines changed: 123 additions & 96 deletions
Original file line numberDiff line numberDiff line change
@@ -64,10 +64,10 @@ namespace o2::aod
6464
{
6565
namespace full
6666
{
67-
DECLARE_SOA_COLUMN(M, m, float); //! Invariant mass of candidate (GeV/c2)
68-
DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c)
69-
DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! ML score of the first configured index
70-
DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! ML score of the second configured index
67+
DECLARE_SOA_COLUMN(M, m, float); //! Invariant mass of candidate (GeV/c2)
68+
DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c)
69+
DECLARE_SOA_COLUMN(MlScore0, mlScore0, float); //! ML score of the first configured index
70+
DECLARE_SOA_COLUMN(MlScore1, mlScore1, float); //! ML score of the second configured index
7171
DECLARE_SOA_COLUMN(ScalarProd, scalarProd, float); //! Scalar product
7272
DECLARE_SOA_COLUMN(Cent, cent, float); //! Centrality
7373
} // namespace full
@@ -125,6 +125,7 @@ struct HfTaskFlowCharmHadrons {
125125
Configurable<float> ptDownSampleMax{"ptDownSampleMax", 10., "Maximum pt for the application of the downsampling factor"};
126126
Configurable<bool> storeResoOccu{"storeResoOccu", false, "Flag to store Occupancy in resolution ThnSparse"};
127127
Configurable<bool> storeEpCosSin{"storeEpCosSin", false, "Flag to store cos and sin of EP angle in ThnSparse"};
128+
Configurable<bool> storeCandEta{"storeCandEta", false, "Flag to store candidates eta"};
128129
Configurable<int> occEstimator{"occEstimator", 0, "Occupancy estimation (0: None, 1: ITS, 2: FT0C)"};
129130
Configurable<bool> saveEpResoHisto{"saveEpResoHisto", false, "Flag to save event plane resolution histogram"};
130131
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
@@ -194,6 +195,7 @@ struct HfTaskFlowCharmHadrons {
194195
ConfigurableAxis thnConfigAxisResoFT0cFV0a{"thnConfigAxisResoFT0cFV0a", {160, -8, 8}, ""};
195196
ConfigurableAxis thnConfigAxisResoFT0cTPCtot{"thnConfigAxisResoFT0cTPCtot", {160, -8, 8}, ""};
196197
ConfigurableAxis thnConfigAxisResoFV0aTPCtot{"thnConfigAxisResoFV0aTPCtot", {160, -8, 8}, ""};
198+
ConfigurableAxis thnConfigAxisCandidateEta{"thnConfigAxisCandidateEta", {100, -5, 5}, ""};
197199

198200
HistogramRegistry registry{"registry", {}};
199201

@@ -214,6 +216,7 @@ struct HfTaskFlowCharmHadrons {
214216
const AxisSpec thnAxisMlTwo{thnConfigAxisMlTwo, "FD score"};
215217
const AxisSpec thnAxisOccupancyITS{thnConfigAxisOccupancyITS, "OccupancyITS"};
216218
const AxisSpec thnAxisOccupancyFT0C{thnConfigAxisOccupancyFT0C, "OccupancyFT0C"};
219+
const AxisSpec thnAxisCandEta{thnConfigAxisCandidateEta, "#eta"};
217220
const AxisSpec thnAxisNoSameBunchPileup{thnConfigAxisNoSameBunchPileup, "NoSameBunchPileup"};
218221
const AxisSpec thnAxisOccupancy{thnConfigAxisOccupancy, "Occupancy"};
219222
const AxisSpec thnAxisNoCollInTimeRangeNarrow{thnConfigAxisNoCollInTimeRangeNarrow, "NoCollInTimeRangeNarrow"};
@@ -231,6 +234,9 @@ struct HfTaskFlowCharmHadrons {
231234
if (storeMl) {
232235
axes.insert(axes.end(), {thnAxisMlOne, thnAxisMlTwo});
233236
}
237+
if (storeCandEta) {
238+
axes.insert(axes.end(), {thnAxisCandEta});
239+
}
234240
if (occEstimator != 0) {
235241
if (occEstimator == 1) {
236242
axes.insert(axes.end(), {thnAxisOccupancyITS, thnAxisNoSameBunchPileup, thnAxisOccupancy,
@@ -317,30 +323,36 @@ struct HfTaskFlowCharmHadrons {
317323
void getQvecDtracks(T1 const& cand,
318324
std::vector<float>& tracksQx,
319325
std::vector<float>& tracksQy,
320-
const float ampl)
326+
const float amplQVec,
327+
QvecEstimator qvecDetector)
321328
{
322-
// TODO: add possibility to consider different weights for the tracks, at the moment only pT is considered;
323-
float const pXTrack0 = cand.pxProng0();
324-
float const pYTrack0 = cand.pyProng0();
325-
float const pTTrack0 = cand.ptProng0();
326-
float const phiTrack0 = std::atan2(pYTrack0, pXTrack0);
327-
float const pXTrack1 = cand.pxProng1();
328-
float const pYTrack1 = cand.pyProng1();
329-
float const pTTrack1 = cand.ptProng1();
330-
float const phiTrack1 = std::atan2(pYTrack1, pXTrack1);
331-
332-
tracksQx.push_back(std::cos(harmonic * phiTrack0) * pTTrack0 / ampl);
333-
tracksQy.push_back(std::sin(harmonic * phiTrack0) * pTTrack0 / ampl);
334-
tracksQx.push_back(std::cos(harmonic * phiTrack1) * pTTrack1 / ampl);
335-
tracksQy.push_back(std::sin(harmonic * phiTrack1) * pTTrack1 / ampl);
329+
auto addProngIfInSubevent = [&](float px, float py, float pz) {
330+
const std::array<float, 3> pVec{px, py, pz};
331+
const float eta = RecoDecay::eta(pVec);
332+
333+
// only subtract daughters that actually contributed to THIS subevent Q
334+
if (qvecDetector == QvecEstimator::TPCPos && eta <= 0.f) {
335+
return;
336+
}
337+
if (qvecDetector == QvecEstimator::TPCNeg && eta >= 0.f) {
338+
return;
339+
}
340+
// for TPCTot: no early return, all prongs contribute
341+
342+
const float pt = RecoDecay::pt(pVec); // or std::hypot(px, py);
343+
const float phi = std::atan2(py, px);
344+
345+
// xQVec,yQVec are normalized with amplQVec, so the daughters' contribution must use the SAME normalization.
346+
tracksQx.push_back(std::cos(harmonic * phi) * pt / amplQVec);
347+
tracksQy.push_back(std::sin(harmonic * phi) * pt / amplQVec);
348+
};
336349

350+
addProngIfInSubevent(cand.pxProng0(), cand.pyProng0(), cand.pzProng0());
351+
addProngIfInSubevent(cand.pxProng1(), cand.pyProng1(), cand.pzProng1());
352+
353+
// 3-prong channels
337354
if constexpr (Channel != DecayChannel::D0ToPiK && Channel != DecayChannel::D0ToKPi) {
338-
float const pXTrack2 = cand.pxProng2();
339-
float const pYTrack2 = cand.pyProng2();
340-
float const pTTrack2 = cand.ptProng2();
341-
float const phiTrack2 = std::atan2(pYTrack2, pXTrack2);
342-
tracksQx.push_back(std::cos(harmonic * phiTrack2) * pTTrack2 / ampl);
343-
tracksQy.push_back(std::sin(harmonic * phiTrack2) * pTTrack2 / ampl);
355+
addProngIfInSubevent(cand.pxProng2(), cand.pyProng2(), cand.pzProng2());
344356
}
345357
}
346358

@@ -353,34 +365,31 @@ struct HfTaskFlowCharmHadrons {
353365
void getQvecXic0Tracks(const T1& cand,
354366
std::vector<float>& tracksQx,
355367
std::vector<float>& tracksQy,
356-
float ampl)
368+
float amplQVec,
369+
QvecEstimator qvecDetector)
357370
{
358-
// add possibility to consider different weights for the tracks, at the moment only pT is considered;
359-
float const pXTrack0 = cand.pxPosV0Dau();
360-
float const pYTrack0 = cand.pyPosV0Dau();
361-
float const pTTrack0 = std::hypot(pXTrack0, pYTrack0);
362-
float const phiTrack0 = std::atan2(pXTrack0, pYTrack0);
363-
float const pXTrack1 = cand.pxNegV0Dau();
364-
float const pYTrack1 = cand.pyNegV0Dau();
365-
float const pTTrack1 = std::hypot(pXTrack1, pYTrack1);
366-
float const phiTrack1 = std::atan2(pXTrack1, pYTrack1);
367-
float const pYTrack2 = cand.pxBachFromCasc();
368-
float const pXTrack2 = cand.pyBachFromCasc();
369-
float const pTTrack2 = std::hypot(pXTrack2, pYTrack2);
370-
float const phiTrack2 = std::atan2(pXTrack2, pYTrack2);
371-
float const pXTrack3 = cand.pxBachFromCharmBaryon();
372-
float const pYTrack3 = cand.pyBachFromCharmBaryon();
373-
float const pTTrack3 = std::hypot(pXTrack3, pYTrack3);
374-
float const phiTrack3 = std::atan2(pXTrack3, pYTrack3);
375-
376-
tracksQx.push_back(std::cos(harmonic * phiTrack0) * pTTrack0 / ampl);
377-
tracksQy.push_back(std::sin(harmonic * phiTrack0) * pTTrack0 / ampl);
378-
tracksQx.push_back(std::cos(harmonic * phiTrack1) * pTTrack1 / ampl);
379-
tracksQy.push_back(std::sin(harmonic * phiTrack1) * pTTrack1 / ampl);
380-
tracksQx.push_back(std::cos(harmonic * phiTrack2) * pTTrack2 / ampl);
381-
tracksQy.push_back(std::sin(harmonic * phiTrack2) * pTTrack2 / ampl);
382-
tracksQx.push_back(std::cos(harmonic * phiTrack3) * pTTrack3 / ampl);
383-
tracksQy.push_back(std::sin(harmonic * phiTrack3) * pTTrack3 / ampl);
371+
auto addProngIfInSubevent = [&](float px, float py, float pz) {
372+
const std::array<float, 3> pVec{px, py, pz};
373+
const float eta = RecoDecay::eta(pVec);
374+
375+
if (qvecDetector == QvecEstimator::TPCPos && eta <= 0.f) {
376+
return;
377+
}
378+
if (qvecDetector == QvecEstimator::TPCNeg && eta >= 0.f) {
379+
return;
380+
}
381+
382+
const float pt = RecoDecay::pt(pVec);
383+
const float phi = std::atan2(py, px);
384+
385+
tracksQx.push_back(std::cos(harmonic * phi) * pt / amplQVec);
386+
tracksQy.push_back(std::sin(harmonic * phi) * pt / amplQVec);
387+
};
388+
389+
addProngIfInSubevent(cand.pxPosV0Dau(), cand.pyPosV0Dau(), cand.pzPosV0Dau());
390+
addProngIfInSubevent(cand.pxNegV0Dau(), cand.pyNegV0Dau(), cand.pzNegV0Dau());
391+
addProngIfInSubevent(cand.pxBachFromCasc(), cand.pyBachFromCasc(), cand.pzBachFromCasc());
392+
addProngIfInSubevent(cand.pxBachFromCharmBaryon(), cand.pyBachFromCharmBaryon(), cand.pzBachFromCharmBaryon());
384393
}
385394
/// Compute the delta psi in the range [0, pi/harmonic]
386395
/// \param psi1 is the first angle
@@ -408,6 +417,7 @@ struct HfTaskFlowCharmHadrons {
408417
/// Fill THnSparse
409418
/// \param mass is the invariant mass of the candidate
410419
/// \param pt is the transverse momentum of the candidate
420+
/// \param eta is the pseudorapidityof the candidate
411421
/// \param cent is the centrality of the collision
412422
/// \param cosNPhi is the cosine of the n*phi angle
413423
/// \param sinNPhi is the sine of the n*phi angle
@@ -418,6 +428,7 @@ struct HfTaskFlowCharmHadrons {
418428
/// \param hfevselflag flag of the collision associated to utilsEvSelHf.h
419429
void fillThn(const float mass,
420430
const float pt,
431+
const float eta,
421432
const float cent,
422433
const float cosNPhi,
423434
const float sinNPhi,
@@ -427,40 +438,50 @@ struct HfTaskFlowCharmHadrons {
427438
const float occupancy,
428439
const o2::hf_evsel::HfCollisionRejectionMask hfevselflag)
429440
{
441+
auto hSparse = registry.get<THnSparse>(HIST("hSparseFlowCharm"));
442+
const int ndim = hSparse->GetNdimensions();
443+
444+
std::vector<double> values;
445+
values.reserve(ndim);
446+
447+
values.push_back(mass);
448+
values.push_back(pt);
449+
values.push_back(cent);
450+
values.push_back(sp);
451+
452+
if (storeEP) {
453+
values.push_back(cosNPhi);
454+
values.push_back(sinNPhi);
455+
values.push_back(cosDeltaPhi);
456+
}
457+
458+
if (storeMl) {
459+
values.push_back(outputMl[0]);
460+
values.push_back(outputMl[1]);
461+
}
462+
463+
if (storeCandEta) {
464+
values.push_back(eta);
465+
}
466+
430467
if (occEstimator != 0) {
431-
std::vector<int> evtSelFlags = getEventSelectionFlags(hfevselflag);
432-
if (storeMl) {
433-
if (storeEP) {
434-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, cosNPhi, sinNPhi, cosDeltaPhi, outputMl[0], outputMl[1], occupancy,
435-
evtSelFlags[0], evtSelFlags[1], evtSelFlags[2], evtSelFlags[3], evtSelFlags[4]);
436-
} else {
437-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, outputMl[0], outputMl[1], occupancy,
438-
evtSelFlags[0], evtSelFlags[1], evtSelFlags[2], evtSelFlags[3], evtSelFlags[4]);
439-
}
440-
} else {
441-
if (storeEP) {
442-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, cosNPhi, sinNPhi, cosDeltaPhi, occupancy,
443-
evtSelFlags[0], evtSelFlags[1], evtSelFlags[2], evtSelFlags[3], evtSelFlags[4]);
444-
} else {
445-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, occupancy,
446-
evtSelFlags[0], evtSelFlags[1], evtSelFlags[2], evtSelFlags[3], evtSelFlags[4]);
447-
}
448-
}
449-
} else {
450-
if (storeMl) {
451-
if (storeEP) {
452-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, cosNPhi, sinNPhi, cosDeltaPhi, outputMl[0], outputMl[1]);
453-
} else {
454-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, outputMl[0], outputMl[1]);
455-
}
456-
} else {
457-
if (storeEP) {
458-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp, cosNPhi, sinNPhi, cosDeltaPhi);
459-
} else {
460-
registry.fill(HIST("hSparseFlowCharm"), mass, pt, cent, sp);
461-
}
462-
}
468+
auto evtSelFlags = getEventSelectionFlags(hfevselflag);
469+
values.push_back(occupancy);
470+
values.push_back(evtSelFlags[0]);
471+
values.push_back(evtSelFlags[1]);
472+
values.push_back(evtSelFlags[2]);
473+
values.push_back(evtSelFlags[3]);
474+
values.push_back(evtSelFlags[4]);
475+
}
476+
477+
if (static_cast<int>(values.size()) != ndim) {
478+
LOGF(fatal,
479+
"hSparseFlowCharm: number of filled dimensions (%d) "
480+
"does not match THnSparse dimensionality (%d).",
481+
static_cast<int>(values.size()), ndim);
463482
}
483+
484+
hSparse->Fill(values.data());
464485
}
465486

466487
/// Check if the collision is selected
@@ -555,7 +576,6 @@ struct HfTaskFlowCharmHadrons {
555576
float yQVec = qVecs[1];
556577
float const amplQVec = qVecs[2];
557578
float const evtPl = epHelper.GetEventPlane(xQVec, yQVec, harmonic);
558-
int nProngs = 3;
559579

560580
for (const auto& candidate : candidates) {
561581
float massCand = 0.;
@@ -590,7 +610,6 @@ struct HfTaskFlowCharmHadrons {
590610
}
591611
}
592612
} else if constexpr (std::is_same_v<T1, CandD0Data> || std::is_same_v<T1, CandD0DataWMl>) {
593-
nProngs = 2;
594613
switch (Channel) {
595614
case DecayChannel::D0ToPiK:
596615
massCand = HfHelper::invMassD0ToPiK(candidate);
@@ -664,27 +683,35 @@ struct HfTaskFlowCharmHadrons {
664683

665684
float ptCand = 0.;
666685
float phiCand = 0.;
686+
float etaCand = 0.;
667687

668688
if constexpr (std::is_same_v<T1, CandXic0Data> || std::is_same_v<T1, CandXic0DataWMl>) {
669689
ptCand = RecoDecay::sqrtSumOfSquares(candidate.pxCharmBaryon(), candidate.pyCharmBaryon());
670690
phiCand = std::atan2(candidate.pxCharmBaryon(), candidate.pyCharmBaryon());
691+
etaCand = candidate.etaBachFromCharmBaryon();
671692
} else {
672693
ptCand = candidate.pt();
673694
phiCand = candidate.phi();
695+
etaCand = candidate.eta();
674696
}
675697

676-
// If TPC is used for the SP estimation, the tracks of the hadron candidate must be removed from the TPC Q vector to avoid double counting
677-
if (qvecDetector == QvecEstimator::TPCNeg || qvecDetector == QvecEstimator::TPCPos) {
678-
float const ampl = amplQVec - static_cast<float>(nProngs);
679-
std::vector<float> tracksQx = {};
680-
std::vector<float> tracksQy = {};
698+
// If TPC is used for the SP estimation, the tracks of the hadron candidate must be removed from the corresponding TPC Q vector to avoid self-correlations
699+
if (qvecDetector == QvecEstimator::TPCNeg ||
700+
qvecDetector == QvecEstimator::TPCPos ||
701+
qvecDetector == QvecEstimator::TPCTot) {
702+
703+
std::vector<float> tracksQx;
704+
std::vector<float> tracksQy;
705+
706+
// IMPORTANT: use the ORIGINAL amplitude used to build this Q-vector
681707
if constexpr (std::is_same_v<T1, CandXic0Data> || std::is_same_v<T1, CandXic0DataWMl>) {
682-
// std::cout<<candidate.pxProng0()<<std::endl;
683-
getQvecXic0Tracks(candidate, tracksQx, tracksQy, ampl);
708+
getQvecXic0Tracks(candidate, tracksQx, tracksQy, amplQVec, static_cast<QvecEstimator>(qvecDetector.value));
684709
} else {
685-
getQvecDtracks<Channel>(candidate, tracksQx, tracksQy, ampl);
710+
getQvecDtracks<Channel>(candidate, tracksQx, tracksQy, amplQVec, static_cast<QvecEstimator>(qvecDetector.value));
686711
}
687-
for (auto iTrack{0u}; iTrack < tracksQx.size(); ++iTrack) {
712+
713+
// subtract daughters' contribution from the (normalized) Q-vector
714+
for (std::size_t iTrack = 0; iTrack < tracksQx.size(); ++iTrack) {
688715
xQVec -= tracksQx[iTrack];
689716
yQVec -= tracksQy[iTrack];
690717
}
@@ -710,7 +737,7 @@ struct HfTaskFlowCharmHadrons {
710737
}
711738
}
712739
if (fillSparse) {
713-
fillThn(massCand, ptCand, cent, cosNPhi, sinNPhi, cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag);
740+
fillThn(massCand, ptCand, etaCand, cent, cosNPhi, sinNPhi, cosDeltaPhi, scalprodCand, outputMl, occupancy, hfevflag);
714741
}
715742
}
716743
}

0 commit comments

Comments
 (0)