Skip to content

Commit 5800a84

Browse files
Update taskCd.cxx
1 parent 2f4650a commit 5800a84

File tree

1 file changed

+215
-13
lines changed

1 file changed

+215
-13
lines changed

PWGHF/D2H/Tasks/taskCd.cxx

Lines changed: 215 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -22,6 +22,7 @@
2222

2323
#include "Common/DataModel/Centrality.h"
2424
#include "Common/DataModel/EventSelection.h"
25+
#include "Common/DataModel/PIDResponseITS.h"
2526

2627
#include <Framework/ASoA.h>
2728
#include <Framework/AnalysisDataModel.h>
@@ -52,18 +53,87 @@ using namespace o2::hf_centrality;
5253
using namespace o2::hf_occupancy;
5354
using namespace o2::hf_evsel;
5455

56+
namespace o2::aod
57+
{
58+
namespace full
59+
{
60+
// Candidate kinematics
61+
DECLARE_SOA_COLUMN(M, m, float); //! Invariant mass of candidate (GeV/c^2)
62+
DECLARE_SOA_COLUMN(Pt, pt, float); //! Transverse momentum of candidate (GeV/c)
63+
DECLARE_SOA_COLUMN(PtProng0, ptProng0, float); //! Transverse momentum of prong 0 (GeV/c)
64+
DECLARE_SOA_COLUMN(PtProng1, ptProng1, float); //! Transverse momentum of prong 1 (GeV/c)
65+
DECLARE_SOA_COLUMN(PtProng2, ptProng2, float); //! Transverse momentum of prong 2 (GeV/c)
66+
DECLARE_SOA_COLUMN(ImpactParameter0, impactParameter0, float); //! Impact parameter (DCA to PV) of prong 0 (cm)
67+
DECLARE_SOA_COLUMN(ImpactParameter1, impactParameter1, float); //! Impact parameter (DCA to PV) of prong 1 (cm)
68+
DECLARE_SOA_COLUMN(ImpactParameter2, impactParameter2, float); //! Impact parameter (DCA to PV) of prong 2 (cm)
69+
DECLARE_SOA_COLUMN(DecayLength, decayLength, float); //! Decay length (3D) of candidate (cm)
70+
DECLARE_SOA_COLUMN(DecayLengthXY, decayLengthXY, float); //! Decay length in transverse plane (cm)
71+
DECLARE_SOA_COLUMN(Cpa, cpa, float); //! Cosine of pointing angle (3D)
72+
DECLARE_SOA_COLUMN(CpaXY, cpaXY, float); //! Cosine of pointing angle in XY plane
73+
DECLARE_SOA_COLUMN(NSigmaTpcDe, nSigmaTpcDe, float); //! TPC nσ for deuteron hypothesis
74+
DECLARE_SOA_COLUMN(NSigmaTpcKa, nSigmaTpcKa, float); //! TPC nσ for kaon hypothesis
75+
DECLARE_SOA_COLUMN(NSigmaTpcPi, nSigmaTpcPi, float); //! TPC nσ for pion hypothesis
76+
DECLARE_SOA_COLUMN(NSigmaItsDe, nSigmaItsDe, float); //! ITS nσ for deuteron hypothesis
77+
DECLARE_SOA_COLUMN(NSigmaTofDe, nSigmaTofDe, float); //! TOF nσ for deuteron hypothesis
78+
DECLARE_SOA_COLUMN(NSigmaTofKa, nSigmaTofKa, float); //! TOF nσ for kaon hypothesis
79+
DECLARE_SOA_COLUMN(NSigmaTofPi, nSigmaTofPi, float); //! TOF nσ for pion hypothesis
80+
DECLARE_SOA_COLUMN(NItsClusters, nItsClusters, int8_t); //! Number of ITS clusters used in the track fit
81+
DECLARE_SOA_COLUMN(NItsNClusterSize, nItsNClusterSize, int8_t); //! Number of ITS clusters size used in the track fit
82+
DECLARE_SOA_COLUMN(NTpcClusters, nTpcClusters, int8_t); //! Number of TPC clusters used in the track fit
83+
DECLARE_SOA_COLUMN(NTpcSignalsDe, nTpcSignalsDe, int8_t); //! Number of TPC signas
84+
DECLARE_SOA_COLUMN(NItsSignalsDe, nItsSignalsDe, int8_t); //! Number of ITS signas
85+
DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int8_t); //! Candidates falg
86+
DECLARE_SOA_COLUMN(Cent, cent, float); //! Centrality
87+
DECLARE_SOA_COLUMN(GIndexCol, gIndexCol, int); //! Global index for the collisionAdd commentMore actions
88+
DECLARE_SOA_COLUMN(TimeStamp, timeStamp, int64_t); //! Timestamp for the collision
89+
} // namespace full
90+
91+
// Full table: include ALL columns declared above
92+
DECLARE_SOA_TABLE(HfCandCd, "AOD", "HFCANDCD",
93+
full::M,
94+
full::Pt,
95+
full::PtProng0,
96+
full::PtProng1,
97+
full::PtProng2,
98+
full::ImpactParameter0,
99+
full::ImpactParameter1,
100+
full::ImpactParameter2,
101+
full::DecayLength,
102+
full::Cpa,
103+
full::NSigmaTpcDe,
104+
full::NSigmaTpcKa,
105+
full::NSigmaTpcPi,
106+
full::NSigmaItsDe,
107+
full::NSigmaTofDe,
108+
full::NSigmaTofKa,
109+
full::NSigmaTofPi,
110+
full::NItsClusters,
111+
full::NItsNClusterSize,
112+
full::NTpcClusters,
113+
full::NTpcSignalsDe,
114+
full::NItsSignalsDe,
115+
full::CandidateSelFlag,
116+
full::Cent,
117+
full::GIndexCol,
118+
full::TimeStamp);
119+
} // namespace o2::aod
120+
55121
struct HfTaskCd {
122+
123+
Produces<o2::aod::HfCandCd> rowCandCd;
56124
Configurable<int> selectionFlagCd{"selectionFlagCd", 1, "Selection Flag for Cd"};
57125
Configurable<std::vector<double>> binsPt{"binsPt", std::vector<double>{hf_cuts_cd_to_de_k_pi::vecBinsPt}, "pT bin limits"};
58126
Configurable<bool> fillTHn{"fillTHn", false, "fill THn"};
127+
Configurable<bool> fillTree{"fillTree", false, "Flag to fill candiates tree"};
59128

60129
SliceCache cache;
61130

62131
using CollisionsWEvSel = soa::Join<aod::Collisions, aod::EvSels>;
63132
using CollisionsWithEvSelFT0C = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs>;
64133
using CollisionsWithEvSelFT0M = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms>;
65134

66-
using CdCandidates = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelCd>>;
135+
using CdCandidates = soa::Filtered<soa::Join<aod::HfCand3Prong, aod::HfSelCd, aod::HfCand3ProngWPidPiKaDe>>;
136+
using HFTracks = soa::Join<aod::FullTracks, aod::TracksDCA, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullDe, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullDe>;
67137

68138
Filter filterSelectCandidates = aod::hf_sel_candidate_cd::isSelCdToDeKPi >= selectionFlagCd;
69139
Preslice<aod::HfCand3Prong> candCdPerCollision = aod::hf_cand::collisionId;
@@ -139,6 +209,11 @@ struct HfTaskCd {
139209
registry.add("Data/hImpParErrProng0", "3-prong candidates;prong 0 impact parameter error (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
140210
registry.add("Data/hImpParErrProng1", "3-prong candidates;prong 1 impact parameter error (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
141211
registry.add("Data/hImpParErrProng2", "3-prong candidates;prong 2 impact parameter error (cm);entries", {HistType::kTH2F, {{100, -1., 1.}, {binsPt, "#it{p}_{T} (GeV/#it{c})"}}});
212+
registry.add("Data/hNsigmaTPCDeVsP", "deuteron;n#sigma^{TPC}_{d};#it{p} (GeV/#it{c})", {HistType::kTH2F, {{200, -10.f, 10.f}, {200, -6.f, 6.f}}});
213+
registry.add("Data/hNsigmaTOFDeVsP", "deuteron;n#sigma^{TOF}_{d};#it{p} (GeV/#it{c})", {HistType::kTH2F, {{200, -10.f, 10.f}, {200, -6.f, 6.f}}});
214+
registry.add("Data/hNsigmaITSDeVsP", "deuteron;n#sigma^{ITS}_{d};#it{p} (GeV/#it{c})", {HistType::kTH2F, {{200, -10.f, 10.f}, {200, -6.f, 6.f}}});
215+
registry.add("Data/hTPCSignalDeVsP", "deuteron;n#sigma^{ITS}_{d};#it{p} (GeV/#it{c})", {HistType::kTH2F, {{200, -10.f, 10.f}, {2000, 0, 2000}}});
216+
registry.add("Data/hITSSignalDeVsP", "deuteron;n#sigma^{ITS}_{d};#it{p} (GeV/#it{c})", {HistType::kTH2F, {{200, -10.f, 10.f}, {20, 0, 20}}});
142217

143218
if (fillTHn) {
144219
const AxisSpec thnAxisMass{thnConfigAxisMass, "inv. mass (de K #pi) (GeV/#it{c}^{2})"};
@@ -156,13 +231,33 @@ struct HfTaskCd {
156231
}
157232
}
158233

234+
// taken from: https://github.com/AliceO2Group/O2Physics/blob/master/EventFiltering/PWGCF/CFFilterAll.cxx
235+
template <typename T>
236+
float itsSignal(T const& track)
237+
{
238+
uint32_t clsizeflag = track.itsClusterSizes();
239+
auto clSizeLayer0 = (clsizeflag >> (0 * 4)) & 0xf;
240+
auto clSizeLayer1 = (clsizeflag >> (1 * 4)) & 0xf;
241+
auto clSizeLayer2 = (clsizeflag >> (2 * 4)) & 0xf;
242+
auto clSizeLayer3 = (clsizeflag >> (3 * 4)) & 0xf;
243+
auto clSizeLayer4 = (clsizeflag >> (4 * 4)) & 0xf;
244+
auto clSizeLayer5 = (clsizeflag >> (5 * 4)) & 0xf;
245+
auto clSizeLayer6 = (clsizeflag >> (6 * 4)) & 0xf;
246+
int numLayers = 7;
247+
int sumClusterSizes = clSizeLayer1 + clSizeLayer2 + clSizeLayer3 + clSizeLayer4 + clSizeLayer5 + clSizeLayer6 + clSizeLayer0;
248+
float cosLamnda = 1. / std::cosh(track.eta());
249+
return (static_cast<float>(sumClusterSizes) / numLayers) * cosLamnda;
250+
};
251+
159252
/// Fill histograms for real data
160-
template <typename CollType, typename CandType>
161-
void fillHistosData(CollType const& collision, CandType const& candidates)
253+
template <typename CollType, typename CandType, typename TrackType>
254+
void fillHistosData(CollType const& collision, CandType const& candidates, TrackType const&)
162255
{
163256
auto thisCollId = collision.globalIndex();
164257
auto groupedCdCandidates = candidates.sliceBy(candCdPerCollision, thisCollId);
165258
auto numPvContributors = collision.numContrib();
259+
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
260+
int64_t timeStamp = bc.timestamp();
166261

167262
for (const auto& candidate : groupedCdCandidates) {
168263
if (!TESTBIT(candidate.hfflag(), aod::hf_cand_3prong::DecayType::CdToDeKPi)) {
@@ -178,6 +273,13 @@ struct HfTaskCd {
178273
const auto chi2PCA = candidate.chi2PCA();
179274
const auto cpa = candidate.cpa();
180275
const auto cpaXY = candidate.cpaXY();
276+
float invMassCd = 0.f;
277+
if (candidate.isSelCdToDeKPi() >= selectionFlagCd) {
278+
invMassCd = HfHelper::invMassCdToDeKPi(candidate);
279+
}
280+
if (candidate.isSelCdToPiKDe() >= selectionFlagCd) {
281+
invMassCd = HfHelper::invMassCdToPiKDe(candidate);
282+
}
181283

182284
if (candidate.isSelCdToDeKPi() >= selectionFlagCd) {
183285
registry.fill(HIST("Data/hMass"), HfHelper::invMassCdToDeKPi(candidate));
@@ -219,8 +321,9 @@ struct HfTaskCd {
219321
registry.fill(HIST("Data/hImpParErrProng1"), candidate.errorImpactParameter1(), pt);
220322
registry.fill(HIST("Data/hImpParErrProng2"), candidate.errorImpactParameter2(), pt);
221323

324+
float const cent = o2::hf_centrality::getCentralityColl(collision);
325+
222326
if (fillTHn) {
223-
float const cent = o2::hf_centrality::getCentralityColl(collision);
224327
double massCd(-1);
225328
if (candidate.isSelCdToDeKPi() >= selectionFlagCd) {
226329
massCd = HfHelper::invMassCdToDeKPi(candidate);
@@ -233,40 +336,139 @@ struct HfTaskCd {
233336
registry.get<THnSparse>(HIST("hnCdVars"))->Fill(valuesToFill.data());
234337
}
235338
}
339+
340+
if (fillTree) {
341+
int candFlag = -999;
342+
343+
float nSigmaTpcDe = 0.f, nSigmaTpcKa = 0.f, nSigmaTpcPi = 0.f;
344+
float nSigmaItsDe = 0.f;
345+
float nSigmaTofDe = 0.f, nSigmaTofKa = 0.f, nSigmaTofPi = 0.f;
346+
347+
int itsNClusterDe = 0;
348+
int itsNClusterSizeDe = 0;
349+
int tpcNClusterDe = 0;
350+
351+
float tpcSignalsDe = 0.f;
352+
float itsSignalsDe = 0.f;
353+
354+
nSigmaTpcKa = candidate.nSigTpcKa1();
355+
nSigmaTofKa = candidate.nSigTofKa1();
356+
357+
const bool selDeKPi = (candidate.isSelCdToDeKPi() >= 1);
358+
const bool selPiKDe = (candidate.isSelCdToPiKDe() >= 1);
359+
360+
TrackType trackDe{};
361+
bool hasDe = false;
362+
363+
if (selDeKPi) {
364+
auto trk0 = candidate.template prong0_as<TrackType>();
365+
trackDe = trk0;
366+
hasDe = true;
367+
368+
candFlag = 1;
369+
nSigmaTpcDe = candidate.nSigTpcDe0();
370+
nSigmaTofDe = candidate.nSigTofDe0();
371+
nSigmaTpcPi = candidate.nSigTpcPi2();
372+
nSigmaTofPi = candidate.nSigTofPi2();
373+
} else if (selPiKDe) {
374+
auto trk2 = candidate.template prong2_as<TrackType>();
375+
trackDe = trk2;
376+
hasDe = true;
377+
378+
candFlag = -1;
379+
nSigmaTpcDe = candidate.nSigTpcDe2();
380+
nSigmaTofDe = candidate.nSigTofDe2();
381+
nSigmaTpcPi = candidate.nSigTpcPi0();
382+
nSigmaTofPi = candidate.nSigTofPi0();
383+
}
384+
385+
if (hasDe) {
386+
const float pSigned = trackDe.p() * trackDe.sign();
387+
388+
nSigmaItsDe = trackDe.itsNSigmaDe();
389+
itsNClusterDe = trackDe.itsNCls();
390+
itsNClusterSizeDe = trackDe.itsClusterSizes();
391+
tpcNClusterDe = trackDe.tpcNClsCrossedRows();
392+
tpcSignalsDe = trackDe.tpcSignal();
393+
itsSignalsDe = itsSignal(trackDe);
394+
395+
// De PID QA
396+
registry.fill(HIST("Data/hNsigmaTPCDeVsP"), pSigned, nSigmaTpcDe);
397+
registry.fill(HIST("Data/hNsigmaTOFDeVsP"), pSigned, nSigmaTofDe);
398+
registry.fill(HIST("Data/hNsigmaITSDeVsP"), pSigned, nSigmaItsDe);
399+
registry.fill(HIST("Data/hTPCSignalDeVsP"), pSigned, tpcSignalsDe);
400+
registry.fill(HIST("Data/hITSSignalDeVsP"), pSigned, itsSignalsDe);
401+
}
402+
403+
rowCandCd(
404+
invMassCd,
405+
pt,
406+
ptProng0,
407+
ptProng1,
408+
ptProng2,
409+
candidate.impactParameter0(),
410+
candidate.impactParameter1(),
411+
candidate.impactParameter2(),
412+
decayLength,
413+
cpa,
414+
nSigmaTpcDe,
415+
nSigmaTpcKa,
416+
nSigmaTpcPi,
417+
nSigmaItsDe,
418+
nSigmaTofDe,
419+
nSigmaTofKa,
420+
nSigmaTofPi,
421+
itsNClusterDe,
422+
itsNClusterSizeDe,
423+
tpcNClusterDe,
424+
tpcSignalsDe,
425+
itsSignalsDe,
426+
candFlag,
427+
cent,
428+
collision.globalIndex(),
429+
timeStamp);
430+
}
236431
}
237432
}
238433
/// Run the analysis on real data
239-
template <typename CollType, typename CandType>
434+
template <typename CollType, typename CandType, typename TrackType>
240435
void runAnalysisPerCollisionData(CollType const& collisions,
241-
CandType const& candidates)
436+
CandType const& candidates,
437+
TrackType const& tracks)
242438
{
243439

244440
for (const auto& collision : collisions) {
245-
fillHistosData(collision, candidates);
441+
fillHistosData(collision, candidates, tracks);
246442
}
247443
}
248444

249445
void processDataStd(CollisionsWEvSel const& collisions,
250446
CdCandidates const& selectedCdCandidates,
251-
aod::Tracks const&)
447+
HFTracks const& tracks)
252448
{
253-
runAnalysisPerCollisionData(collisions, selectedCdCandidates);
449+
// inlcude ITS PID information
450+
auto tracksWithItsPid = soa::Attach<HFTracks, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe>(tracks);
451+
runAnalysisPerCollisionData(collisions, selectedCdCandidates, tracksWithItsPid);
254452
}
255453
PROCESS_SWITCH(HfTaskCd, processDataStd, "Process Data with the standard method", true);
256454

257455
void processDataStdWithFT0C(CollisionsWithEvSelFT0C const& collisions,
258456
CdCandidates const& selectedCdCandidates,
259-
aod::Tracks const&)
457+
HFTracks const& tracks)
260458
{
261-
runAnalysisPerCollisionData(collisions, selectedCdCandidates);
459+
// inlcude ITS PID information
460+
auto tracksWithItsPid = soa::Attach<HFTracks, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe>(tracks);
461+
runAnalysisPerCollisionData(collisions, selectedCdCandidates, tracksWithItsPid);
262462
}
263463
PROCESS_SWITCH(HfTaskCd, processDataStdWithFT0C, "Process real data with the standard method and with FT0C centrality", false);
264464

265465
void processDataStdWithFT0M(CollisionsWithEvSelFT0M const& collisions,
266466
CdCandidates const& selectedCdCandidates,
267-
aod::Tracks const&)
467+
HFTracks const& tracks)
268468
{
269-
runAnalysisPerCollisionData(collisions, selectedCdCandidates);
469+
// inlcude ITS PID information
470+
auto tracksWithItsPid = soa::Attach<HFTracks, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe>(tracks);
471+
runAnalysisPerCollisionData(collisions, selectedCdCandidates, tracksWithItsPid);
270472
}
271473
PROCESS_SWITCH(HfTaskCd, processDataStdWithFT0M, "Process real data with the standard method and with FT0M centrality", false);
272474
};

0 commit comments

Comments
 (0)