Skip to content

Add QA, datalink and Mc labels for 3body decay #3235

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mpuccio merged 7 commits into from
Aug 14, 2023
Merged

Add QA, datalink and Mc labels for 3body decay #3235

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
a2433d8
Add QA, Decay3Body datalink and Mc labels
wang-yuanzhe Aug 1, 2023
2f08658
Merge branch 'AliceO2Group:master' into yuanzhe_dev
wang-yuanzhe Aug 1, 2023
b871272
Please consider the following formatting changes
alibuild Aug 1, 2023
fcaea81
Merge pull request #4 from alibuild/alibot-cleanup-3235
wang-yuanzhe Aug 1, 2023
e9caad1
Add MC Dalitz plot
wang-yuanzhe Aug 1, 2023
0a46d00
Change getMassPDG to physics constant
wang-yuanzhe Aug 1, 2023
d7138c4
Change getMassPDG to physics constant
wang-yuanzhe Aug 1, 2023
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
29 changes: 25 additions & 4 deletions PWGLF/DataModel/Vtx3BodyTables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ DECLARE_SOA_INDEX_COLUMN_FULL(Track0, track0, int, Tracks, "_0"); //!
DECLARE_SOA_INDEX_COLUMN_FULL(Track1, track1, int, Tracks, "_1"); //!
DECLARE_SOA_INDEX_COLUMN_FULL(Track2, track2, int, Tracks, "_2"); //!
DECLARE_SOA_INDEX_COLUMN(Collision, collision); //!
DECLARE_SOA_INDEX_COLUMN(Decay3Body, decay3body); //!

// General 3 body Vtx properties: position, momentum
DECLARE_SOA_COLUMN(PxTrack0, pxtrack0, float); //! track0 px at min
Expand Down Expand Up @@ -109,7 +110,7 @@ DECLARE_SOA_EXPRESSION_COLUMN(Pz, pz, //! 3 body vtx pz
} // namespace vtx3body

DECLARE_SOA_TABLE_FULL(StoredVtx3BodyDatas, "Vtx3BodyDatas", "AOD", "Vtx3BodyDATA", //!
o2::soa::Index<>, vtx3body::Track0Id, vtx3body::Track1Id, vtx3body::Track2Id, vtx3body::CollisionId,
o2::soa::Index<>, vtx3body::Track0Id, vtx3body::Track1Id, vtx3body::Track2Id, vtx3body::CollisionId, vtx3body::Decay3BodyId,
vtx3body::X, vtx3body::Y, vtx3body::Z,
vtx3body::PxTrack0, vtx3body::PyTrack0, vtx3body::PzTrack0,
vtx3body::PxTrack1, vtx3body::PyTrack1, vtx3body::PzTrack1,
Expand Down Expand Up @@ -152,11 +153,31 @@ namespace vtx3body
DECLARE_SOA_INDEX_COLUMN(Vtx3BodyData, vtx3BodyData); //! Index to Vtx3BodyData entry
}

DECLARE_SOA_TABLE(Vtx3BodyDataLink, "AOD", "VtxDATALINK", //! Joinable table with Decay3bodys which links to Vtx3BodyData which is not produced for all entries
DECLARE_SOA_TABLE(Decay3BodyDataLink, "AOD", "DECAY3BODYLINK", //! Joinable table with Decay3bodys which links to Vtx3BodyData which is not produced for all entries
vtx3body::Vtx3BodyDataId);

using Vtxs3BodyLinked = soa::Join<Decay3Bodys, Vtx3BodyDataLink>;
using Vtx3BodyLinked = Vtxs3BodyLinked::iterator;
using Decay3BodysLinked = soa::Join<Decay3Bodys, Decay3BodyDataLink>;
using Decay3BodyLinked = Decay3BodysLinked::iterator;

// Definition of labels for Vtx3BodyDatas
namespace mcvtx3bodylabel
{
DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); //! MC particle for Vtx3BodyDatas
} // namespace mcvtx3bodylabel

DECLARE_SOA_TABLE(McVtx3BodyLabels, "AOD", "MCVTXLABEL", //! Table joinable with Vtx3BodyData containing the MC labels
mcvtx3bodylabel::McParticleId);
using McVtx3BodyLabel = McVtx3BodyLabels::iterator;

// Definition of labels for Decay3Bodys // Full table, joinable with Decay3Bodys (CAUTION: NOT WITH Vtx3BodyDATA)
namespace mcfullvtx3bodylabel
{
DECLARE_SOA_INDEX_COLUMN(McParticle, mcParticle); //! MC particle for Decay3Bodys
} // namespace mcfullvtx3bodylabel

DECLARE_SOA_TABLE(McFullVtx3BodyLabels, "AOD", "MCFULLVTXLABEL", //! Table joinable with Decay3Bodys
mcfullvtx3bodylabel::McParticleId);
using McFullVtx3BodyLabel = McFullVtx3BodyLabels::iterator;

} // namespace o2::aod
#endif // PWGLF_DATAMODEL_VTX3BODYTABLES_H_
75 changes: 54 additions & 21 deletions PWGLF/TableProducer/hypertriton3bodybuilder.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,7 +49,7 @@ using MyTracksIU = FullTracksExtIU;
// in case requested
using LabeledTracks = soa::Join<FullTracksExtIU, aod::McTrackLabels>;

struct hypertriton3bodybuilder {
struct hypertriton3bodyBuilder {

Produces<aod::StoredVtx3BodyDatas> vtx3bodydata;
Service<o2::ccdb::BasicCCDBManager> ccdb;
Expand Down Expand Up @@ -188,7 +188,7 @@ struct hypertriton3bodybuilder {
initCCDB(bc);
registry.fill(HIST("hEventCounter"), 0.5);

for (auto& vtx3body : decay3bodys) { // FIXME: turn into combination(...)
for (auto& vtx3body : decay3bodys) {

registry.fill(HIST("hVtx3BodyCounter"), 0.5);

Expand Down Expand Up @@ -264,7 +264,7 @@ struct hypertriton3bodybuilder {
registry.fill(HIST("hVtx3BodyCounter"), 6.5);

vtx3bodydata(
t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(),
t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(), vtx3body.globalIndex(),
pos[0], pos[1], pos[2],
p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1], p2[2],
fitter3body.getChi2AtPCACandidate(),
Expand All @@ -273,9 +273,30 @@ struct hypertriton3bodybuilder {
}
};

struct hypertriton3bodyDataLinkBuilder {
Produces<aod::Decay3BodyDataLink> vtxdataLink;

void init(InitContext const&) {}

void process(aod::Decay3Bodys const& decay3bodytable, aod::Vtx3BodyDatas const& vtxdatatable)
{
std::vector<int> lIndices;
lIndices.reserve(decay3bodytable.size());
for (int ii = 0; ii < decay3bodytable.size(); ii++)
lIndices[ii] = -1;
for (auto& vtxdata : vtxdatatable) {
lIndices[vtxdata.decay3bodyId()] = vtxdata.globalIndex();
}
for (int ii = 0; ii < decay3bodytable.size(); ii++) {
vtxdataLink(lIndices[ii]);
}
}
};

struct hypertriton3bodyLabelBuilder {

Produces<aod::McV0Labels> vtxlabels;
Produces<aod::McVtx3BodyLabels> vtxlabels;
Produces<aod::McFullVtx3BodyLabels> vtxfulllabels;

// for bookkeeping purposes: how many V0s come from same mother etc
HistogramRegistry registry{
Expand Down Expand Up @@ -314,32 +335,38 @@ struct hypertriton3bodyLabelBuilder {
}
PROCESS_SWITCH(hypertriton3bodyLabelBuilder, processDoNotBuildLabels, "Do not produce MC label tables", true);

void processBuildLabels(aod::Collisions::iterator const& collision, aod::Decay3Bodys const& decay3bodys, LabeledTracks const&, aod::McParticles const& particlesMC)
void processBuildLabels(aod::Decay3BodysLinked const& decay3bodys, aod::Vtx3BodyDatas const& vtx3bodydatas, LabeledTracks const&, aod::McParticles const& particlesMC)
{
for (auto& vtx : decay3bodys) {
std::vector<int> lIndices;
lIndices.reserve(vtx3bodydatas.size());
for (int ii = 0; ii < vtx3bodydatas.size(); ii++) {
lIndices[ii] = -1;
}

for (auto& decay3body : decay3bodys) {

int lLabel = -1;
int lPDG = -1;
float lPt = -1;
double MClifetime = -1;
bool is3bodyDecay = false; // all considered V0s

auto lTrack0 = vtx.track0_as<LabeledTracks>();
auto lTrack1 = vtx.track1_as<LabeledTracks>();
auto lTrack2 = vtx.track2_as<LabeledTracks>();
auto lTrack0 = decay3body.track0_as<LabeledTracks>();
auto lTrack1 = decay3body.track1_as<LabeledTracks>();
auto lTrack2 = decay3body.track2_as<LabeledTracks>();
registry.fill(HIST("hLabelCounter"), 0.5);

// Association check
// There might be smarter ways of doing this in the future
if (!lTrack0.has_mcParticle() || !lTrack1.has_mcParticle() || !lTrack2.has_mcParticle()) {
vtxlabels(-1);
vtxfulllabels(-1);
continue;
}
auto lMCTrack0 = lTrack0.mcParticle_as<aod::McParticles>();
auto lMCTrack1 = lTrack1.mcParticle_as<aod::McParticles>();
auto lMCTrack2 = lTrack2.mcParticle_as<aod::McParticles>();
if (!lMCTrack0.has_mothers() || !lMCTrack1.has_mothers() || !lMCTrack2.has_mothers()) {
vtxlabels(-1);
vtxfulllabels(-1);
continue;
}

Expand All @@ -350,22 +377,22 @@ struct hypertriton3bodyLabelBuilder {
lLabel = lMother1.globalIndex();
lPt = lMother1.pt();
lPDG = lMother1.pdgCode();
MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - lMother2.vx(), lMCTrack2.vy() - lMother2.vy(), lMCTrack2.vz() - lMother2.vz()) * 2.991 / lMother2.p();
MClifetime = RecoDecay::sqrtSumOfSquares(lMCTrack2.vx() - lMother2.vx(), lMCTrack2.vy() - lMother2.vy(), lMCTrack2.vz() - lMother2.vz()) * o2::constants::physics::MassHyperTriton / lMother2.p();
is3bodyDecay = true; // vtxs with the same mother
}
}
}
} // end association check
if (!is3bodyDecay) {
vtxlabels(-1);
vtxfulllabels(-1);
continue;
}
registry.fill(HIST("hLabelCounter"), 1.5);

// Intended for cross-checks only
// N.B. no rapidity cut!
if (lPDG == 1010010030 && lMCTrack0.pdgCode() == 2212 && lMCTrack1.pdgCode() == -211 && lMCTrack2.pdgCode() == 1000010020) {
double hypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kPiPlus), 1.87561});
double hypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged, o2::constants::physics::MassDeuteron});
registry.fill(HIST("hHypertritonCounter"), 0.5);
registry.fill(HIST("hHypertriton"), lPt);
registry.fill(HIST("hHypertritonMass"), hypertritonMCMass);
Expand All @@ -384,7 +411,7 @@ struct hypertriton3bodyLabelBuilder {
}
}
if (lPDG == -1010010030 && lMCTrack0.pdgCode() == 211 && lMCTrack1.pdgCode() == -2212 && lMCTrack2.pdgCode() == -1000010020) {
double antiHypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kProton), 1.87561});
double antiHypertritonMCMass = RecoDecay::m(array{array{lMCTrack0.px(), lMCTrack0.py(), lMCTrack0.pz()}, array{lMCTrack1.px(), lMCTrack1.py(), lMCTrack1.pz()}, array{lMCTrack2.px(), lMCTrack2.py(), lMCTrack2.pz()}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron});
registry.fill(HIST("hHypertritonCounter"), 2.5);
registry.fill(HIST("hAntiHypertriton"), lPt);
registry.fill(HIST("hAntiHypertritonMass"), antiHypertritonMCMass);
Expand All @@ -403,24 +430,30 @@ struct hypertriton3bodyLabelBuilder {
}
}

// Construct label table (note: this will be joinable with V0Datas)
vtxlabels(
lLabel);
// Construct label table
vtxfulllabels(lLabel);
if (decay3body.vtx3BodyDataId() != -1) {
lIndices[decay3body.vtx3BodyDataId()] = lLabel;
}
}
for (int ii = 0; ii < vtx3bodydatas.size(); ii++) {
vtxlabels(lIndices[ii]);
}
}
PROCESS_SWITCH(hypertriton3bodyLabelBuilder, processBuildLabels, "Produce MC label tables", false);
};

struct hypertriton3bodyinitializer {
struct hypertriton3bodyInitializer {
Spawns<aod::Vtx3BodyDatas> vtx3bodydatas;
void init(InitContext const&) {}
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
return WorkflowSpec{
adaptAnalysisTask<hypertriton3bodyinitializer>(cfgc),
adaptAnalysisTask<hypertriton3bodyBuilder>(cfgc),
adaptAnalysisTask<hypertriton3bodyDataLinkBuilder>(cfgc),
adaptAnalysisTask<hypertriton3bodyLabelBuilder>(cfgc),
adaptAnalysisTask<hypertriton3bodybuilder>(cfgc),
adaptAnalysisTask<hypertriton3bodyInitializer>(cfgc),
};
}
8 changes: 4 additions & 4 deletions PWGLF/TableProducer/hypertriton3bodyfinder.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -496,10 +496,10 @@ struct hypertriton3bodyFinder {
float p2V0 = pt2V0 + pV0[2] * pV0[2], ptV0 = std::sqrt(pt2V0);
// apply mass selections
// float p2Pos = pP[0] * pP[0] + pP[1] * pP[1] + pP[2] * pP[2], p2Neg = pN[0] * pN[0] + pN[1] * pN[1] + pN[2] * pN[2];
float massV0LambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{RecoDecay::getMassPDG(kProton), RecoDecay::getMassPDG(kPiPlus)});
float massV0AntiLambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{RecoDecay::getMassPDG(kPiPlus), RecoDecay::getMassPDG(kProton)});
float massV0LambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{o2::constants::physics::MassProton, o2::constants::physics::MassPionCharged});
float massV0AntiLambdaHyp = RecoDecay::m(array{array{pP[0], pP[1], pP[2]}, array{pN[0], pN[1], pN[2]}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton});
float massMargin = 20 * (0.001 * (1. + 0.5 * ptV0)) + 0.07;
if (massV0LambdaHyp - RecoDecay::getMassPDG(kLambda0) > massMargin && massV0AntiLambdaHyp - RecoDecay::getMassPDG(kLambda0) > massMargin) {
if (massV0LambdaHyp - o2::constants::physics::MassLambda > massMargin && massV0AntiLambdaHyp - o2::constants::physics::MassLambda > massMargin) {
continue;
}
registry.fill(HIST("hV0Counter"), 7.5);
Expand Down Expand Up @@ -642,7 +642,7 @@ struct hypertriton3bodyFinder {

// Fix: Reconstruction H3L DCA Check
vtx3bodydata(
t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(),
t0.globalIndex(), t1.globalIndex(), t2.globalIndex(), collision.globalIndex(), 0,
vertexXYZ[0], vertexXYZ[1], vertexXYZ[2],
p0[0], p0[1], p0[2], p1[0], p1[1], p1[2], p2[0], p2[1], p2[2],
fitter3body.getChi2AtPCACandidate(),
Expand Down
Loading