From 16f84fe1f5e3fa4146f1f3235107d0e1bd5e7258 Mon Sep 17 00:00:00 2001 From: Paul Buehler Date: Mon, 30 Oct 2023 10:33:45 +0100 Subject: [PATCH] PWG-UD: UDTurials03[a,b] and UDTutorials04 (#3722) * veto only on forward tracks with good timeing * add function to detect rho-> mu mu Mc events * protect cases with missing associations * improve processing of Mc truth * minor updates * replace UDTutorials03 with UDTutorials03a and UDTutorials03b * clang-format * clang-format * remove unused variable * clang-format --- PWGUD/Core/DGSelector.h | 16 +- PWGUD/Core/UDHelpers.h | 43 ++ PWGUD/TableProducer/DGBCCandProducer.cxx | 32 +- PWGUD/TableProducer/DGCandProducer.cxx | 319 +++++++++----- PWGUD/Tasks/DGCandAnalyzer.cxx | 9 +- Tutorials/PWGUD/CMakeLists.txt | 15 + Tutorials/PWGUD/UDTutorial_03a.cxx | 380 ++++++++++++++++ Tutorials/PWGUD/UDTutorial_03b.cxx | 525 +++++++++++++++++++++++ Tutorials/PWGUD/UDTutorial_04.cxx | 478 +++++++++++++++++++++ 9 files changed, 1695 insertions(+), 122 deletions(-) create mode 100644 Tutorials/PWGUD/UDTutorial_03a.cxx create mode 100644 Tutorials/PWGUD/UDTutorial_03b.cxx create mode 100644 Tutorials/PWGUD/UDTutorial_04.cxx diff --git a/PWGUD/Core/DGSelector.h b/PWGUD/Core/DGSelector.h index 12a2bd038c4..c354a0c7297 100644 --- a/PWGUD/Core/DGSelector.h +++ b/PWGUD/Core/DGSelector.h @@ -54,9 +54,9 @@ class DGSelector // forward tracks LOGF(debug, "FwdTracks %i", fwdtracks.size()); if (!diffCuts.withFwdTracks()) { - // only consider tracks with MID (good timing) for (auto& fwdtrack : fwdtracks) { LOGF(info, " %i / %f / %f / %f / %f", fwdtrack.trackType(), fwdtrack.eta(), fwdtrack.pt(), fwdtrack.p(), fwdtrack.trackTimeRes()); + // only consider tracks with MID (good timing) if (fwdtrack.trackType() == 0 || fwdtrack.trackType() == 3) { return 2; } @@ -158,12 +158,14 @@ class DGSelector } // no activity in muon arm - LOGF(debug, "FwdTracks %i", fwdtracks.size()); - for (auto& fwdtrack : fwdtracks) { - LOGF(debug, " %i / %f / %f / %f", fwdtrack.trackType(), fwdtrack.eta(), fwdtrack.pt(), fwdtrack.p()); - } - if (fwdtracks.size() > 0) { - return 2; + if (!diffCuts.withFwdTracks()) { + for (auto& fwdtrack : fwdtracks) { + LOGF(info, " %i / %f / %f / %f / %f", fwdtrack.trackType(), fwdtrack.eta(), fwdtrack.pt(), fwdtrack.p(), fwdtrack.trackTimeRes()); + // only consider tracks with MID (good timing) + if (fwdtrack.trackType() == 0 || fwdtrack.trackType() == 3) { + return 2; + } + } } // number of tracks diff --git a/PWGUD/Core/UDHelpers.h b/PWGUD/Core/UDHelpers.h index fa60436f7f7..8917d400bf6 100644 --- a/PWGUD/Core/UDHelpers.h +++ b/PWGUD/Core/UDHelpers.h @@ -596,6 +596,25 @@ bool isPythiaCDE(T MCparts) return false; } +// ----------------------------------------------------------------------------- +// In rho -> mu+ + mu- events generated with STARlight the stack starts with +// 443013, 13, -13 or 443013, -13, 13 +template +bool isSTARLightRhomumu(T MCparts) +{ + if (MCparts.size() < 3) { + return false; + } else { + if (MCparts.iteratorAt(0).pdgCode() != 443013) + return false; + if (abs(MCparts.iteratorAt(1).pdgCode()) != 13) + return false; + if (MCparts.iteratorAt(2).pdgCode() != -MCparts.iteratorAt(1).pdgCode()) + return false; + } + return true; +} + // ----------------------------------------------------------------------------- // In pp events produced with GRANIITTI the stack starts with // 22212/22212/99/22212/2212/99/90 @@ -630,6 +649,30 @@ bool isGraniittiCDE(T MCparts) return true; } +// ----------------------------------------------------------------------------- +// function to select MC events of interest +template +int isOfInterest(T MCparts) +{ + + // PYTHIA CDE + if (isPythiaCDE(MCparts)) { + return 1; + } + + // GRANIITTI CDE + if (isGraniittiCDE(MCparts)) { + return 2; + } + + // STARLIGHT rho -> mu+ + mu- + if (isSTARLightRhomumu(MCparts)) { + return 3; + } + + return 0; +} + // ----------------------------------------------------------------------------- // Invariant mass of GRANIITTI generated event template diff --git a/PWGUD/TableProducer/DGBCCandProducer.cxx b/PWGUD/TableProducer/DGBCCandProducer.cxx index fb89dfad037..03d5c0a8347 100644 --- a/PWGUD/TableProducer/DGBCCandProducer.cxx +++ b/PWGUD/TableProducer/DGBCCandProducer.cxx @@ -62,6 +62,9 @@ struct tracksWGTInBCs { void processBarrel(BCs const& bcs, CCs const& collisions, TCs const& tracks, ATs const& ambTracks) { // run number + if (bcs.size() <= 0) { + return; + } int rnum = bcs.iteratorAt(0).runNumber(); // container to sort tracks with good timing according to their matching/closest BC @@ -73,7 +76,7 @@ struct tracksWGTInBCs { for (auto const& track : tracks) { registry.get(HIST("barrel/Tracks"))->Fill(0., 1.); - // is this track aPV track? + // is this track a PV track? if (track.isPVContributor()) { registry.get(HIST("barrel/Tracks"))->Fill(1., 1.); } @@ -96,17 +99,23 @@ struct tracksWGTInBCs { // compute the BC closest in time auto firstCompatibleBC = ambTracksSlice.begin().bc().begin().globalBC(); - LOGF(debug, "Track time %f", track.trackTime()); + LOGF(debug, "First compatible BC %d Track time %f", firstCompatibleBC, track.trackTime()); closestBC = (uint64_t)(firstCompatibleBC + (track.trackTime() / o2::constants::lhc::LHCBunchSpacingNS)); } else { // this track is not ambiguous, has hence a unique association to a collision/BC + if (!track.has_collision()) { + continue; + } auto collision = track.collision_as(); + if (!collision.has_foundBC()) { + continue; + } closestBC = collision.foundBC_as().globalBC(); } // update tracksInBCList - LOGF(debug, "Closest BC %d", closestBC); + LOGF(debug, "Updating tracksInBCList with %d", closestBC); tracksInBCList[closestBC].emplace_back((int32_t)track.globalIndex()); } } @@ -129,8 +138,8 @@ struct tracksWGTInBCs { break; } } - tracksWGTInBCs(indBCToSave, rnum, tracksInBC.first, tracksInBC.second); LOGF(debug, " BC %i/%u with %i tracks with good timing", indBCToSave, tracksInBC.first, tracksInBC.second.size()); + tracksWGTInBCs(indBCToSave, rnum, tracksInBC.first, tracksInBC.second); } LOGF(debug, "barrel done"); } @@ -142,6 +151,9 @@ struct tracksWGTInBCs { void processForward(BCs& bcs, CCs& collisions, aod::FwdTracks& fwdTracks, aod::AmbiguousFwdTracks& ambFwdTracks) { // run number + if (bcs.size() <= 0) { + return; + } int rnum = bcs.iteratorAt(0).runNumber(); // container to sort forward tracks according to their matching/closest BC @@ -176,11 +188,18 @@ struct tracksWGTInBCs { (fwdTrack.trackTime() / o2::constants::lhc::LHCBunchSpacingNS)); } else { // this track is not ambiguous, has hence a unique association to a collision/BC + if (!fwdTrack.has_collision()) { + continue; + } auto collision = fwdTrack.collision_as(); - closestBC = collision.bc_as().globalBC(); + if (!collision.has_foundBC()) { + continue; + } + closestBC = collision.foundBC_as().globalBC(); } // update tracksInBCList + LOGF(debug, "Updating fwdTracksInBCList with %d", closestBC); fwdTracksInBCList[closestBC].emplace_back((int32_t)fwdTrack.globalIndex()); } } @@ -205,7 +224,7 @@ struct tracksWGTInBCs { fwdTracksWGTInBCs(indBCToSave, rnum, fwdTracksInBC.first, fwdTracksInBC.second); LOGF(debug, " BC %i/%u with %i forward tracks with good timing", indBCToSave, fwdTracksInBC.first, fwdTracksInBC.second.size()); } - LOGF(debug, "fwd done"); + LOGF(debug, "forward done"); } PROCESS_SWITCH(tracksWGTInBCs, processForward, "Process forward tracks", false); @@ -539,6 +558,7 @@ struct DGBCCandProducer { bcnum = tibc.bcnum(); } } + bool withCollision = false; while (bc2go || tibc2go) { LOGF(debug, "Testing bc %d/%d/%d", bcnum, bc.globalBC(), tibc.bcnum()); diff --git a/PWGUD/TableProducer/DGCandProducer.cxx b/PWGUD/TableProducer/DGCandProducer.cxx index c2a0d8eb8f2..3dd6f35f5c3 100644 --- a/PWGUD/TableProducer/DGCandProducer.cxx +++ b/PWGUD/TableProducer/DGCandProducer.cxx @@ -136,8 +136,8 @@ struct DGCandProducer { // add histograms for the different process functions registry.add("reco/Stat", "Cut statistics; Selection criterion; Collisions", {HistType::kTH1F, {{14, -0.5, 13.5}}}); registry.add("reco/pt1Vspt2", "2 prong events, p_{T} versus p_{T}", {HistType::kTH2F, {{100, -3., 3.}, {100, -3., 3.0}}}); - registry.add("reco/TPCsignal1", "2 prong events, TPC signal of particle 1", {HistType::kTH2F, {{200, -3., 3.}, {200, 0., 100.0}}}); - registry.add("reco/TPCsignal2", "2 prong events, TPC signal of particle 2", {HistType::kTH2F, {{200, -3., 3.}, {200, 0., 100.0}}}); + registry.add("reco/TPCsignal1", "2 prong events, TPC signal versus p_{T} of particle 1", {HistType::kTH2F, {{200, -3., 3.}, {200, 0., 100.0}}}); + registry.add("reco/TPCsignal2", "2 prong events, TPC signal versus p_{T} of particle 2", {HistType::kTH2F, {{200, -3., 3.}, {200, 0., 100.0}}}); registry.add("reco/sig1VsSig2TPC", "2 prong events, TPC signal versus TPC signal", {HistType::kTH2F, {{100, 0., 100.}, {100, 0., 100.}}}); } @@ -208,7 +208,7 @@ struct DGCandProducer { } // produce TPC signal histograms for 2-track events - LOGF(info, "DG candidate: number of PV tracks %d", collision.numContrib()); + LOGF(debug, "DG candidate: number of PV tracks %d", collision.numContrib()); if (collision.numContrib() == 2) { auto cnt = 0; float pt1 = 0., pt2 = 0.; @@ -261,13 +261,8 @@ struct McDGCandProducer { "registry", {}}; - // this function properly updates UDMcCollisions and UDMcParticles and returns the value - // deltaIndex, which is needed to correct the McParticles indices - // For a given McCollision all associated McParticles are saved - template - void updateMcUDTables(TMcCollision const& mccol, - TMcParticles const& McParts, - int64_t& deltaIndex) + template + void updateUDMcCollisions(TMcCollision const& mccol) { // save mccol outputMcCollisions(mccol.bcId(), @@ -278,7 +273,36 @@ struct McDGCandProducer { mccol.t(), mccol.weight(), mccol.impactParameter()); + } + + template + void updateUDMcParticle(TMcParticle const& McPart, + int64_t& deltaIndex) + { + // save McPart + // mother and daughter indices are set to -1 + // ATTENTION: this can be improved to also include mother and daughter indices + std::vector newmids; + int32_t newdids[2] = {-1, -1}; + + // update UDMcParticles + outputMcParticles(outputMcCollisions.lastIndex(), + McPart.pdgCode(), + McPart.statusCode(), + McPart.flags(), + newmids, + newdids, + McPart.weight(), + McPart.px(), + McPart.py(), + McPart.pz(), + McPart.e()); + } + template + void updateUDMcParticles(TMcParticles const& McParts, + int64_t& deltaIndex) + { // save McParts // calculate conversion from old indices to new indices // old = mcpart.globalIndex() @@ -324,22 +348,93 @@ struct McDGCandProducer { } } + template + void updateUDMcTrackLabel(TTrack const& udtrack, int64_t const& deltaIndex) + { + // udtrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) + auto trackId = udtrack.trackId(); + if (trackId >= 0) { + auto track = udtrack.template track_as(); + auto mcTrackId = track.mcParticleId(); + if (mcTrackId >= 0) { + auto udMcTrackId = mcTrackId + deltaIndex; + outputMcTrackLabels(udMcTrackId, track.mcMask()); + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, -1); + } + } + + template + void updateUDMcTrackLabels(TTrack const& udtracks, int64_t const& deltaIndex) + { + // loop over all tracks + for (auto udtrack : udtracks) { + // udtrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) + auto trackId = udtrack.trackId(); + if (trackId >= 0) { + auto track = udtrack.template track_as(); + auto mcTrackId = track.mcParticleId(); + if (mcTrackId >= 0) { + auto udMcTrackId = mcTrackId + deltaIndex; + outputMcTrackLabels(udMcTrackId, track.mcMask()); + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, -1); + } + } + } + void init(InitContext& context) { // add histograms for the different process functions if (context.mOptions.get("processMC")) { registry.add("mcTruth/collisions", "Number of associated collisions", {HistType::kTH1F, {{11, -0.5, 10.5}}}); - registry.add("mcTruth/collType", "Collision type", {HistType::kTH1F, {{4, -0.5, 3.5}}}); + registry.add("mcTruth/collType", "Collision type", {HistType::kTH1F, {{5, -0.5, 4.5}}}); registry.add("mcTruth/IVMpt", "Invariant mass versus p_{T}", {HistType::kTH2F, {{150, 0.0, 3.0}, {150, 0.0, 3.0}}}); } } // process function for MC data - // save all GRANIITTI diffractive events and the MC truth of the DG events + // save the MC truth of all events of interest and of the DG events void processMC(aod::McCollisions const& mccols, aod::McParticles const& mcparts, UDCCs const& dgcands, UDTCs const& udtracks, CCs const& collisions, BCs const& bcs, TCs const& tracks) { + LOGF(info, "Number of McCollisions %d", mccols.size()); + LOGF(info, "Number of DG candidates %d", dgcands.size()); + LOGF(info, "Number of UD tracks %d", udtracks.size()); + + /* + example code for a list with dgcnad and mccol + + std::vector> qq; + qq.push_back(std::vector{1,1}); + qq.push_back(std::vector{2,2}); + qq.push_back(std::vector{3,-1}); + qq.push_back(std::vector{5,-1}); + qq.push_back(std::vector{7,9}); + qq.push_back(std::vector{-1,3}); + qq.push_back(std::vector{-1,4}); + qq.push_back(std::vector{-1,5}); + + std::sort(qq.begin(), qq.end(), + [](const std::vector& a, const std::vector& b) { + if (a[1] != -1 && b[1] != -1) { + return a[1] < b[1]; + } else { + return a[0] < b[0]; + } + }); + + for (auto q : qq) { + LOGF(info, " q[0] %d q [1] %d", q[0], q[1]); + } + */ // loop over McCollisions and UDCCs simultaneously auto mccol = mccols.iteratorAt(0); @@ -349,114 +444,124 @@ struct McDGCandProducer { int64_t lastSaved = -1; int64_t deltaIndex = 0; - int64_t firstIndex = 0, lastIndex = 0; - while (true) { - // determine the next dgcand with an associated collision - while (!dgcand.has_collision() && dgcand != lastdgcand) { - outputMcCollsLabels(-1); - dgcand++; - } + + // advance dgcand and mccol until both are AtEnd + int64_t mccolId = mccol.globalIndex(); + int64_t mcdgId = -1; + auto dgcandAtEnd = dgcand == lastdgcand; + auto mccolAtEnd = mccol == lastmccol; + bool goon = true; + while (goon) { + // check if dgcand has an associated McCollision if (!dgcand.has_collision()) { - // no dgcand left - outputMcCollsLabels(-1); - break; + mcdgId = -1; + } else { + auto dgcandCol = dgcand.collision_as(); + if (!dgcandCol.has_mcCollision()) { + mcdgId = -1; + } else { + mcdgId = dgcandCol.mcCollision().globalIndex(); + } } + LOGF(info, "\nstart of loop mcdgId %d mccolId %d", mcdgId, mccolId); - // related mc truth - auto dgcandCol = dgcand.collision_as(); - if (!dgcandCol.has_mcCollision()) { - // this collision has no MC truth - outputMcCollsLabels(-1); - continue; - } - auto mcdg = dgcandCol.mcCollision(); - auto dgmcId = mcdg.globalIndex(); - - // save also all GRANIITTI diffractive events - // keep UD Mc sorted according to AOD Mc - auto mccolId = mccol.globalIndex(); - while (mccolId <= dgmcId) { - auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccol.globalIndex()); - bool isGraniittiDiff = udhelpers::isGraniittiCDE(mcPartsSlice); - bool isPythiaDiff = udhelpers::isPythiaCDE(mcPartsSlice); - registry.get(HIST("mcTruth/collType"))->Fill(0., 1.); - registry.get(HIST("mcTruth/collType"))->Fill(1., (!isPythiaDiff && !isGraniittiDiff) * 1.); - registry.get(HIST("mcTruth/collType"))->Fill(2., isPythiaDiff * 1.); - registry.get(HIST("mcTruth/collType"))->Fill(3., isGraniittiDiff * 1.); - - if (isGraniittiDiff || isPythiaDiff) { - firstIndex = outputMcParticles.lastIndex() + 1; - updateMcUDTables(mccol, mcPartsSlice, deltaIndex); - lastSaved = mccolId; + // get dgcand tracks + auto dgTracks = udtracks.sliceByCached(aod::udtrack::udCollisionId, dgcand.globalIndex(), cache); - auto ivm = udhelpers::ivmGraniittiCDE(mcPartsSlice); - registry.get(HIST("mcTruth/IVMpt"))->Fill(ivm.M(), ivm.Perp()); - } - if (mccol == lastmccol) { - break; - } - mccol++; - mccolId = mccol.globalIndex(); - } + // two cases to consider + // 1. the event to process is a dgcand. In this case the Mc tables as well as the McLabel tables are updated + // 2. the event to process is an event of interest. In this case only the Mc tables are updated + if ((!dgcandAtEnd && !mccolAtEnd && (mcdgId <= mccolId)) || mccolAtEnd) { + // this is case 1. + LOGF(info, "doing case 1 with mcdgId %d", mcdgId); - // save the MC truth of the actual dgcand - // but check if this has not been added to the table yet - if (lastSaved != dgmcId) { - auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccol.globalIndex()); - firstIndex = outputMcParticles.lastIndex() + 1; - updateMcUDTables(mccol, mcPartsSlice, deltaIndex); - lastSaved = mccolId; + // update UDMcCollisions and UDMcColsLabels (for each UDCollision -> UDMcCollisions) + // update UDMcParticles and UDMcTrackLabels (for each UDTrack -> UDMcParticles) + + // If the dgcand has an associated McCollision then the McCollision and all associated + // McParticles are saved + if (mcdgId != lastSaved) { + if (mcdgId >= 0) { + LOGF(info, " saving McCollision %d", mcdgId); + // update UDMcCollisions + auto dgcandMcCol = dgcand.collision_as().mcCollision(); + updateUDMcCollisions(dgcandMcCol); + + // update UDMcColsLabels (for each UDCollision -> UDMcCollisions) + outputMcCollsLabels(outputMcCollisions.lastIndex()); + + // update lastSaved + lastSaved = mcdgId; + + // update UDMcParticles + auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mcdgId); + updateUDMcParticles(mcPartsSlice, deltaIndex); + + // update UDMcTrackLabels (for each UDTrack -> UDMcParticles) + updateUDMcTrackLabels(dgTracks, deltaIndex); - auto ivm = udhelpers::ivmGraniittiCDE(mcPartsSlice); - registry.get(HIST("mcTruth/IVMpt"))->Fill(ivm.M(), ivm.Perp()); - } - outputMcCollsLabels(outputMcCollisions.lastIndex()); - - // save the mclabels of the related tracks into outputMcTrackLabels - lastIndex = outputMcParticles.lastIndex(); - auto colTracks = udtracks.sliceByCached(aod::udtrack::udCollisionId, dgcand.globalIndex(), cache); - for (auto colTrack : colTracks) { - // colTrack (UDTCs) -> track (TCs) -> mcTrack (McParticles) -> udMcTrack (UDMcParticles) - auto trackId = colTrack.trackId(); - if (trackId >= 0) { - auto track = colTrack.track_as(); - auto mcTrackId = track.mcParticleId(); - if (mcTrackId >= 0) { - auto udMcTrackId = mcTrackId + deltaIndex; - outputMcTrackLabels(udMcTrackId, track.mcMask()); } else { - outputMcTrackLabels(-1, track.mcMask()); + // If the dgcand has no associated McCollision then only the McParticles which are associated + // with the tracks of the dgcand are saved + LOGF(info, " saving McCollision %d", -1); + + // update UDMcColsLabels (for each UDCollision -> UDMcCollisions) + outputMcCollsLabels(-1); + + // update UDMcParticles and UDMcTrackLabels (for each UDTrack -> UDMcParticles) + // loop over tracks of dgcand + for (auto dgtrack : dgTracks) { + if (dgtrack.has_track()) { + auto track = dgtrack.track_as(); + if (track.has_mcParticle()) { + auto mcPart = track.mcParticle(); + updateUDMcParticle(mcPart, deltaIndex); + updateUDMcTrackLabel(dgtrack, deltaIndex); + } else { + outputMcTrackLabels(-1, track.mcMask()); + } + } else { + outputMcTrackLabels(-1, -1); + } + } } + } + // advance dgcand + if (dgcand != lastdgcand) { + dgcand++; } else { - outputMcTrackLabels(-1, -1); + dgcandAtEnd = true; } - } + } else { + // this is case 2. + LOGF(info, "doing case 2"); - // next dg candidate - if (dgcand == lastdgcand) { - break; - } - dgcand++; - } + // update UDMcCollisions and UDMcParticles + if (mccolId != lastSaved) { + LOGF(info, " saving McCollision %d", mccolId); + + // update UDMcCollisions + updateUDMcCollisions(mccol); + + // update UDMcParticles + auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mcdgId); + updateUDMcParticles(mcPartsSlice, deltaIndex); - // save remaining GRANIITTI diffractive events - while (mccol != lastmccol) { - auto mcPartsSlice = mcparts.sliceBy(mcPartsPerMcCollision, mccol.globalIndex()); - bool isGraniittiDiff = udhelpers::isGraniittiCDE(mcPartsSlice); - bool isPythiaDiff = udhelpers::isPythiaCDE(mcPartsSlice); - registry.get(HIST("mcTruth/collType"))->Fill(0., 1.); - registry.get(HIST("mcTruth/collType"))->Fill(1., (!isPythiaDiff && !isGraniittiDiff) * 1.); - registry.get(HIST("mcTruth/collType"))->Fill(2., isPythiaDiff * 1.); - registry.get(HIST("mcTruth/collType"))->Fill(3., isGraniittiDiff * 1.); - - if (isGraniittiDiff || isPythiaDiff) { - updateMcUDTables(mccol, mcPartsSlice, deltaIndex); - - // update IVM versus pT - auto ivm = udhelpers::ivmGraniittiCDE(mcPartsSlice); - registry.get(HIST("mcTruth/IVMpt"))->Fill(ivm.M(), ivm.Perp()); + // update lastSaved + lastSaved = mccolId; + } + + // advance mccol + if (mccol != lastmccol) { + mccol++; + mccolId = mccol.globalIndex(); + } else { + mccolAtEnd = true; + } } - mccol++; + + goon = !dgcandAtEnd || !mccolAtEnd; + LOGF(info, "end of loop mcdgId %d mccolId %d", mcdgId, mccolId); } } PROCESS_SWITCH(McDGCandProducer, processMC, "Produce MC tables", false); diff --git a/PWGUD/Tasks/DGCandAnalyzer.cxx b/PWGUD/Tasks/DGCandAnalyzer.cxx index 8b328bd044d..9cabbc71a69 100644 --- a/PWGUD/Tasks/DGCandAnalyzer.cxx +++ b/PWGUD/Tasks/DGCandAnalyzer.cxx @@ -33,7 +33,7 @@ using namespace o2::framework::expressions; struct DGCandAnalyzer { // configurables - Configurable verbose{"Verbose", {}, "Additional print outs"}; + Configurable verbose{"Verbose", {}, "Additional printouts"}; Configurable candCaseSel{"CandCase", {}, "0: all Cands, 1: only ColCands,2: only BCCands"}; Configurable goodRunsFile{"goodRunsFile", {}, "json with list of good runs"}; @@ -144,6 +144,7 @@ struct DGCandAnalyzer { registry.add("stat/candCaseAll", "Types of all DG candidates", {HistType::kTH1F, {{5, -0.5, 4.5}}}); registry.add("stat/candCaseSel", "Types of all selectedDG candidates", {HistType::kTH1F, {{5, -0.5, 4.5}}}); registry.add("stat/nDGperRun", "Number of DG collisions per run", {HistType::kTH1D, {{1, 0, 1}}}); + registry.add("stat/nPVtracks", "Number of PV tracks of analyzed collisions", {HistType::kTH1D, {{51, -0.5, 50.5}}}); registry.add("tracks/nSigmaTPCPEl", "nSigma TPC for electrons", {HistType::kTH2F, {axispt, {100, -20.0, 20.0}}}); registry.add("tracks/nSigmaTPCPPi", "nSigma TPC for pions", {HistType::kTH2F, {axispt, {100, -20.0, 20.0}}}); @@ -223,6 +224,7 @@ struct DGCandAnalyzer { { // count collisions registry.fill(HIST("stat/candCaseAll"), 0., 1.); + registry.fill(HIST("stat/nPVtracks"), dgcand.numContrib(), 1.); // accept only selected run numbers int run = dgcand.runNumber(); @@ -239,6 +241,7 @@ struct DGCandAnalyzer { lastRun = run; LOGF(info, "done!"); } + registry.fill(HIST("stat/candCaseAll"), 1, 1.); // is BB bunch? auto bcnum = dgcand.globalBC(); @@ -257,11 +260,13 @@ struct DGCandAnalyzer { candCase = 2; } else if (dgcand.posX() == -2. && dgcand.posY() == 2. && dgcand.posZ() == -2.) { candCase = 3; + } else if (dgcand.posX() == -3. && dgcand.posY() == 3. && dgcand.posZ() == -3.) { + candCase = 3; } if (candCaseSel > 0 && candCase != candCaseSel) { return; } - registry.fill(HIST("stat/candCaseAll"), candCase, 1.); + registry.fill(HIST("stat/candCaseAll"), 2, 1.); // fill FIT amplitude histograms registry.fill(HIST("FIT/FT0AAmplitude"), dgcand.totalFT0AmplitudeA(), 1.); diff --git a/Tutorials/PWGUD/CMakeLists.txt b/Tutorials/PWGUD/CMakeLists.txt index ca6f47a141a..99f2d7b345a 100644 --- a/Tutorials/PWGUD/CMakeLists.txt +++ b/Tutorials/PWGUD/CMakeLists.txt @@ -23,3 +23,18 @@ o2physics_add_dpl_workflow(udtutorial-02b SOURCES UDTutorial_02b.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::DGPIDSelector COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-03a + SOURCES UDTutorial_03a.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-03b + SOURCES UDTutorial_03b.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(udtutorial-04 + SOURCES UDTutorial_04.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/Tutorials/PWGUD/UDTutorial_03a.cxx b/Tutorials/PWGUD/UDTutorial_03a.cxx new file mode 100644 index 00000000000..831abeb61c3 --- /dev/null +++ b/Tutorials/PWGUD/UDTutorial_03a.cxx @@ -0,0 +1,380 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// \brief UD tutorial +// \author Paul Buehler, paul.buehler@oeaw.ac.at +// \since October 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "TDatabasePDG.h" +#include "TLorentzVector.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/Core/UDHelpers.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UDTutorial03a { + SliceCache cache; + + // a pdg object + TDatabasePDG* pdg = nullptr; + + // get a DGCutparHolder + Configurable verbosity{"Verbosity", 0, "Determines level of verbosity"}; + + // initialize HistogramRegistry + HistogramRegistry registry{ + "registry", + {}}; + + // define abbreviations + using CCs = soa::Join; + using CC = CCs::iterator; + using TCs = soa::Join; + using TC = TCs::iterator; + + void init(InitContext& context) + { + // PDG + pdg = TDatabasePDG::Instance(); + + // add histograms for the different process functions + if (context.mOptions.get("processMCTruth")) { + registry.add("MC/Stat", "Count generated events; ; Entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}); + registry.add("MC/nParts", "Number of McParticles per collision; Number of McParticles; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("MC/genEtaPt", "Generated events; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/genRap", "Generated events; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/genMPt", "Generated events; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("MC/accEtaPt", "Generated events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/accRap", "Generated events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/accMPt", "Generated events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + } + + if (context.mOptions.get("processReco")) { + registry.add("Reco/Stat", "Count reconstruted events; ; Entries", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + registry.add("Reco/nTracks", "Number of reconstructed tracks per collision; Number of reconstructed tracks; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("Reco/nPVContributors", "Number of PV contributors per collision; Number of PV contributors; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("Reco/selEtaPt", "Selected events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("Reco/selRap", "Selected events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("Reco/selMPt", "Reconstructed events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + } + } + + Preslice partPerMcCollision = aod::mcparticle::mcCollisionId; + + // retrieve particle mass (GeV/c^2) from TDatabasePDG + float particleMass(TDatabasePDG* pdg, int pid) + { + auto mass = 0.; + TParticlePDG* pdgparticle = pdg->GetParticle(pid); + if (pdgparticle != nullptr) { + mass = pdgparticle->Mass(); + } + return mass; + } + + // check if a reconstructed track is a muon candidate + bool isMuonCandidate_rec(TC track) + { + if (abs(track.tpcNSigmaMu()) > 3.) { + return false; + } + return true; + } + + // check if a generated event is of the type rho0 -> mu+ + mu- using the MC particle stack + template + std::vector getDaughterParts_gen(MCTrack const& parts) + { + std::vector selectedParts; + + // in this case we expect the data files to contain events of the type rho0 -> mu+ + mu- + if (udhelpers::isSTARLightRhomumu(parts)) { + selectedParts.push_back(1); + selectedParts.push_back(2); + } + return selectedParts; + } + + // retrieve the two muon candidates of a given reconstructed collision + std::vector getDaughterTracks_rec(CC const& collision, TCs const& tracks) + { + // return a vector of track indices + std::vector emptySelection; + std::vector selectedTracks; + + // the collision must have exactly 2 PV tracks + if (collision.numContrib() != 2) { + return emptySelection; + } + + // the 2 PV tracks must have opposite charge signs + // and be muon candidates + int netCharge = 0; + int ind = -1; + for (auto track : tracks) { + ind++; + if (track.isPVContributor()) { + if (!isMuonCandidate_rec(track)) { + return emptySelection; + } + netCharge += track.sign(); + selectedTracks.push_back(ind); + } + } + if (netCharge != 0) { + return emptySelection; + } + + // the two PV contributors are both muon candidates + return selectedTracks; + } + + // compute the 2-track invariant mass using muon hypothesis + template + bool computeIVM_gen(TMcPart const& parts, std::vector partIds, TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + if (partIds.size() != 2) { + return false; + } + + // first particle + auto part1 = parts.iteratorAt(partIds.at(0)); + auto m1 = particleMass(pdg, 13); + auto ene1 = sqrt(pow(part1.px(), 2.) + pow(part1.py(), 2.) + pow(part1.pz(), 2.) + pow(m1, 2.)); + *lv1 = TLorentzVector(part1.px(), part1.py(), part1.pz(), ene1); + + // second particle + auto part2 = parts.iteratorAt(partIds.at(1)); + auto m2 = particleMass(pdg, 13); + auto ene2 = sqrt(pow(part2.px(), 2.) + pow(part2.py(), 2.) + pow(part2.pz(), 2.) + pow(m2, 2.)); + *lv2 = TLorentzVector(part2.px(), part2.py(), part2.pz(), ene2); + + // system + *lv = *lv1 + *lv2; + + return true; + } + + bool computeIVM_rec(TCs const& tracks, std::vector trackIds, TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + if (trackIds.size() != 2) { + return false; + } + + // first track + auto tr1 = tracks.iteratorAt(trackIds.at(0)); + auto m1 = particleMass(pdg, 13); + auto ene1 = sqrt(pow(tr1.px(), 2.) + pow(tr1.py(), 2.) + pow(tr1.pz(), 2.) + pow(m1, 2.)); + *lv1 = TLorentzVector(tr1.px(), tr1.py(), tr1.pz(), ene1); + + // second track + auto tr2 = tracks.iteratorAt(trackIds.at(1)); + auto m2 = particleMass(pdg, 13); + auto ene2 = sqrt(pow(tr2.px(), 2.) + pow(tr2.py(), 2.) + pow(tr2.pz(), 2.) + pow(m2, 2.)); + *lv2 = TLorentzVector(tr2.px(), tr2.py(), tr2.pz(), ene2); + + // system + *lv = *lv1 + *lv2; + + return true; + } + + // check a pair of particles to be in the acceptance of the detector + bool isInAcceptance(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + // the eta of the two muons to be in [-0.9, 0.9] + // the pT of the two muons to be in [0.1, infty] + // the y of the (mu+ + mu-)-system to be in [-0.9, 0.9] + + // check first track + if (abs(lv1->Eta()) > 0.9 || lv1->Pt() < 0.1) { + return false; + } + + // check second track + if (abs(lv2->Eta()) > 0.9 || lv2->Pt() < 0.1) { + return false; + } + + // check system + if (abs(lv->Rapidity()) > 0.9) { + return false; + } else { + return true; + } + } + + // check a reconstructed pair of tracks to be a candidate for an event of the type rho0 -> mu+ + mu- + bool isSelected_rec(TCs const& tracks, std::vector const& trackIds) + { + // tracks is expected to contain two tracks + if (trackIds.size() != 2) { + return false; + } + + // first track + auto tr1 = tracks.iteratorAt(trackIds.at(0)); + if (!tr1.isPVContributor()) { + return false; + } + + // second track + auto tr2 = tracks.iteratorAt(trackIds.at(1)); + if (!tr2.isPVContributor()) { + return false; + } + + // apply selection criteria + TLorentzVector* lv1 = new TLorentzVector(); + TLorentzVector* lv2 = new TLorentzVector(); + TLorentzVector* lv = new TLorentzVector(); + if (!computeIVM_rec(tracks, trackIds, lv1, lv2, lv)) { + return false; + } + // select generated events which are in the acceptance + // (!isInAcceptance(lv1, lv2, lv)) { + // return false; + //} + + // check tracks to be muon candidates + if (!isMuonCandidate_rec(tr1)) { + return false; + } + if (!isMuonCandidate_rec(tr2)) { + return false; + } + + // more selection criteria can be applied here ... + + // has passed all selection crtieria + return true; + } + + // ............................................................................................................... + void processMCTruth(aod::McCollisions const& mccollisions, CCs const& collisions, aod::McParticles const& McParts, TCs const& tracks) + { + // number of McCollisions in DF + if (verbosity > 0) { + LOGF(info, "Number of MC collisions %d", mccollisions.size()); + } + + // some variables + TLorentzVector* lv1_gen = new TLorentzVector(); + TLorentzVector* lv2_gen = new TLorentzVector(); + TLorentzVector* lv_gen = new TLorentzVector(); + + // loop over all genererated collisions + for (auto mccollision : mccollisions) { + registry.get(HIST("MC/Stat"))->Fill(0., 1.); + + // get McParticles which belong to mccollision + auto partSlice = McParts.sliceBy(partPerMcCollision, mccollision.globalIndex()); + registry.get(HIST("MC/nParts"))->Fill(partSlice.size(), 1.); + if (verbosity > 0) { + LOGF(info, "Number of McParts %d", partSlice.size()); + } + + // compute M versus Pt distributions for 2 cases + // 1. MC/genMPt - using all generated events, McTruth values + // 2. MC/accMPt - generated events which are within detector acceptance, McTruth values + + // select generated events of interest + // retrieve the index of the daughter particles in partSlice + // we expect there to be exactly 2 + auto daughterPartIds = getDaughterParts_gen(partSlice); + if (daughterPartIds.size() != 2) { + continue; + } + + // use the two selected McParticles to compute the invariant mass + // lv*_gen are TLorentzVectors of the two tracks and the sum + if (!computeIVM_gen(partSlice, daughterPartIds, lv1_gen, lv2_gen, lv_gen)) { + continue; + } + registry.get(HIST("MC/genEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("MC/genEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("MC/genRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("MC/genMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + if (verbosity > 0) { + LOGF(info, "IVMgen %f pT %f", lv_gen->M(), lv_gen->Pt()); + } + + // select generated events which are in the acceptance + if (!isInAcceptance(lv1_gen, lv2_gen, lv_gen)) { + continue; + } + registry.get(HIST("MC/accEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("MC/accEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("MC/accRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("MC/accMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + } + } + PROCESS_SWITCH(UDTutorial03a, processMCTruth, "Process MC truth", true); + + // ............................................................................................................... + void processReco(CC const& collision, TCs const& tracks, aod::McCollisions const& mccollisions, aod::McParticles const& McParts) + { + registry.get(HIST("Reco/Stat"))->Fill(0., 1.); + registry.get(HIST("Reco/nTracks"))->Fill(tracks.size(), 1.); + int nContributors = 0; + for (auto track : tracks) { + if (track.isPVContributor()) { + nContributors++; + } + } + registry.get(HIST("Reco/nPVContributors"))->Fill(nContributors, 1.); + + // some variables + TLorentzVector* lv1_rec = new TLorentzVector(); + TLorentzVector* lv2_rec = new TLorentzVector(); + TLorentzVector* lv_rec = new TLorentzVector(); + + // select 2 muon candidates using the reconstructed information + // MC truth is not considered in this step + auto daughterTrackIds = getDaughterTracks_rec(collision, tracks); + if (daughterTrackIds.size() != 2) { + if (verbosity > 0) { + LOGF(info, "Found %d daughter tracks.", daughterTrackIds.size()); + } + return; + } + registry.get(HIST("Reco/Stat"))->Fill(1., 1.); + + // apply all selection cuts on the selected reconstructed tracks + if (!isSelected_rec(tracks, daughterTrackIds)) { + return; + } + registry.get(HIST("Reco/Stat"))->Fill(2., 1.); + + // compute the invariant mass using the reconstructed information + if (!computeIVM_rec(tracks, daughterTrackIds, lv1_rec, lv2_rec, lv_rec)) { + return; + } + registry.get(HIST("Reco/selEtaPt"))->Fill(lv1_rec->Eta(), lv1_rec->Pt(), 1.); + registry.get(HIST("Reco/selEtaPt"))->Fill(lv2_rec->Eta(), lv2_rec->Pt(), 1.); + registry.get(HIST("Reco/selRap"))->Fill(lv_rec->Rapidity(), 1.); + registry.get(HIST("Reco/selMPt"))->Fill(lv_rec->M(), lv_rec->Pt(), 1.); + } + PROCESS_SWITCH(UDTutorial03a, processReco, "Process reconstructed data", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"udtutorial03a"})}; +} diff --git a/Tutorials/PWGUD/UDTutorial_03b.cxx b/Tutorials/PWGUD/UDTutorial_03b.cxx new file mode 100644 index 00000000000..58e8c1e26ed --- /dev/null +++ b/Tutorials/PWGUD/UDTutorial_03b.cxx @@ -0,0 +1,525 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// \brief UD tutorial +// \author Paul Buehler, paul.buehler@oeaw.ac.at +// \since October 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "TDatabasePDG.h" +#include "TLorentzVector.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/PIDResponse.h" +#include "PWGUD/Core/UDHelpers.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UDTutorial03b { + SliceCache cache; + + // a pdg object + TDatabasePDG* pdg = nullptr; + + // get a DGCutparHolder + Configurable verbosity{"Verbosity", 0, "Determines level of verbosity"}; + + // initialize HistogramRegistry + HistogramRegistry registry{ + "registry", + {}}; + + // define abbreviations + using CCs = soa::Join; + using CC = CCs::iterator; + using TCs = soa::Join; + using TC = TCs::iterator; + + void init(InitContext& context) + { + // PDG + pdg = TDatabasePDG::Instance(); + + // add histograms for the different process functions + if (context.mOptions.get("processMCTruth")) { + registry.add("MC/Stat", "Count generated events; ; Entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}); + registry.add("MC/recCols", "Number of reconstructed collisions; Number of reconstructed collisions; Entries", {HistType::kTH1F, {{31, -0.5, 30.5}}}); + registry.add("MC/nParts", "Number of McParticles per collision; Number of McParticles; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("MC/nRecTracks", "Number of reconstructed tracks per McParticle; Number of reconstructed tracks per McParticle; Entries", {HistType::kTH1F, {{11, -0.5, 10.5}}}); + registry.add("MC/genEtaPt", "Generated events; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/genRap", "Generated events; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/genMPt", "Generated events; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("MC/accEtaPt", "Generated events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/accRap", "Generated events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/accMPt", "Generated events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("MC/selEtaPt", "Selected events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/selRap", "Selected events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/selMPt", "Selected events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("MC/pDiff", "McTruth - reconstructed track momentum; McTruth - reconstructed track momentum; Entries", {HistType::kTH2F, {{240, -6., 6.}, {3, -1.5, 1.5}}}); + } + + if (context.mOptions.get("processReco")) { + registry.add("Reco/Stat", "Count reconstruted events; ; Entries", {HistType::kTH1F, {{5, -0.5, 4.5}}}); + registry.add("Reco/nTracks", "Number of reconstructed tracks per collision; Number of reconstructed tracks; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("Reco/nPVContributors", "Number of PV contributors per collision; Number of PV contributors; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("Reco/selEtaPt", "Selected events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("Reco/selRap", "Selected events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("Reco/selMPt", "Reconstructed events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("Reco/mcEtaPt", "Generated events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("Reco/mcRap", "Generated events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("Reco/mcMPt", "Generated events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("Reco/pDiff", "McTruth - reconstructed track momentum; McTruth - reconstructed track momentum; Entries", {HistType::kTH2F, {{240, -6., 6.}, {3, -1.5, 1.5}}}); + } + } + + Preslice partPerMcCollision = aod::mcparticle::mcCollisionId; + PresliceUnsorted colPerMcCollision = aod::mccollisionlabel::mcCollisionId; + PresliceUnsorted trackPerMcParticle = aod::mctracklabel::mcParticleId; + + // retrieve particle mass (GeV/c^2) from TDatabasePDG + float particleMass(TDatabasePDG* pdg, int pid) + { + auto mass = 0.; + TParticlePDG* pdgparticle = pdg->GetParticle(pid); + if (pdgparticle != nullptr) { + mass = pdgparticle->Mass(); + } + return mass; + } + + // check if a reconstructed track represents a muon candidate + bool isMuonCandidate_rec(TC track) + { + if (abs(track.tpcNSigmaMu()) > 3.) { + return false; + } + return true; + } + + // check if a generated event is of the type rho0 -> mu+ + mu- using the MC particle stack + template + std::vector getDaughterParts_gen(MCTrack const& parts) + { + std::vector selectedParts; + + // in this case we expect the data files to contain events of the type rho0 -> mu+ + mu- + if (udhelpers::isSTARLightRhomumu(parts)) { + selectedParts.push_back(1); + selectedParts.push_back(2); + } + return selectedParts; + } + + // find the McParticles belongin to given tracks + template + std::vector getDaughterParts_rec(TCs const& tracks, std::vector trackIds, MCTrack const& parts) + { + std::vector emptySelection; + std::vector selectedParts; + + for (auto trackId : trackIds) { + auto tr = tracks.iteratorAt(trackId); + if (!tr.has_mcParticle()) { + return emptySelection; + } else { + selectedParts.push_back(tr.mcParticle().globalIndex()); + } + } + return selectedParts; + } + + // retrieve the reconstructed tracks which are associated with the given McParticles + template + std::vector getDaughterTracks_gen(MCTrack const& parts, std::vector partIds, TCs const& tracks) + { + // return a vector of track indices + std::vector emptySelection; + std::vector selectedTracks; + + for (auto partId : partIds) { + auto part = parts.iteratorAt(partId); + auto trs = tracks.sliceBy(trackPerMcParticle, part.globalIndex()); + if (trs.size() == 0) { + return emptySelection; + } + if (trs.size() > 1) { + LOGF(info, "%d tracks belong to same McParticle!", trs.size()); + } + for (auto tr : trs) { + selectedTracks.push_back(tr.globalIndex()); + } + } + return selectedTracks; + } + + // retrieve the two muon candidates of a given reconstructed collision + std::vector getDaughterTracks_rec(CC const& collision, TCs const& tracks) + { + // return a vector of track indices + std::vector emptySelection; + std::vector selectedTracks; + + // the collision must have exactly 2 PV tracks + if (collision.numContrib() != 2) { + return emptySelection; + } + + // the 2 PV tracks must have opposite charge signs + // and be muon candidates + int netCharge = 0; + int ind = -1; + for (auto track : tracks) { + ind++; + if (track.isPVContributor()) { + if (!isMuonCandidate_rec(track)) { + return emptySelection; + } + netCharge += track.sign(); + selectedTracks.push_back(ind); + } + } + if (netCharge != 0) { + return emptySelection; + } + + // the two PV contributors are both muon candidates + return selectedTracks; + } + + // compute the 2-track invariant mass using muon hypothesis + template + bool computeIVM(TTrack const& tracks, std::vector trackIds, TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + if (trackIds.size() != 2) { + return false; + } + + // first track + auto tr1 = tracks.iteratorAt(trackIds.at(0)); + auto m1 = particleMass(pdg, 13); + auto ene1 = sqrt(pow(tr1.px(), 2.) + pow(tr1.py(), 2.) + pow(tr1.pz(), 2.) + pow(m1, 2.)); + *lv1 = TLorentzVector(tr1.px(), tr1.py(), tr1.pz(), ene1); + + // second track + auto tr2 = tracks.iteratorAt(trackIds.at(1)); + auto m2 = particleMass(pdg, 13); + auto ene2 = sqrt(pow(tr2.px(), 2.) + pow(tr2.py(), 2.) + pow(tr2.pz(), 2.) + pow(m2, 2.)); + *lv2 = TLorentzVector(tr2.px(), tr2.py(), tr2.pz(), ene2); + + // system + *lv = *lv1 + *lv2; + + return true; + } + + // check a pair of particles to be in the acceptance of the detector + bool isInAcceptance(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + // the eta of the two muons to be in [-0.9, 0.9] + // the pT of the two muons to be in [0.1, infty] + // the y of the (mu+ + mu-)-system to be in [-0.9, 0.9] + + // check first track + if (abs(lv1->Eta()) > 0.9 || lv1->Pt() < 0.1) { + return false; + } + + // check second track + if (abs(lv2->Eta()) > 0.9 || lv2->Pt() < 0.1) { + return false; + } + + // check system + if (abs(lv->Rapidity()) > 0.9) { + return false; + } else { + return true; + } + } + + // check a reconstructed pair of tracks to be a candidate for an event of the type rho0 -> mu+ + mu- + bool isSelected_rec(TCs const& tracks, std::vector const& trackIds) + { + // tracks is expected to contain two tracks + if (trackIds.size() != 2) { + return false; + } + + // first track + auto tr1 = tracks.iteratorAt(trackIds.at(0)); + if (!tr1.isPVContributor()) { + return false; + } + + // second track + auto tr2 = tracks.iteratorAt(trackIds.at(1)); + if (!tr2.isPVContributor()) { + return false; + } + + // apply selection criteria + TLorentzVector* lv1 = new TLorentzVector(); + TLorentzVector* lv2 = new TLorentzVector(); + TLorentzVector* lv = new TLorentzVector(); + if (!computeIVM(tracks, trackIds, lv1, lv2, lv)) { + return false; + } + // select generated events which are in the acceptance + // (!isInAcceptance(lv1, lv2, lv)) { + // return false; + //} + + // check tracks to be muon candidates + if (!isMuonCandidate_rec(tr1)) { + return false; + } + if (!isMuonCandidate_rec(tr2)) { + return false; + } + + // more selection criteria can be applied here ... + + // has passed all selection crtieria + return true; + } + + // ............................................................................................................... + void processMCTruth(aod::McCollisions const& mccollisions, CCs const& collisions, aod::McParticles const& McParts, TCs const& tracks) + { + // number of McCollisions in DF + if (verbosity > 0) { + LOGF(info, "Number of MC collisions %d", mccollisions.size()); + } + + // some variables + TLorentzVector* lv1_gen = new TLorentzVector(); + TLorentzVector* lv2_gen = new TLorentzVector(); + TLorentzVector* lv_gen = new TLorentzVector(); + TLorentzVector* lv1_rec = new TLorentzVector(); + TLorentzVector* lv2_rec = new TLorentzVector(); + TLorentzVector* lv_rec = new TLorentzVector(); + + // loop over all genererated collisions + for (auto mccollision : mccollisions) { + registry.get(HIST("MC/Stat"))->Fill(0., 1.); + + // get reconstructed collision which belongs to mccollision + auto colSlice = collisions.sliceBy(colPerMcCollision, mccollision.globalIndex()); + registry.get(HIST("MC/recCols"))->Fill(colSlice.size(), 1.); + + // get McParticles which belong to mccollision + auto partSlice = McParts.sliceBy(partPerMcCollision, mccollision.globalIndex()); + registry.get(HIST("MC/nParts"))->Fill(partSlice.size(), 1.); + if (verbosity > 0) { + LOGF(info, "Number of McParts %d", partSlice.size()); + } + + // compute M versus Pt distributions for 3 cases + // 1. MC/genMPt - using all generated events, McTruth values + // 2. MC/accMPt - generated events which are within detector acceptance, McTruth values + // 3. MC/selMPt - events which path selection criteria, reconstructed values + + // select generated events of interest + // retrieve the index of the daughter particles in partSlice + // we expect there to be exactly 2 + auto daughterPartIds = getDaughterParts_gen(partSlice); + if (daughterPartIds.size() != 2) { + continue; + } + + // use the two selected McParticles to compute the invariant mass + // lv*_gen are TLorentzVectors of the two tracks and the sum + if (!computeIVM(partSlice, daughterPartIds, lv1_gen, lv2_gen, lv_gen)) { + continue; + } + registry.get(HIST("MC/genEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("MC/genEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("MC/genRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("MC/genMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + if (verbosity > 0) { + LOGF(info, "IVMgen %f pT %f", lv_gen->M(), lv_gen->Pt()); + } + + // select generated events which are in the acceptance + if (!isInAcceptance(lv1_gen, lv2_gen, lv_gen)) { + continue; + } + registry.get(HIST("MC/accEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("MC/accEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("MC/accRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("MC/accMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + + // now obtain the reconstructed tracks + // which are the tracks associated with the McParticles + // we expect there to be exactly 2 + auto daughterTrackIds = getDaughterTracks_gen(partSlice, daughterPartIds, tracks); + if (daughterTrackIds.size() != 2) { + continue; + } + + // apply all selection cuts on the selected reconstructed tracks + if (!isSelected_rec(tracks, daughterTrackIds)) { + continue; + } + + // compute the invariant mass using the reconstructed tracks + if (!computeIVM(tracks, daughterTrackIds, lv1_rec, lv2_rec, lv_rec)) { + continue; + } + registry.get(HIST("MC/selEtaPt"))->Fill(lv1_rec->Eta(), lv1_rec->Pt(), 1.); + registry.get(HIST("MC/selEtaPt"))->Fill(lv2_rec->Eta(), lv2_rec->Pt(), 1.); + registry.get(HIST("MC/selRap"))->Fill(lv_rec->Rapidity(), 1.); + registry.get(HIST("MC/selMPt"))->Fill(lv_rec->M(), lv_rec->Pt(), 1.); + + // compute the difference between generated and reconstructed particle momentum + for (auto McPart : partSlice) { + // get track which corresponds to McPart + auto trackSlice = tracks.sliceBy(trackPerMcParticle, McPart.globalIndex()); + registry.get(HIST("MC/nRecTracks"))->Fill(trackSlice.size(), 1.); + + // are there reconstructed tracks? + if (trackSlice.size() > 0) { + for (auto track : trackSlice) { + auto pTrack = track.p(); + auto pPart = McPart.p(); + auto pDiff = pTrack - pPart; + registry.get(HIST("MC/pDiff"))->Fill(pDiff, track.isPVContributor(), 1.); + if (verbosity > 0) { + LOGF(info, " PID: %d Generated: %d Process: %d PV contributor: %d dP: %f", McPart.pdgCode(), McPart.producedByGenerator(), McPart.getProcess(), track.isPVContributor(), pDiff); + } + } + } else { + registry.get(HIST("MC/pDiff"))->Fill(-5.9, -1, 1.); + if (verbosity > 0) { + LOGF(info, " PID: %d Generated: %d Process: %d PV contributor: No dP: nan", McPart.pdgCode(), McPart.producedByGenerator(), McPart.getProcess()); + } + } + } + if (verbosity > 0) { + LOGF(info, ""); + } + } + } + PROCESS_SWITCH(UDTutorial03b, processMCTruth, "Process MC truth", true); + + // ............................................................................................................... + void processReco(CC const& collision, TCs const& tracks, aod::McCollisions const& mccollisions, aod::McParticles const& McParts) + { + registry.get(HIST("Reco/Stat"))->Fill(0., 1.); + registry.get(HIST("Reco/nTracks"))->Fill(tracks.size(), 1.); + int nContributors = 0; + for (auto track : tracks) { + if (track.isPVContributor()) { + nContributors++; + } + } + registry.get(HIST("Reco/nPVContributors"))->Fill(nContributors, 1.); + + // some variables + TLorentzVector* lv1_rec = new TLorentzVector(); + TLorentzVector* lv2_rec = new TLorentzVector(); + TLorentzVector* lv_rec = new TLorentzVector(); + TLorentzVector* lv1_gen = new TLorentzVector(); + TLorentzVector* lv2_gen = new TLorentzVector(); + TLorentzVector* lv_gen = new TLorentzVector(); + + // select 2 muon candidates using the reconstructed information + // MC truth is not considered in this step + auto daughterTrackIds = getDaughterTracks_rec(collision, tracks); + if (daughterTrackIds.size() != 2) { + if (verbosity > 0) { + LOGF(info, "Found %d daughter tracks.", daughterTrackIds.size()); + } + return; + } + registry.get(HIST("Reco/Stat"))->Fill(1., 1.); + + // apply all selection cuts on the selected reconstructed tracks + if (!isSelected_rec(tracks, daughterTrackIds)) { + return; + } + registry.get(HIST("Reco/Stat"))->Fill(2., 1.); + + // compute the invariant mass using the reconstructed information + if (!computeIVM(tracks, daughterTrackIds, lv1_rec, lv2_rec, lv_rec)) { + return; + } + registry.get(HIST("Reco/selEtaPt"))->Fill(lv1_rec->Eta(), lv1_rec->Pt(), 1.); + registry.get(HIST("Reco/selEtaPt"))->Fill(lv2_rec->Eta(), lv2_rec->Pt(), 1.); + registry.get(HIST("Reco/selRap"))->Fill(lv_rec->Rapidity(), 1.); + registry.get(HIST("Reco/selMPt"))->Fill(lv_rec->M(), lv_rec->Pt(), 1.); + + // now access the McTruth information + // get McCollision belonging to collision + if (collision.has_mcCollision()) { + // auto mccollision = collision.mcCollision(); + registry.get(HIST("Reco/Stat"))->Fill(3., 1.); + } else { + if (verbosity > 0) { + LOGF(info, "This collision has no associated McCollision"); + } + } + + // compute the difference between generated and reconstructed momentum + for (auto track : tracks) { + // is there an associated McParticle? + if (track.has_mcParticle()) { + auto pTrack = track.p(); + auto McPart = track.mcParticle(); + auto pPart = McPart.p(); + auto pDiff = pTrack - pPart; + registry.get(HIST("Reco/pDiff"))->Fill(pDiff, track.isPVContributor(), 1.); + if (verbosity > 0) { + LOGF(info, " PID: %d Generated: %d Process: %d PV contributor: %d dP: %f", McPart.pdgCode(), McPart.producedByGenerator(), McPart.getProcess(), track.isPVContributor(), pDiff); + } + } else { + registry.get(HIST("Reco/pDiff"))->Fill(-5.9, -1, 1.); + if (verbosity > 0) { + LOGF(info, " This track has no associated McParticle"); + } + } + } + if (verbosity > 0) { + LOGF(info, ""); + } + + // both muon candidates are required to have an associated McParticle + auto daughterPartIds = getDaughterParts_rec(tracks, daughterTrackIds, McParts); + if (daughterPartIds.size() != 2) { + if (verbosity > 0) { + LOGF(info, "Found %d daughter McParticles.", daughterPartIds.size()); + } + return; + } else { + registry.get(HIST("Reco/Stat"))->Fill(4., 1.); + + // compute invariant mass using McTruth + if (!computeIVM(McParts, daughterPartIds, lv1_gen, lv2_gen, lv_gen)) { + return; + } + registry.get(HIST("Reco/mcEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("Reco/mcEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("Reco/mcRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("Reco/mcMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + } + } + PROCESS_SWITCH(UDTutorial03b, processReco, "Process reconstructed data", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"udtutorial03b"})}; +} diff --git a/Tutorials/PWGUD/UDTutorial_04.cxx b/Tutorials/PWGUD/UDTutorial_04.cxx new file mode 100644 index 00000000000..4e4ea19a1b7 --- /dev/null +++ b/Tutorials/PWGUD/UDTutorial_04.cxx @@ -0,0 +1,478 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// \brief UD tutorial +// \author Paul Buehler, paul.buehler@oeaw.ac.at +// \since October 2023 + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" + +#include "TDatabasePDG.h" +#include "TLorentzVector.h" +#include "PWGUD/DataModel/UDTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +struct UDTutorial04 { + SliceCache cache; + + // a pdg object + TDatabasePDG* pdg = nullptr; + + // get a DGCutparHolder + Configurable verbosity{"Verbosity", 0, "Determines level of verbosity"}; + + // initialize HistogramRegistry + HistogramRegistry registry{ + "registry", + {}}; + + // define abbreviations + using CCs = soa::Join; + using CC = CCs::iterator; + using TCs = soa::Join; + using TC = TCs::iterator; + + void init(InitContext& context) + { + // PDG + pdg = TDatabasePDG::Instance(); + + // add histograms for the different process functions + if (context.mOptions.get("processMCTruth")) { + registry.add("MC/Stat", "Count generated events; ; Entries", {HistType::kTH1F, {{1, -0.5, 0.5}}}); + registry.add("MC/recCols", "Number of reconstructed collisions; Number of reconstructed collisions; Entries", {HistType::kTH1F, {{31, -0.5, 30.5}}}); + registry.add("MC/nParts", "Number of McParticles per collision; Number of McParticles; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("MC/nRecTracks", "Number of reconstructed tracks per McParticle; Number of reconstructed tracks per McParticle; Entries", {HistType::kTH1F, {{11, -0.5, 10.5}}}); + registry.add("MC/genEtaPt", "Generated events; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/genRap", "Generated events; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/genMPt", "Generated events; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("MC/accEtaPt", "Generated events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/accRap", "Generated events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/accMPt", "Generated events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + registry.add("MC/selEtaPt", "Selected events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("MC/selRap", "Selected events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("MC/selMPt", "Selected events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + + registry.add("Reco/nTracks", "Number of reconstructed tracks per collision; Number of reconstructed tracks; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("Reco/nPVContributors", "Number of PV contributors per collision; Number of PV contributors; Entries", {HistType::kTH1F, {{51, -0.5, 50.5}}}); + registry.add("Reco/selEtaPt", "Selected events in acceptance; eta (1); Pt (GeV/c)", {HistType::kTH2F, {{300, -1.5, 1.5}, {250, 0.0, 5.0}}}); + registry.add("Reco/selRap", "Selected events in acceptance; Rapidity (1)", {HistType::kTH1F, {{300, -1.5, 1.5}}}); + registry.add("Reco/selMPt", "Reconstructed events in acceptance; Mass (GeV/c^2); Pt (GeV/c)", {HistType::kTH2F, {{250, 2.5, 5.0}, {100, 0.0, 1.0}}}); + + registry.add("general/pDiff", "McTruth - reconstructed track momentum; McTruth - reconstructed track momentum; Entries", {HistType::kTH2F, {{240, -6., 6.}, {3, -1.5, 1.5}}}); + } + } + + // retrieve particle mass (GeV/c^2) from TDatabasePDG + float particleMass(TDatabasePDG* pdg, int pid) + { + auto mass = 0.; + TParticlePDG* pdgparticle = pdg->GetParticle(pid); + if (pdgparticle != nullptr) { + mass = pdgparticle->Mass(); + } + return mass; + } + + // check if a reconstructed track represents a muon candidate + bool isMuonCandidate_rec(TC track) + { + if (abs(track.tpcNSigmaMu()) > 3.) { + return false; + } + return true; + } + + // check if a generated event is of the type rho0 -> mu+ + mu- using the MC particle stack + template + std::vector getDaughterParts_gen(MCTrack const& parts) + { + std::vector selectedParts; + + // in this case we expect the data files to contain events of the type rho0 -> mu+ + mu- + // the first 3 particles in the particle stack are expected to have a pid of: 443013, +-(13, -13) + if (parts.size() > 3) { + if (parts.iteratorAt(0).pdgCode() == 443013 && abs(parts.iteratorAt(1).pdgCode()) == 13 && parts.iteratorAt(2).pdgCode() == -parts.iteratorAt(1).pdgCode()) { + selectedParts.push_back(1); + selectedParts.push_back(2); + } + } + return selectedParts; + } + + // find the McParticles belongin to given tracks + template + std::vector getDaughterParts_rec(TCs const& tracks, std::vector trackIds, MCTrack const& parts) + { + std::vector emptySelection; + std::vector selectedParts; + + for (auto trackId : trackIds) { + auto tr = tracks.iteratorAt(trackId); + if (!tr.has_udMcParticle()) { + return emptySelection; + } else { + selectedParts.push_back(tr.udMcParticle().globalIndex()); + } + } + return selectedParts; + } + + // retrieve the reconstructed tracks which are associated with the given McParticles + template + std::vector getDaughterTracks_gen(McPart const& parts, std::vector partIds, TCs const& tracks) + { + // it is expected that getDaughterParts_gen(parts) returns exactly 2 indices + std::vector emptySelection; + std::vector selectedTracks; + + for (auto partId : partIds) { + auto part = parts.iteratorAt(partId); + auto trs = tracks.sliceBy(trackPerMcParticle, part.globalIndex()); + if (trs.size() == 0) { + return emptySelection; + } + if (trs.size() > 1) { + LOGF(info, "%d tracks belong to same McParticle!", trs.size()); + } + for (auto tr : trs) { + selectedTracks.push_back(tr.globalIndex()); + } + } + return selectedTracks; + } + + // retrieve the two muon candidates of a given reconstructed collision + std::vector getDaughterTracks_rec(CC const& collision, TCs const& tracks) + { + // return a vector of track indices + std::vector emptySelection; + std::vector selectedTracks; + + // the collision must have exactly 2 PV tracks + if (collision.numContrib() != 2) { + return emptySelection; + } + + // the 2 PV tracks must have opposite charge signs + // and be muon candidates + int netCharge = 0; + int ind = -1; + for (auto track : tracks) { + ind++; + if (track.isPVContributor()) { + if (!isMuonCandidate_rec(track)) { + return emptySelection; + } + netCharge += track.sign(); + selectedTracks.push_back(ind); + } + } + if (netCharge != 0) { + return emptySelection; + } + + // the two PV contributors are both muon candidates + return selectedTracks; + } + + // compute the 2-track invariant mass using muon hypothesis + template + bool computeIVM(TTrack const& tracks, std::vector trackIds, TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + if (trackIds.size() != 2) { + return false; + } + + // first track + auto tr1 = tracks.iteratorAt(trackIds.at(0)); + auto m1 = particleMass(pdg, 13); + auto ene1 = sqrt(pow(tr1.px(), 2.) + pow(tr1.py(), 2.) + pow(tr1.pz(), 2.) + pow(m1, 2.)); + *lv1 = TLorentzVector(tr1.px(), tr1.py(), tr1.pz(), ene1); + + // second track + auto tr2 = tracks.iteratorAt(trackIds.at(1)); + auto m2 = particleMass(pdg, 13); + auto ene2 = sqrt(pow(tr2.px(), 2.) + pow(tr2.py(), 2.) + pow(tr2.pz(), 2.) + pow(m2, 2.)); + *lv2 = TLorentzVector(tr2.px(), tr2.py(), tr2.pz(), ene2); + + // system + *lv = *lv1 + *lv2; + + return true; + } + + // check a pair of particles to be in the acceptance of the detector + bool isInAcceptance(TLorentzVector* lv1, TLorentzVector* lv2, TLorentzVector* lv) + { + // the eta of the two muons to be in [-0.9, 0.9] + // the pT of the two muons to be in [0.1, infty] + // the y of the (mu+ + mu-)-system to be in [-0.9, 0.9] + + // check first track + if (abs(lv1->Eta()) > 0.9 || lv1->Pt() < 0.1) { + return false; + } + + // check second track + if (abs(lv2->Eta()) > 0.9 || lv2->Pt() < 0.1) { + return false; + } + + // check system + if (abs(lv->Rapidity()) > 0.9) { + return false; + } else { + return true; + } + } + + // check a reconstructed pair of tracks to be a candidate for an event of the type rho0 -> mu+ + mu- + template + bool isSelected_rec(TTrack const& tracks, std::vector const& trackIds) + { + // tracks is expected to contain two tracks + if (trackIds.size() != 2) { + return false; + } + + // first track + auto tr1 = tracks.iteratorAt(trackIds.at(0)); + if (!tr1.isPVContributor()) { + return false; + } + + // second track + auto tr2 = tracks.iteratorAt(trackIds.at(1)); + if (!tr2.isPVContributor()) { + return false; + } + + // apply selection criteria + TLorentzVector* lv1 = new TLorentzVector(); + TLorentzVector* lv2 = new TLorentzVector(); + TLorentzVector* lv = new TLorentzVector(); + if (!computeIVM(tracks, trackIds, lv1, lv2, lv)) { + return false; + } + // select generated events which are in the acceptance + if (!isInAcceptance(lv1, lv2, lv)) { + return false; + } + + // check tracks to be muon candidates + if (!isMuonCandidate_rec(tr1)) { + return false; + } + if (!isMuonCandidate_rec(tr2)) { + return false; + } + + // more selection criteria can be applied here ... + + // has passed all selection crtieria + return true; + } + + Preslice partPerMcCollision = aod::udmcparticle::udMcCollisionId; + PresliceUnsorted colPerMcCollision = aod::udcollision::udMcCollisionId; + PresliceUnsorted trackPerMcParticle = aod::udmctracklabel::udMcParticleId; + + // ............................................................................................................... + void processMCTruth(aod::UDMcCollisions const& mccollisions, CCs const& collisions, aod::UDMcParticles const& McParts, TCs const& tracks) + { + // number of McCollisions in DF + if (verbosity > 0) { + LOGF(info, "Number of MC collisions %d", mccollisions.size()); + } + + // some variables + TLorentzVector* lv1_gen = new TLorentzVector(); + TLorentzVector* lv2_gen = new TLorentzVector(); + TLorentzVector* lv_gen = new TLorentzVector(); + TLorentzVector* lv1_rec = new TLorentzVector(); + TLorentzVector* lv2_rec = new TLorentzVector(); + TLorentzVector* lv_rec = new TLorentzVector(); + + // loop over all generated collisions + for (auto mccollision : mccollisions) { + registry.get(HIST("MC/Stat"))->Fill(0., 1.); + + // get reconstructed collision which belongs to mccollision + auto colSlice = collisions.sliceBy(colPerMcCollision, mccollision.globalIndex()); + registry.get(HIST("MC/recCols"))->Fill(colSlice.size(), 1.); + + // get McParticles which belong to mccollision + auto partSlice = McParts.sliceBy(partPerMcCollision, mccollision.globalIndex()); + registry.get(HIST("MC/nParts"))->Fill(partSlice.size(), 1.); + if (verbosity > 0) { + LOGF(info, "Number of McParts %d", partSlice.size()); + } + + // loop over McParticles + for (auto McPart : partSlice) { + // get track which corresponds to McPart + auto trackSlice = tracks.sliceBy(trackPerMcParticle, McPart.globalIndex()); + registry.get(HIST("MC/nRecTracks"))->Fill(trackSlice.size(), 1.); + + // compute momentum difference between MCTruth and Reconstruction + if (trackSlice.size() > 0) { + for (auto track : trackSlice) { + auto pTrack = sqrt(track.px() * track.px() + track.py() * track.py() + track.pz() * track.pz()); + auto pPart = sqrt(McPart.px() * McPart.px() + McPart.py() * McPart.py() + McPart.pz() * McPart.pz()); + auto pDiff = pTrack - pPart; + registry.get(HIST("general/pDiff"))->Fill(pDiff, track.isPVContributor(), 1.); + if (verbosity > 0) { + LOGF(info, " PID: %d Generated: %d Process: %d PV contributor: %d dP: %f", McPart.pdgCode(), McPart.producedByGenerator(), McPart.getProcess(), track.isPVContributor(), pDiff); + } + } + } else { + registry.get(HIST("general/pDiff"))->Fill(-5.9, -1, 1.); + if (verbosity > 0) { + LOGF(info, " PID: %d Generated: %d Process: %d PV contributor: No dP: nan", McPart.pdgCode(), McPart.producedByGenerator(), McPart.getProcess()); + } + } + } + if (verbosity > 0) { + LOGF(info, ""); + } + + // compute M versus Pt distributions for 3 cases + // 1. MC/genMPt - using all generated events and McTruth values + // 2. MC/accMPt - generated events which are within detector acceptance and McTruth values + // 3. MC/selMPt - events which path selection criteria and reconstructed values + + // select generated events of interest + // retrieve the index of the daughter particles in partSlice + // we expect there to be exactly 2 + auto daughterPartIds = getDaughterParts_gen(partSlice); + if (daughterPartIds.size() != 2) { + continue; + } + + // use the two selected McParticles to compute the invariant mass + // lv*_gen are TLorentzVectors of the two tracks and the sum + if (!computeIVM(partSlice, daughterPartIds, lv1_gen, lv2_gen, lv_gen)) { + continue; + } + registry.get(HIST("MC/genEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("MC/genEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("MC/genRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("MC/genMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + if (verbosity > 0) { + LOGF(info, "IVMgen %f pT %f", lv_gen->M(), lv_gen->Pt()); + } + + // select generated events which are in the acceptance + if (!isInAcceptance(lv1_gen, lv2_gen, lv_gen)) { + continue; + } + registry.get(HIST("MC/accEtaPt"))->Fill(lv1_gen->Eta(), lv1_gen->Pt(), 1.); + registry.get(HIST("MC/accEtaPt"))->Fill(lv2_gen->Eta(), lv2_gen->Pt(), 1.); + registry.get(HIST("MC/accRap"))->Fill(lv_gen->Rapidity(), 1.); + registry.get(HIST("MC/accMPt"))->Fill(lv_gen->M(), lv_gen->Pt(), 1.); + + // now obtain the reconstructed tracks + // which are the tracks associated with the McParticles + auto daughterTrackIds = getDaughterTracks_gen(partSlice, daughterPartIds, tracks); + if (daughterTrackIds.size() != 2) { + continue; + } + + // apply all selection cuts on the selected reconstructed tracks + if (!isSelected_rec(tracks, daughterTrackIds)) { + continue; + } + + // compute the invariant mass using the reconstructed tracks + if (!computeIVM(tracks, daughterTrackIds, lv1_rec, lv2_rec, lv_rec)) { + continue; + } + registry.get(HIST("MC/selEtaPt"))->Fill(lv1_rec->Eta(), lv1_rec->Pt(), 1.); + registry.get(HIST("MC/selEtaPt"))->Fill(lv2_rec->Eta(), lv2_rec->Pt(), 1.); + registry.get(HIST("MC/selRap"))->Fill(lv_rec->Rapidity(), 1.); + registry.get(HIST("MC/selMPt"))->Fill(lv_rec->M(), lv_rec->Pt(), 1.); + + if (colSlice.size() == 1) { + LOGF(info, "This is a nice collision."); + LOGF(info, " McCollision: %d", mccollision.globalIndex()); + LOGF(info, " Collision: %d", colSlice.iteratorAt(0).globalIndex()); + } + } + } + PROCESS_SWITCH(UDTutorial04, processMCTruth, "Process MC truth", true); + + // ............................................................................................................... + void processReco(CC const& collision, TCs const& tracks, aod::UDMcCollisions const& mccollisions, aod::UDMcParticles const& McParts) + { + registry.get(HIST("Reco/nTracks"))->Fill(tracks.size(), 1.); + Partition PVContributors = aod::udtrack::isPVContributor == true; + PVContributors.bindTable(tracks); + registry.get(HIST("Reco/nPVContributors"))->Fill(PVContributors.size(), 1.); + + // some variables + TLorentzVector* lv1_rec = new TLorentzVector(); + TLorentzVector* lv2_rec = new TLorentzVector(); + TLorentzVector* lv_rec = new TLorentzVector(); + + // get McCollision belonging to collision + if (!collision.has_udMcCollision()) { + if (verbosity > 0) { + LOGF(info, "This collision has no associated McCollision"); + } + } + // auto mccollision = collision.udMcCollision(); + + // select 2 muon candidates using the reconstructed information + // MC truth is not considered in this step + auto daughterTrackIds = getDaughterTracks_rec(collision, tracks); + if (daughterTrackIds.size() != 2) { + if (verbosity > 0) { + LOGF(info, "Found %d daughter tracks.", daughterTrackIds.size()); + } + return; + } + + // both muon candidates are required to have an associated McParticle + auto daughterPartIds = getDaughterParts_rec(tracks, daughterTrackIds, McParts); + if (daughterPartIds.size() != 2) { + if (verbosity > 0) { + LOGF(info, "Found %d daughter McParticles.", daughterPartIds.size()); + } + return; + } + + // compute the invariant mass using the reconstructed information + if (!computeIVM(tracks, daughterTrackIds, lv1_rec, lv2_rec, lv_rec)) { + return; + } + + // apply all selection cuts on the selected reconstructed tracks + if (!isSelected_rec(tracks, daughterTrackIds)) { + return; + } + registry.get(HIST("Reco/selEtaPt"))->Fill(lv1_rec->Eta(), lv1_rec->Pt(), 1.); + registry.get(HIST("Reco/selEtaPt"))->Fill(lv2_rec->Eta(), lv2_rec->Pt(), 1.); + registry.get(HIST("Reco/selRap"))->Fill(lv_rec->Rapidity(), 1.); + registry.get(HIST("Reco/selMPt"))->Fill(lv_rec->M(), lv_rec->Pt(), 1.); + } + PROCESS_SWITCH(UDTutorial04, processReco, "Process reconstructed data", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"udtutorial04"})}; +}