From 0818ac236dce5ad69575eeaf62c54432ffc62e9e Mon Sep 17 00:00:00 2001 From: motomioya Date: Thu, 26 Oct 2023 05:44:43 +0900 Subject: [PATCH 1/9] Add search window for matching --- Common/TableProducer/mftmchMatchingML.cxx | 40 +++++++++++++---------- 1 file changed, 23 insertions(+), 17 deletions(-) diff --git a/Common/TableProducer/mftmchMatchingML.cxx b/Common/TableProducer/mftmchMatchingML.cxx index eb632e1a060..a6e1a9279dc 100644 --- a/Common/TableProducer/mftmchMatchingML.cxx +++ b/Common/TableProducer/mftmchMatchingML.cxx @@ -68,6 +68,8 @@ struct mftmchMatchingML { Configurable cfgModelDir{"ccdb-path", "Users/m/mooya/models", "base path to the ONNX models"}; Configurable cfgModelName{"ccdb-file", "model_LHC22o.onnx", "name of ONNX model file"}; Configurable cfgThrScore{"threshold-score", 0.5, "Threshold value for matching score"}; + Configurable cfgColWindow{"collision-window", 1, "Search window (collision ID) for MFT track"}; + Configurable cfgXYWindow{"XY-window", 3, "Search window (delta XY) for MFT track"}; Ort::Env env{ORT_LOGGING_LEVEL_WARNING, "model-explorer"}; Ort::SessionOptions session_options; @@ -159,7 +161,7 @@ struct mftmchMatchingML { std::vector input_tensor_values; input_tensor_values = getVariables(fwdtrack, mfttrack); - if (input_tensor_values[8] < 3) { + if (input_tensor_values[8] < cfgXYWindow) { std::vector input_tensors; input_tensors.push_back(Ort::Experimental::Value::CreateTensor(input_tensor_values.data(), input_tensor_values.size(), input_shape)); @@ -199,26 +201,30 @@ struct mftmchMatchingML { } } - void process(aod::Collisions::iterator const& collision, soa::Filtered const& fwdtracks, aod::MFTTracks const& mfttracks) + void process(aod::Collisions const& collisions, soa::Filtered const& fwdtracks, aod::MFTTracks const& mfttracks) { for (auto& [fwdtrack, mfttrack] : combinations(CombinationsFullIndexPolicy(fwdtracks, mfttracks))) { if (fwdtrack.trackType() == aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - double result = matchONNX(fwdtrack, mfttrack); - if (result > cfgThrScore) { - double mftchi2 = mfttrack.chi2(); - SMatrix5 mftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); - std::vector mftv1; - SMatrix55 mftcovs(mftv1.begin(), mftv1.end()); - o2::track::TrackParCovFwd mftpars1{mfttrack.z(), mftpars, mftcovs, mftchi2}; - mftpars1.propagateToZlinear(collision.posZ()); - - float dcaX = (mftpars1.getX() - collision.posX()); - float dcaY = (mftpars1.getY() - collision.posY()); - double px = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * cos(mfttrack.phi()); - double py = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi()); - double pz = fwdtrack.p() * cos(M_PI / 2 - atan(mfttrack.tgl())); - fwdtrackml(fwdtrack.collisionId(), 0, mfttrack.x(), mfttrack.y(), mfttrack.z(), mfttrack.phi(), mfttrack.tgl(), fwdtrack.sign() / std::sqrt(std::pow(px, 2) + std::pow(py, 2)), fwdtrack.nClusters(), -1, -1, -1, -1, -1, result, mfttrack.globalIndex(), fwdtrack.globalIndex(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mfttrack.trackTime(), mfttrack.trackTimeRes(), mfttrack.eta(), std::sqrt(std::pow(px, 2) + std::pow(py, 2)), std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)), dcaX, dcaY); + if (fwdtrack.has_collision() && mfttrack.has_collision()){ + if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow){ + double result = matchONNX(fwdtrack, mfttrack); + if (result > cfgThrScore) { + double mftchi2 = mfttrack.chi2(); + SMatrix5 mftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); + std::vector mftv1; + SMatrix55 mftcovs(mftv1.begin(), mftv1.end()); + o2::track::TrackParCovFwd mftpars1{mfttrack.z(), mftpars, mftcovs, mftchi2}; + mftpars1.propagateToZlinear(mfttrack.collision().posZ()); + + float dcaX = (mftpars1.getX() - mfttrack.collision().posX()); + float dcaY = (mftpars1.getY() - mfttrack.collision().posY()); + double px = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * cos(mfttrack.phi()); + double py = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi()); + double pz = fwdtrack.p() * cos(M_PI / 2 - atan(mfttrack.tgl())); + fwdtrackml(fwdtrack.collisionId(), 0, mfttrack.x(), mfttrack.y(), mfttrack.z(), mfttrack.phi(), mfttrack.tgl(), fwdtrack.sign() / std::sqrt(std::pow(px, 2) + std::pow(py, 2)), fwdtrack.nClusters(), fwdtrack.pDca(), fwdtrack.rAtAbsorberEnd(), 0, 0, 0, result, mfttrack.globalIndex(), fwdtrack.globalIndex(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mfttrack.trackTime(), mfttrack.trackTimeRes(), mfttrack.eta(), std::sqrt(std::pow(px, 2) + std::pow(py, 2)), std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)), dcaX, dcaY); + } + } } } } From 229096cd6cee48e5e054daf0dbfd84d4785f0335 Mon Sep 17 00:00:00 2001 From: motomioya Date: Thu, 26 Oct 2023 05:46:08 +0900 Subject: [PATCH 2/9] Fixed bugs in filter --- Common/TableProducer/mftmchMatchingML.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/TableProducer/mftmchMatchingML.cxx b/Common/TableProducer/mftmchMatchingML.cxx index a6e1a9279dc..9c69590a213 100644 --- a/Common/TableProducer/mftmchMatchingML.cxx +++ b/Common/TableProducer/mftmchMatchingML.cxx @@ -59,7 +59,7 @@ struct mftmchMatchingML { float chi2up = 1000000; float chi2MatchMCHMIDup = 1000000; - Filter etaFilter = ((etalow < aod::fwdtrack::eta) && (etaup < aod::fwdtrack::eta)); + Filter etaFilter = ((etalow < aod::fwdtrack::eta) && (aod::fwdtrack::eta < etaup)); Filter pDcaFilter = (((pDCAcutrAtBsorberEndlow1 < aod::fwdtrack::rAtAbsorberEnd) && (aod::fwdtrack::rAtAbsorberEnd < pDCAcutrAtBsorberEndup1) && (aod::fwdtrack::pDca < pDCAcutdcaup1)) || ((pDCAcutrAtBsorberEndlow2 < aod::fwdtrack::rAtAbsorberEnd) && (aod::fwdtrack::rAtAbsorberEnd < pDCAcutrAtBsorberEndup2) && (aod::fwdtrack::pDca < pDCAcutdcaup2))); Filter chi2Filter = (aod::fwdtrack::chi2 < chi2up); Filter chi2MatchFilter = (aod::fwdtrack::chi2MatchMCHMID < chi2MatchMCHMIDup); From b61e43134f6fa45013d69e51368ee52704a39ccb Mon Sep 17 00:00:00 2001 From: motomioya Date: Thu, 26 Oct 2023 15:59:08 +0900 Subject: [PATCH 3/9] Modify for formatting --- Common/TableProducer/mftmchMatchingML.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Common/TableProducer/mftmchMatchingML.cxx b/Common/TableProducer/mftmchMatchingML.cxx index 9c69590a213..f828dad2af0 100644 --- a/Common/TableProducer/mftmchMatchingML.cxx +++ b/Common/TableProducer/mftmchMatchingML.cxx @@ -206,8 +206,8 @@ struct mftmchMatchingML { for (auto& [fwdtrack, mfttrack] : combinations(CombinationsFullIndexPolicy(fwdtracks, mfttracks))) { if (fwdtrack.trackType() == aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - if (fwdtrack.has_collision() && mfttrack.has_collision()){ - if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow){ + if (fwdtrack.has_collision() && mfttrack.has_collision()) { + if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow) { double result = matchONNX(fwdtrack, mfttrack); if (result > cfgThrScore) { double mftchi2 = mfttrack.chi2(); From 5d93541ae7e8bdb07559b892ab54094e77f5ae52 Mon Sep 17 00:00:00 2001 From: motomioya Date: Thu, 30 Nov 2023 22:51:12 +0900 Subject: [PATCH 4/9] Fixed to select one mfttrack per fwdtrack --- Common/TableProducer/mftmchMatchingML.cxx | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/Common/TableProducer/mftmchMatchingML.cxx b/Common/TableProducer/mftmchMatchingML.cxx index f828dad2af0..6320ec7b625 100644 --- a/Common/TableProducer/mftmchMatchingML.cxx +++ b/Common/TableProducer/mftmchMatchingML.cxx @@ -203,13 +203,24 @@ struct mftmchMatchingML { void process(aod::Collisions const& collisions, soa::Filtered const& fwdtracks, aod::MFTTracks const& mfttracks) { - for (auto& [fwdtrack, mfttrack] : combinations(CombinationsFullIndexPolicy(fwdtracks, mfttracks))) { - + for (auto& fwdtrack : fwdtracks) { if (fwdtrack.trackType() == aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - if (fwdtrack.has_collision() && mfttrack.has_collision()) { - if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow) { - double result = matchONNX(fwdtrack, mfttrack); - if (result > cfgThrScore) { + double bestscore = 0; + int bestmfttrackid = -1; + for (auto& mfttrack : mfttracks) { + if (fwdtrack.has_collision() && mfttrack.has_collision()) { + if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow) { + double result = matchONNX(fwdtrack, mfttrack); + if (result > cfgThrScore) { + bestscore = result; + bestmfttrackid = mfttrack.globalIndex(); + } + } + } + } + if (bestmfttrackid != -1) { + for (auto& mfttrack : mfttracks) { + if (mfttrack.globalIndex() == bestmfttrackid) { double mftchi2 = mfttrack.chi2(); SMatrix5 mftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); std::vector mftv1; From 4eed1aecd7753da9068bfcac5ed969c1f3cdf6af Mon Sep 17 00:00:00 2001 From: motomioya Date: Fri, 1 Dec 2023 14:28:12 +0900 Subject: [PATCH 5/9] Fixed MCH track to match single MFT track --- Common/TableProducer/mftmchMatchingML.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/TableProducer/mftmchMatchingML.cxx b/Common/TableProducer/mftmchMatchingML.cxx index 6320ec7b625..c8e66d4a247 100644 --- a/Common/TableProducer/mftmchMatchingML.cxx +++ b/Common/TableProducer/mftmchMatchingML.cxx @@ -233,7 +233,7 @@ struct mftmchMatchingML { double px = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * cos(mfttrack.phi()); double py = fwdtrack.p() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi()); double pz = fwdtrack.p() * cos(M_PI / 2 - atan(mfttrack.tgl())); - fwdtrackml(fwdtrack.collisionId(), 0, mfttrack.x(), mfttrack.y(), mfttrack.z(), mfttrack.phi(), mfttrack.tgl(), fwdtrack.sign() / std::sqrt(std::pow(px, 2) + std::pow(py, 2)), fwdtrack.nClusters(), fwdtrack.pDca(), fwdtrack.rAtAbsorberEnd(), 0, 0, 0, result, mfttrack.globalIndex(), fwdtrack.globalIndex(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mfttrack.trackTime(), mfttrack.trackTimeRes(), mfttrack.eta(), std::sqrt(std::pow(px, 2) + std::pow(py, 2)), std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)), dcaX, dcaY); + fwdtrackml(fwdtrack.collisionId(), 0, mfttrack.x(), mfttrack.y(), mfttrack.z(), mfttrack.phi(), mfttrack.tgl(), fwdtrack.sign() / std::sqrt(std::pow(px, 2) + std::pow(py, 2)), fwdtrack.nClusters(), fwdtrack.pDca(), fwdtrack.rAtAbsorberEnd(), 0, 0, 0, bestscore, mfttrack.globalIndex(), fwdtrack.globalIndex(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mfttrack.trackTime(), mfttrack.trackTimeRes(), mfttrack.eta(), std::sqrt(std::pow(px, 2) + std::pow(py, 2)), std::sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2)), dcaX, dcaY); } } } From 139244799ba45c8fd8dace304dec83c8679fde45 Mon Sep 17 00:00:00 2001 From: motomioya Date: Fri, 1 Dec 2023 14:54:42 +0900 Subject: [PATCH 6/9] Modified white space --- Common/TableProducer/mftmchMatchingML.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/TableProducer/mftmchMatchingML.cxx b/Common/TableProducer/mftmchMatchingML.cxx index c8e66d4a247..13d4140156e 100644 --- a/Common/TableProducer/mftmchMatchingML.cxx +++ b/Common/TableProducer/mftmchMatchingML.cxx @@ -207,7 +207,7 @@ struct mftmchMatchingML { if (fwdtrack.trackType() == aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { double bestscore = 0; int bestmfttrackid = -1; - for (auto& mfttrack : mfttracks) { + for (auto& mfttrack : mfttracks) { if (fwdtrack.has_collision() && mfttrack.has_collision()) { if (0 <= fwdtrack.collisionId() - mfttrack.collisionId() && fwdtrack.collisionId() - mfttrack.collisionId() < cfgColWindow) { double result = matchONNX(fwdtrack, mfttrack); From 81a7c67727172352a23b8b07e55f2f4d594f6f15 Mon Sep 17 00:00:00 2001 From: motomioya Date: Wed, 8 May 2024 16:07:07 +0900 Subject: [PATCH 7/9] Add process to associated muons --- PWGDQ/TableProducer/tableMakerMC.cxx | 580 ++++++++++++++++++++++++++- 1 file changed, 579 insertions(+), 1 deletion(-) diff --git a/PWGDQ/TableProducer/tableMakerMC.cxx b/PWGDQ/TableProducer/tableMakerMC.cxx index 7e73247701b..4a900066b3f 100644 --- a/PWGDQ/TableProducer/tableMakerMC.cxx +++ b/PWGDQ/TableProducer/tableMakerMC.cxx @@ -39,6 +39,7 @@ #include "PWGDQ/Core/MCSignalLibrary.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/CollisionAssociationTables.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" #include "Field/MagneticField.h" @@ -75,6 +76,8 @@ using MyBarrelTracksWithCov = soa::Join; using MyMuons = soa::Join; using MyMuonsWithCov = soa::Join; +using MyMuonsColl = soa::Join; +using MyMuonsCollWithCov = soa::Join; using MyEvents = soa::Join; using MyEventsWithMults = soa::Join; @@ -227,7 +230,9 @@ struct TableMakerMC { context.mOptions.get("processBarrelOnlyWithMults") || context.mOptions.get("processAmbiguousBarrelOnly")); bool enableMuonHistos = (context.mOptions.get("processFull") || context.mOptions.get("processFullWithCov") || context.mOptions.get("processMuonOnlyWithCent") || context.mOptions.get("processMuonOnlyWithCov") || - context.mOptions.get("processAmbiguousMuonOnlyWithCov") || context.mOptions.get("processAmbiguousMuonOnly")); + context.mOptions.get("processAmbiguousMuonOnlyWithCov") || context.mOptions.get("processAmbiguousMuonOnly") || + context.mOptions.get("processAssociatedMuonOnly") || context.mOptions.get("processAssociatedMuonOnlyWithCov") || + context.mOptions.get("processAssociatedMuonOnlyWithCovAndCent") || context.mOptions.get("processAssociatedMuonOnlyWithCovAndMults"); // TODO: switch on/off histogram classes depending on which process function we run if (enableBarrelHistos) { if (fDoDetailedQA) { @@ -302,6 +307,8 @@ struct TableMakerMC { Preslice perMcCollision = aod::mcparticle::mcCollisionId; Preslice perCollisionTracks = aod::track::collisionId; Preslice perCollisionMuons = aod::fwdtrack::collisionId; + Preslice trackIndicesPerCollision = aod::track_association::collisionId; + Preslice fwdtrackIndicesPerCollision = aod::track_association::collisionId; // Templated function instantianed for all of the process functions template @@ -877,6 +884,549 @@ struct TableMakerMC { fEventLabels.clear(); } + // Templated function instantianed for all of the process functions + template + void fullSkimmingIndices(TEvent const& collisions, aod::BCsWithTimestamps const&, TTracks const& tracksBarrel, TMuons const& tracksMuon, aod::McCollisions const& /*mcEvents*/, aod::McParticles_001 const& mcTracks, AssocTracks const& trackIndices, AssocMuons const& fwdtrackIndices) + { + // Loop over collisions and produce skimmed data tables for: + // 1) all selected collisions + // 2) all MC collisions pointed to by selected collisions + // 2.1) MC collision labels + // 3) all selected tracks + // 4) MC track labels + // Maps of indices are stored for all selected MC tracks (selected MC signals + MC truth particles for selected tracks) + // The MC particles table is generated in a separate loop over McParticles, after the indices for all MC particles we want + // to store are fed into the ordered std maps + + // temporary variables used for the indexing of the skimmed MC stack + std::map fNewLabels; + std::map fNewLabelsReversed; + std::map fMCFlags; + std::map fEventIdx; + std::map fEventLabels; + int fCounters[2] = {0, 0}; //! [0] - particle counter, [1] - event counter + + uint16_t mcflags = 0; + uint64_t trackFilteringTag = 0; + uint8_t trackTempFilterMap = 0; + + for (auto& collision : collisions) { + // TODO: investigate the collisions without corresponding mcCollision + if (!collision.has_mcCollision()) { + continue; + } + auto bc = collision.template bc_as(); + if (fCurrentRun != bc.runNumber()) { + if (fIsRun2 == true) { + grpmagrun2 = fCCDB->getForTimeStamp(grpmagPathRun2, bc.timestamp()); + if (grpmagrun2 != nullptr) { + o2::base::Propagator::initFieldFromGRP(grpmagrun2); + } + } else { + grpmag = fCCDB->getForTimeStamp(grpmagPath, bc.timestamp()); + if (grpmag != nullptr) { + o2::base::Propagator::initFieldFromGRP(grpmag); + } + if (fPropMuon) { + VarManager::SetupMuonMagField(); + } + } + fCurrentRun = bc.runNumber(); + } + + // get the trigger aliases + uint32_t triggerAliases = collision.alias_raw(); + // store the selection decisions + uint64_t tag = collision.selection_raw(); + if (collision.sel7()) { + tag |= (uint64_t(1) << evsel::kNsel); //! SEL7 stored at position kNsel in the tag bit map + } + + auto mcCollision = collision.mcCollision(); + VarManager::ResetValues(0, VarManager::kNEventWiseVariables); + VarManager::fgValues[VarManager::kRunNo] = bc.runNumber(); + VarManager::fgValues[VarManager::kBC] = bc.globalBC(); + VarManager::fgValues[VarManager::kTimestamp] = bc.timestamp(); + VarManager::FillEvent(collision); // extract event information and place it in the fValues array + VarManager::FillEvent(mcCollision); + + if (fDoDetailedQA) { + fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); + } + // fill stats information, before selections + for (int i = 0; i < kNaliases; i++) { + if (triggerAliases & (uint32_t(1) << i)) { + (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(i)); + } + } + (reinterpret_cast(fStatsList->At(0)))->Fill(2.0, static_cast(kNaliases)); + + if (!fEventCut->IsSelected(VarManager::fgValues)) { + continue; + } + + if (fConfigQA) { + fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + } + + // fill stats information, after selections + for (int i = 0; i < kNaliases; i++) { + if (triggerAliases & (uint32_t(1) << i)) { + (reinterpret_cast(fStatsList->At(0)))->Fill(3.0, static_cast(i)); + } + } + (reinterpret_cast(fStatsList->At(0)))->Fill(3.0, static_cast(kNaliases)); + + event(tag, bc.runNumber(), collision.posX(), collision.posY(), collision.posZ(), collision.numContrib(), collision.collisionTime(), collision.collisionTimeRes()); + if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionMult) > 0 && (TEventFillMap & VarManager::ObjTypes::CollisionCent) > 0) { + eventExtended(collision.bc().globalBC(), collision.bc().triggerMask(), 0, triggerAliases, VarManager::fgValues[VarManager::kCentVZERO], + collision.multTPC(), collision.multFV0A(), collision.multFV0C(), collision.multFT0A(), collision.multFT0C(), + collision.multFDDA(), collision.multFDDC(), collision.multZNA(), collision.multZNC(), collision.multTracklets(), collision.multNTracksPV(), + collision.centFT0C()); + } else if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionMult) > 0) { + eventExtended(bc.globalBC(), bc.triggerMask(), bc.timestamp(), triggerAliases, VarManager::fgValues[VarManager::kCentVZERO], + collision.multTPC(), collision.multFV0A(), collision.multFV0C(), collision.multFT0A(), collision.multFT0C(), + collision.multFDDA(), collision.multFDDC(), collision.multZNA(), collision.multZNC(), collision.multTracklets(), collision.multNTracksPV(), + -1); + } else if constexpr ((TEventFillMap & VarManager::ObjTypes::CollisionCent) > 0) { + eventExtended(bc.globalBC(), bc.triggerMask(), bc.timestamp(), triggerAliases, VarManager::fgValues[VarManager::kCentVZERO], + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, collision.centFT0C()); + } else { + eventExtended(bc.globalBC(), bc.triggerMask(), bc.timestamp(), triggerAliases, VarManager::fgValues[VarManager::kCentVZERO], -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1); + } + eventVtxCov(collision.covXX(), collision.covXY(), collision.covXZ(), collision.covYY(), collision.covYZ(), collision.covZZ(), collision.chi2()); + // make an entry for this MC event only if it was not already added to the table + if (!(fEventLabels.find(mcCollision.globalIndex()) != fEventLabels.end())) { + eventMC(mcCollision.generatorsID(), mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), + mcCollision.t(), mcCollision.weight(), mcCollision.impactParameter()); + fEventLabels[mcCollision.globalIndex()] = fCounters[1]; + fCounters[1]++; + } + eventMClabels(fEventLabels.find(mcCollision.globalIndex())->second, collision.mcMask()); + + // loop over the MC truth tracks and find those that need to be written + auto groupedMcTracks = mcTracks.sliceBy(perMcCollision, mcCollision.globalIndex()); + for (auto& mctrack : groupedMcTracks) { + // check all the requested MC signals and fill a decision bit map + mcflags = 0; + int i = 0; + for (auto& sig : fMCSignals) { + bool checked = false; + if constexpr (soa::is_soa_filtered_v) { + auto mctrack_raw = groupedMcTracks.rawIteratorAt(mctrack.globalIndex()); + checked = sig.CheckSignal(false, mctrack_raw); + } else { + checked = sig.CheckSignal(false, mctrack); + } + if (checked) { + mcflags |= (uint16_t(1) << i); + } + i++; + } + if (mcflags == 0) { + continue; + } + + if (!(fNewLabels.find(mctrack.globalIndex()) != fNewLabels.end())) { + fNewLabels[mctrack.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mctrack.globalIndex(); + fMCFlags[mctrack.globalIndex()] = mcflags; + fEventIdx[mctrack.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; + fCounters[0]++; + + // if any of the MC signals was matched, then fill histograms and write that MC particle into the new stack + // fill histograms for each of the signals, if found + if (fConfigQA) { + VarManager::FillTrackMC(mcTracks, mctrack); + int j = 0; + for (auto signal = fMCSignals.begin(); signal != fMCSignals.end(); signal++, j++) { + if (mcflags & (uint16_t(1) << j)) { + fHistMan->FillHistClass(Form("MCTruth_%s", (*signal).GetName()), VarManager::fgValues); + } + } + } + } + } // end loop over mc stack + + int isAmbiguous = 0; + // loop over reconstructed tracks + if constexpr (static_cast(TTrackFillMap)) { + trackBasic.reserve(tracksBarrel.size()); + trackBarrel.reserve(tracksBarrel.size()); + trackBarrelPID.reserve(tracksBarrel.size()); + trackBarrelLabels.reserve(tracksBarrel.size()); + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::TrackCov)) { + trackBarrelCov.reserve(tracksBarrel.size()); + } + + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); + // loop over tracks + for (auto& trackId : trackIdsThisCollision) { + auto track = trackId.template track_as(); + if constexpr ((TTrackFillMap & VarManager::ObjTypes::AmbiTrack) > 0) { + if (fIsAmbiguous) { + isAmbiguous = (track.compatibleCollIds().size() != 1); + } + } + trackFilteringTag = uint64_t(0); + trackTempFilterMap = uint8_t(0); + VarManager::FillTrack(track); + // If no MC particle is found, skip the track + if (!track.has_mcParticle()) { + continue; + } + auto mctrack = track.template mcParticle_as(); + VarManager::FillTrackMC(mcTracks, mctrack); + + if (fDoDetailedQA) { + fHistMan->FillHistClass("TrackBarrel_BeforeCuts", VarManager::fgValues); + if (fIsAmbiguous && isAmbiguous == 1) { + fHistMan->FillHistClass("Ambiguous_TrackBarrel_BeforeCuts", VarManager::fgValues); + } + } + // apply track cuts and fill stats histogram + int i = 0; + for (auto& cut : fTrackCuts) { + if (cut.IsSelected(VarManager::fgValues)) { + trackTempFilterMap |= (uint8_t(1) << i); + if (fConfigQA) { + fHistMan->FillHistClass(Form("TrackBarrel_%s", cut.GetName()), VarManager::fgValues); // fill the reconstructed truth + if (fIsAmbiguous && isAmbiguous == 1) { + fHistMan->FillHistClass(Form("Ambiguous_TrackBarrel_%s", cut.GetName()), VarManager::fgValues); + } + } + (reinterpret_cast(fStatsList->At(1)))->Fill(static_cast(i)); + } + i++; + } + if (!trackTempFilterMap) { + continue; + } + + // store filtering information + if (track.isGlobalTrack()) { + trackFilteringTag |= (uint64_t(1) << 0); // BIT0: global track + } + if (track.isGlobalTrackSDD()) { + trackFilteringTag |= (uint64_t(1) << 1); // BIT1: global track SSD + } + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::TrackV0Bits)) { // BIT2-6: V0Bits + trackFilteringTag |= (uint64_t(track.pidbit()) << 2); + for (int iv0 = 0; iv0 < 5; iv0++) { + if (track.pidbit() & (uint8_t(1) << iv0)) { + (reinterpret_cast(fStatsList->At(1)))->Fill(fTrackCuts.size() + static_cast(iv0)); + } + } + } + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::DalitzBits)) { + trackFilteringTag |= (uint64_t(track.dalitzBits()) << 7); // BIT7-14: Dalitz + } + trackFilteringTag |= (uint64_t(trackTempFilterMap) << 15); // BIT15-...: user track filters + + mcflags = 0; + i = 0; // runs over the MC signals + int j = 0; // runs over the track cuts + // check all the specified signals and fill histograms for MC truth matched tracks + for (auto& sig : fMCSignals) { + if (sig.CheckSignal(true, mctrack)) { + mcflags |= (uint16_t(1) << i); + if (fDoDetailedQA) { + j = 0; + for (auto& cut : fTrackCuts) { + if (trackTempFilterMap & (uint8_t(1) << j)) { + fHistMan->FillHistClass(Form("TrackBarrel_%s_%s", cut.GetName(), sig.GetName()), VarManager::fgValues); // fill the reconstructed truth + } + j++; + } + } + } + i++; + } + + // if the MC truth particle corresponding to this reconstructed track is not already written, + // add it to the skimmed stack + if (!(fNewLabels.find(mctrack.globalIndex()) != fNewLabels.end())) { + fNewLabels[mctrack.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mctrack.globalIndex(); + fMCFlags[mctrack.globalIndex()] = mcflags; + fEventIdx[mctrack.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; + fCounters[0]++; + } + + trackBasic(event.lastIndex(), trackFilteringTag, track.pt(), track.eta(), track.phi(), track.sign(), isAmbiguous); + trackBarrel(track.x(), track.alpha(), track.y(), track.z(), track.snp(), track.tgl(), track.signed1Pt(), + track.tpcInnerParam(), track.flags(), track.itsClusterMap(), track.itsChi2NCl(), + track.tpcNClsFindable(), track.tpcNClsFindableMinusFound(), track.tpcNClsFindableMinusCrossedRows(), + track.tpcNClsShared(), track.tpcChi2NCl(), + track.trdChi2(), track.trdPattern(), track.tofChi2(), + track.length(), track.dcaXY(), track.dcaZ(), + track.trackTime(), track.trackTimeRes(), track.tofExpMom(), + track.detectorMap()); + trackBarrelPID(track.tpcSignal(), + track.tpcNSigmaEl(), track.tpcNSigmaMu(), + track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr(), + track.beta(), + track.tofNSigmaEl(), track.tofNSigmaMu(), + track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr(), + track.trdSignal()); + trackBarrelLabels(fNewLabels.find(mctrack.index())->second, track.mcMask(), mcflags); + if constexpr (static_cast(TTrackFillMap & VarManager::ObjTypes::TrackCov)) { + trackBarrelCov(track.cYY(), track.cZY(), track.cZZ(), track.cSnpY(), track.cSnpZ(), + track.cSnpSnp(), track.cTglY(), track.cTglZ(), track.cTglSnp(), track.cTglTgl(), + track.c1PtY(), track.c1PtZ(), track.c1PtSnp(), track.c1PtTgl(), track.c1Pt21Pt2()); + } + } // end loop over reconstructed tracks + } // end if constexpr (static_cast(TTrackFillMap)) + + if constexpr (static_cast(TMuonFillMap)) { + // build the muon tables + muonBasic.reserve(tracksMuon.size()); + muonExtra.reserve(tracksMuon.size()); + if constexpr (static_cast(TMuonFillMap & VarManager::ObjTypes::MuonCov)) { + muonCov.reserve(tracksMuon.size()); + } + muonLabels.reserve(tracksMuon.size()); // TODO: enable this once we have fwdtrack labels + + auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + // loop over muons + + // first we need to get the correct indices + int nDel = 0; + int idxPrev = -1; + std::map newEntryNb; + std::map newMatchIndex; + + for (auto& muonId : fwdtrackIdsThisCollision) { + auto muon = muonId.template fwdtrack_as(); + trackFilteringTag = uint64_t(0); + trackTempFilterMap = uint8_t(0); + + if (!muon.has_mcParticle()) { + continue; + } + + VarManager::FillTrack(muon); + if (fPropMuon) { + VarManager::FillPropagateMuon(muon, collision); + } + + if (muon.index() > idxPrev + 1) { // checks if some muons are filtered even before the skimming function + nDel += muon.index() - (idxPrev + 1); + } + idxPrev = muon.index(); + + // check the cuts and filters + int i = 0; + for (auto& cut : fMuonCuts) { + if (cut.IsSelected(VarManager::fgValues)) { + trackTempFilterMap |= (uint8_t(1) << i); + if (fConfigQA) { + fHistMan->FillHistClass(Form("Muons_%s", cut.GetName()), VarManager::fgValues); + } + (reinterpret_cast(fStatsList->At(2)))->Fill(static_cast(i)); + } + i++; + } + if (!trackTempFilterMap) { // does not pass the cuts + nDel++; + } else { // it passes the cuts and will be saved in the tables + newEntryNb[muon.index()] = muon.index() - nDel; + } + } + + // now let's save the muons with the correct indices and matches + for (auto& muonId : fwdtrackIdsThisCollision) { + auto muon = muonId.template fwdtrack_as(); + if constexpr ((TMuonFillMap & VarManager::ObjTypes::AmbiMuon) > 0) { + if (fIsAmbiguous) { + isAmbiguous = (muon.compatibleCollIds().size() != 1); + } + } + + trackFilteringTag = uint64_t(0); + trackTempFilterMap = uint8_t(0); + + if (!muon.has_mcParticle()) { + continue; + } + auto mctrack = muon.template mcParticle_as(); + VarManager::FillTrack(muon); + // recalculte pDca for global muon tracks + if (static_cast(muon.trackType()) < 2) { + auto const& matchMCH = tracksMuon.rawIteratorAt(static_cast(muon.matchMCHTrackId())); + VarManager::FillMuonPDca(matchMCH, collision); + } else if (static_cast(muon.trackType()) > 2) { + VarManager::FillMuonPDca(muon, collision); + } + if (fPropMuon) { + VarManager::FillPropagateMuon(muon, collision); + } + VarManager::FillTrackMC(mcTracks, mctrack); + + if (fDoDetailedQA) { + fHistMan->FillHistClass("Muons_BeforeCuts", VarManager::fgValues); + if (fIsAmbiguous && isAmbiguous == 1) { + fHistMan->FillHistClass("Ambiguous_Muons_BeforeCuts", VarManager::fgValues); + } + } + // apply the muon selection cuts and fill the stats histogram + int i = 0; + for (auto& cut : fMuonCuts) { + if (cut.IsSelected(VarManager::fgValues)) { + trackTempFilterMap |= (uint8_t(1) << i); + fHistMan->FillHistClass(Form("Muons_%s", cut.GetName()), VarManager::fgValues); + if (fIsAmbiguous && isAmbiguous == 1) { + fHistMan->FillHistClass(Form("Ambiguous_Muons_%s", cut.GetName()), VarManager::fgValues); + } + (reinterpret_cast(fStatsList->At(2)))->Fill(static_cast(i)); + } + i++; + } + if (!trackTempFilterMap) { + continue; + } + // store the cut decisions + trackFilteringTag |= uint64_t(trackTempFilterMap); // BIT0-7: user selection cuts + + mcflags = 0; + i = 0; // runs over the MC signals + int j = 0; // runs over the track cuts + // check all the specified signals and fill histograms for MC truth matched tracks + for (auto& sig : fMCSignals) { + if (sig.CheckSignal(true, mctrack)) { + mcflags |= (uint16_t(1) << i); + if (fDoDetailedQA) { + for (auto& cut : fMuonCuts) { + if (trackTempFilterMap & (uint8_t(1) << j)) { + fHistMan->FillHistClass(Form("Muons_%s_%s", cut.GetName(), sig.GetName()), VarManager::fgValues); // fill the reconstructed truth + } + j++; + } + } + } + i++; + } + + // if the MC truth particle corresponding to this reconstructed track is not already written, + // add it to the skimmed stack + if (!(fNewLabels.find(mctrack.globalIndex()) != fNewLabels.end())) { + fNewLabels[mctrack.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mctrack.globalIndex(); + fMCFlags[mctrack.globalIndex()] = mcflags; + fEventIdx[mctrack.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; + fCounters[0]++; + } + + // update the matching MCH/MFT index + + if (static_cast(muon.trackType()) == 0 || static_cast(muon.trackType()) == 2) { // MCH-MFT or GLB track + int matchIdx = muon.matchMCHTrackId() - muon.offsets(); + if (newEntryNb.count(matchIdx) > 0) { // if the key exists which means the match will not get deleted + newMatchIndex[muon.index()] = newEntryNb[matchIdx]; // update the match for this muon to the updated entry of the match + newMatchIndex[muon.index()] += muonBasic.lastIndex() + 1 - newEntryNb[muon.index()]; // adding the offset of muons, muonBasic.lastIndex() start at -1 + if (static_cast(muon.trackType()) == 0) { // for now only do this to global tracks + newMatchIndex[matchIdx] = newEntryNb[muon.index()]; // add the updated index of this muon as a match to mch track + newMatchIndex[matchIdx] += muonBasic.lastIndex() + 1 - newEntryNb[muon.index()]; // adding the offset, muonBasic.lastIndex() start at -1 + } + } else { + newMatchIndex[muon.index()] = -1; + } + } else if (static_cast(muon.trackType() == 4)) { // an MCH track + // in this case the matches should be filled from the other types but we need to check + if (newMatchIndex.count(muon.index()) == 0) { + newMatchIndex[muon.index()] = -1; + } + } + + muonBasic(event.lastIndex(), newMatchIndex.find(muon.index())->second, -1, trackFilteringTag, VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], muon.sign(), isAmbiguous); + if constexpr (static_cast(TMuonFillMap & VarManager::ObjTypes::MuonCov)) { + if (fPropMuon) { + muonExtra(muon.nClusters(), VarManager::fgValues[VarManager::kMuonPDca], VarManager::fgValues[VarManager::kMuonRAtAbsorberEnd], + muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(), + muon.matchScoreMCHMFT(), muon.mchBitMap(), muon.midBitMap(), + muon.midBoards(), muon.trackType(), VarManager::fgValues[VarManager::kMuonDCAx], VarManager::fgValues[VarManager::kMuonDCAy], + muon.trackTime(), muon.trackTimeRes()); + } + muonCov(VarManager::fgValues[VarManager::kX], VarManager::fgValues[VarManager::kY], VarManager::fgValues[VarManager::kZ], VarManager::fgValues[VarManager::kPhi], VarManager::fgValues[VarManager::kTgl], muon.sign() / VarManager::fgValues[VarManager::kPt], + VarManager::fgValues[VarManager::kMuonCXX], VarManager::fgValues[VarManager::kMuonCXY], VarManager::fgValues[VarManager::kMuonCYY], VarManager::fgValues[VarManager::kMuonCPhiX], VarManager::fgValues[VarManager::kMuonCPhiY], VarManager::fgValues[VarManager::kMuonCPhiPhi], + VarManager::fgValues[VarManager::kMuonCTglX], VarManager::fgValues[VarManager::kMuonCTglY], VarManager::fgValues[VarManager::kMuonCTglPhi], VarManager::fgValues[VarManager::kMuonCTglTgl], VarManager::fgValues[VarManager::kMuonC1Pt2X], VarManager::fgValues[VarManager::kMuonC1Pt2Y], + VarManager::fgValues[VarManager::kMuonC1Pt2Phi], VarManager::fgValues[VarManager::kMuonC1Pt2Tgl], VarManager::fgValues[VarManager::kMuonC1Pt21Pt2]); + } else { + muonExtra(muon.nClusters(), muon.pDca(), muon.rAtAbsorberEnd(), + muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(), + muon.matchScoreMCHMFT(), muon.mchBitMap(), muon.midBitMap(), + muon.midBoards(), muon.trackType(), muon.fwdDcaX(), muon.fwdDcaY(), + muon.trackTime(), muon.trackTimeRes()); + } + muonLabels(fNewLabels.find(mctrack.index())->second, muon.mcMask(), mcflags); + } + } // end if constexpr (static_cast(TMuonFillMap)) + } // end loop over collisions + + // Loop over the label map, create the mother/daughter relationships if these exist and write the skimmed MC stack + for (const auto& [newLabel, oldLabel] : fNewLabelsReversed) { + auto mctrack = mcTracks.iteratorAt(oldLabel); + uint16_t mcflags = fMCFlags.find(oldLabel)->second; + + std::vector mothers; + if (mctrack.has_mothers()) { + for (auto& m : mctrack.mothersIds()) { + if (m < mcTracks.size()) { // protect against bad mother indices + if (fNewLabels.find(m) != fNewLabels.end()) { + mothers.push_back(fNewLabels.find(m)->second); + } + } else { + cout << "Mother label (" << m << ") exceeds the McParticles size (" << mcTracks.size() << ")" << endl; + cout << " Check the MC generator" << endl; + } + } + } + + // TODO: Check that the daughter slice in the skimmed table works as expected + // Note that not all daughters from the original table are preserved in the skimmed MC stack + std::vector daughters; + if (mctrack.has_daughters()) { + for (int d = mctrack.daughtersIds()[0]; d <= mctrack.daughtersIds()[1]; ++d) { + // TODO: remove this check as soon as issues with MC production are fixed + if (d < mcTracks.size()) { // protect against bad daughter indices + if (fNewLabels.find(d) != fNewLabels.end()) { + daughters.push_back(fNewLabels.find(d)->second); + } + } else { + cout << "Daughter label (" << d << ") exceeds the McParticles size (" << mcTracks.size() << ")" << endl; + cout << " Check the MC generator" << endl; + } + } + } + int daughterRange[2] = {-1, -1}; + if (daughters.size() > 0) { + daughterRange[0] = daughters[0]; + daughterRange[1] = daughters[daughters.size() - 1]; + } + + trackMC(fEventIdx.find(oldLabel)->second, mctrack.pdgCode(), mctrack.statusCode(), mctrack.flags(), + mothers, daughterRange, + mctrack.weight(), mctrack.pt(), mctrack.eta(), mctrack.phi(), mctrack.e(), + mctrack.vx(), mctrack.vy(), mctrack.vz(), mctrack.vt(), mcflags); + for (unsigned int isig = 0; isig < fMCSignals.size(); isig++) { + if (mcflags & (uint16_t(1) << isig)) { + (reinterpret_cast(fStatsList->At(3)))->Fill(static_cast(isig)); + } + } + if (mcflags == 0) { + (reinterpret_cast(fStatsList->At(3)))->Fill(static_cast(fMCSignals.size())); + } + } // end loop over labels + + fCounters[0] = 0; + fCounters[1] = 0; + fNewLabels.clear(); + fNewLabelsReversed.clear(); + fMCFlags.clear(); + fEventIdx.clear(); + fEventLabels.clear(); + } + void DefineHistograms(TString histClasses) { std::unique_ptr objArray(histClasses.Tokenize(";")); @@ -1040,6 +1590,30 @@ struct TableMakerMC { { fullSkimming(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, nullptr); } + // Produce muon tables only based on track-collision association tables -------------------------------------------------------------------------------------- + void processAssociatedMuonOnly(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + { + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + } + + void processAssociatedMuonOnlyWithCov(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + { + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + } + + void processAssociatedMuonOnlyWithCovAndCent(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + { + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + } + + void processAssociatedMuonOnlyWithCovAndMults(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + { + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + } // Produce muon tables only for ambiguous tracks studies -------------------------------------------------------------------------------------- void processAmbiguousMuonOnly(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, soa::Filtered const& tracksMuon, @@ -1085,6 +1659,10 @@ struct TableMakerMC { PROCESS_SWITCH(TableMakerMC, processMuonOnlyWithCov, "Produce muon skims, with muon covariance matrix", false); PROCESS_SWITCH(TableMakerMC, processMuonOnlyWithCent, "Produce muon skims, w/ centrality", false); PROCESS_SWITCH(TableMakerMC, processOnlyBCs, "Analyze the BCs to store sampled lumi", false); + PROCESS_SWITCH(TableMakerMC, processAssociatedMuonOnly, "Produce muon skims using track-collision association tables", false); + PROCESS_SWITCH(TableMakerMC, processAssociatedMuonOnlyWithCov, "Produce muon skims using track-collision association tables, with muon covariance matrix", false); + PROCESS_SWITCH(TableMakerMC, processAssociatedMuonOnlyWithCovAndCent, "Produce muon skims using track-collision association tables, w/ centrality", false); + PROCESS_SWITCH(TableMakerMC, processAssociatedMuonOnlyWithCovAndMults, "Produce muon skims using track-collision association tables, w/ mults", false); PROCESS_SWITCH(TableMakerMC, processAmbiguousMuonOnly, "Build muon-only DQ skimmed data model with QA plots for ambiguous muons", false); PROCESS_SWITCH(TableMakerMC, processAmbiguousMuonOnlyWithCov, "Build muon-only with cov DQ skimmed data model with QA plots for ambiguous muons", false); PROCESS_SWITCH(TableMakerMC, processAmbiguousBarrelOnly, "Build barrel-only DQ skimmed data model with QA plots for ambiguous tracks", false); From fd3fe9e15521390d11410b60c3a8bf88956211ec Mon Sep 17 00:00:00 2001 From: motomioya Date: Sun, 19 May 2024 00:41:16 +0900 Subject: [PATCH 8/9] Modified to follow conventions --- PWGDQ/TableProducer/tableMakerMC.cxx | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/PWGDQ/TableProducer/tableMakerMC.cxx b/PWGDQ/TableProducer/tableMakerMC.cxx index 4a900066b3f..7ad83d257ac 100644 --- a/PWGDQ/TableProducer/tableMakerMC.cxx +++ b/PWGDQ/TableProducer/tableMakerMC.cxx @@ -232,7 +232,7 @@ struct TableMakerMC { context.mOptions.get("processMuonOnlyWithCent") || context.mOptions.get("processMuonOnlyWithCov") || context.mOptions.get("processAmbiguousMuonOnlyWithCov") || context.mOptions.get("processAmbiguousMuonOnly") || context.mOptions.get("processAssociatedMuonOnly") || context.mOptions.get("processAssociatedMuonOnlyWithCov") || - context.mOptions.get("processAssociatedMuonOnlyWithCovAndCent") || context.mOptions.get("processAssociatedMuonOnlyWithCovAndMults"); + context.mOptions.get("processAssociatedMuonOnlyWithCovAndCent") || context.mOptions.get("processAssociatedMuonOnlyWithCovAndMults")); // TODO: switch on/off histogram classes depending on which process function we run if (enableBarrelHistos) { if (fDoDetailedQA) { @@ -1592,25 +1592,29 @@ struct TableMakerMC { } // Produce muon tables only based on track-collision association tables -------------------------------------------------------------------------------------- void processAssociatedMuonOnly(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + soa::Filtered const& tracksMuon, + aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } void processAssociatedMuonOnlyWithCov(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + soa::Filtered const& tracksMuon, + aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } void processAssociatedMuonOnlyWithCovAndCent(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + soa::Filtered const& tracksMuon, + aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } void processAssociatedMuonOnlyWithCovAndMults(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) + soa::Filtered const& tracksMuon, + aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } From 4f76a96521c6f981463f9a5dc94c49ad1112a748 Mon Sep 17 00:00:00 2001 From: motomioya Date: Wed, 22 May 2024 13:22:22 +0900 Subject: [PATCH 9/9] Add DCA recalculation for muons --- PWGDQ/Core/VarManager.h | 12 ++++++++++++ PWGDQ/TableProducer/tableMaker.cxx | 14 +++++++------- PWGDQ/TableProducer/tableMakerMC.cxx | 22 +++++++++++----------- 3 files changed, 30 insertions(+), 18 deletions(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index c313adaf552..f5779cb204d 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -2049,6 +2049,18 @@ void VarManager::FillTrackCollision(T const& track, C const& collision, float* v } } } + if constexpr ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0) { + + o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(track, collision); + o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(track, collision, kToDCA); + + float dcaX = (propmuonAtDCA.getX() - collision.posX()); + float dcaY = (propmuonAtDCA.getY() - collision.posY()); + float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + values[kMuonPDca] = track.p() * dcaXY; + values[kMuonDCAx] = dcaX; + values[kMuonDCAy] = dcaY; + } } template diff --git a/PWGDQ/TableProducer/tableMaker.cxx b/PWGDQ/TableProducer/tableMaker.cxx index fdf51557de7..9f96e5c5f8b 100644 --- a/PWGDQ/TableProducer/tableMaker.cxx +++ b/PWGDQ/TableProducer/tableMaker.cxx @@ -668,9 +668,9 @@ struct TableMaker { // recalculte pDca for global muon tracks if (static_cast(muon.trackType()) < 2) { auto const& matchMCH = muon.template matchMCHTrack_as(); - VarManager::FillMuonPDca(matchMCH, collision); + VarManager::FillTrackCollision(matchMCH, collision); } else if (static_cast(muon.trackType()) > 2) { - VarManager::FillMuonPDca(muon, collision); + VarManager::FillTrackCollision(muon, collision); } if (fPropMuon) { @@ -1039,9 +1039,9 @@ struct TableMaker { // recalculte pDca for global muon tracks if (static_cast(muon.trackType()) < 2) { auto const& matchMCH = tracksMuon.rawIteratorAt(static_cast(muon.matchMCHTrackId())); - VarManager::FillMuonPDca(matchMCH, collision); + VarManager::FillTrackCollision(matchMCH, collision); } else if (static_cast(muon.trackType()) > 2) { - VarManager::FillMuonPDca(muon, collision); + VarManager::FillTrackCollision(muon, collision); } if (fPropMuon) { @@ -1545,7 +1545,7 @@ struct TableMaker { { for (auto& collision : collisions) { auto muonIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); - fullSkimmingIndices(collision, bcs, nullptr, tracksMuon, nullptr, muonIdsThisCollision); + fullSkimmingIndices(collision, bcs, nullptr, tracksMuon, nullptr, muonIdsThisCollision); } } @@ -1554,7 +1554,7 @@ struct TableMaker { { for (auto& collision : collisions) { auto muonIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); - fullSkimmingIndices(collision, bcs, nullptr, tracksMuon, nullptr, muonIdsThisCollision); + fullSkimmingIndices(collision, bcs, nullptr, tracksMuon, nullptr, muonIdsThisCollision); } } @@ -1563,7 +1563,7 @@ struct TableMaker { { for (auto& collision : collisions) { auto muonIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); - fullSkimmingIndices(collision, bcs, nullptr, tracksMuon, nullptr, muonIdsThisCollision); + fullSkimmingIndices(collision, bcs, nullptr, tracksMuon, nullptr, muonIdsThisCollision); } } diff --git a/PWGDQ/TableProducer/tableMakerMC.cxx b/PWGDQ/TableProducer/tableMakerMC.cxx index 7ad83d257ac..df6e47f8818 100644 --- a/PWGDQ/TableProducer/tableMakerMC.cxx +++ b/PWGDQ/TableProducer/tableMakerMC.cxx @@ -1254,9 +1254,9 @@ struct TableMakerMC { // recalculte pDca for global muon tracks if (static_cast(muon.trackType()) < 2) { auto const& matchMCH = tracksMuon.rawIteratorAt(static_cast(muon.matchMCHTrackId())); - VarManager::FillMuonPDca(matchMCH, collision); + VarManager::FillTrackCollision(matchMCH, collision); } else if (static_cast(muon.trackType()) > 2) { - VarManager::FillMuonPDca(muon, collision); + VarManager::FillTrackCollision(muon, collision); } if (fPropMuon) { VarManager::FillPropagateMuon(muon, collision); @@ -1595,28 +1595,28 @@ struct TableMakerMC { soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { - fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } void processAssociatedMuonOnlyWithCov(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { - fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } - void processAssociatedMuonOnlyWithCovAndCent(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, + void processAssociatedMuonOnlyWithCovAndCent(MyEventsWithCent const& collisions, aod::BCsWithTimestamps const& bcs, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { - fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } - void processAssociatedMuonOnlyWithCovAndMults(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, - soa::Filtered const& tracksMuon, + void processAssociatedMuonOnlyWithCovAndMults(MyEventsWithCentAndMults const& collisions, aod::BCsWithTimestamps const& bcs, + soa::Filtered const& tracksMuon, aod::McCollisions const& mcEvents, aod::McParticles_001 const& mcTracks, aod::FwdTrackAssoc const& fwdtrackIndices) { - fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); + fullSkimmingIndices(collisions, bcs, nullptr, tracksMuon, mcEvents, mcTracks, nullptr, fwdtrackIndices); } // Produce muon tables only for ambiguous tracks studies -------------------------------------------------------------------------------------- void processAmbiguousMuonOnly(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs,