diff --git a/PWGEM/Dilepton/DataModel/lmeeMLTables.h b/PWGEM/Dilepton/DataModel/lmeeMLTables.h index 3036164b9ad..4bb0032c934 100644 --- a/PWGEM/Dilepton/DataModel/lmeeMLTables.h +++ b/PWGEM/Dilepton/DataModel/lmeeMLTables.h @@ -153,19 +153,21 @@ DECLARE_SOA_COLUMN(YMFTatMP, yMFTatMP, float); //! y of MFT track in MFT-M DECLARE_SOA_COLUMN(XErrMFTatMP, xErrMFTatMP, float); //! x error of MFT track in MFT-MCH-MID track at matching plane DECLARE_SOA_COLUMN(YErrMFTatMP, yErrMFTatMP, float); //! y error of MFT track in MFT-MCH-MID track at matching plane -DECLARE_SOA_COLUMN(Chi2MFT, chi2MFT, float); //! chi2/ndf of MFT track -DECLARE_SOA_COLUMN(Chi2MCHMID, chi2MCHMID, float); //! chi2/ndf of MCH-MID track -DECLARE_SOA_COLUMN(Chi2MFTMCHMID, chi2MFTMCHMID, float); //! chi2/ndf of MFT-MCH-MID track -DECLARE_SOA_COLUMN(NClustersMFT, nClustersMFT, uint8_t); //! -DECLARE_SOA_COLUMN(IsPrimary, isPrimary, bool); //! -DECLARE_SOA_COLUMN(IsCorrectMatch, isCorrectMatch, bool); //! +DECLARE_SOA_COLUMN(Chi2MFT, chi2MFT, float); //! chi2/ndf of MFT track +DECLARE_SOA_COLUMN(Chi2MCHMID, chi2MCHMID, float); //! chi2/ndf of MCH-MID track +DECLARE_SOA_COLUMN(Chi2MFTMCHMID, chi2MFTMCHMID, float); //! chi2/ndf of MFT-MCH-MID track +DECLARE_SOA_COLUMN(NClustersMFT, nClustersMFT, uint8_t); //! +DECLARE_SOA_COLUMN(IsPrimaryMFT, isPrimaryMFT, bool); //! +DECLARE_SOA_COLUMN(IsPrimaryMCHMID, isPrimaryMCHMID, bool); //! +DECLARE_SOA_COLUMN(IsCorrectMatch, isCorrectMatch, bool); //! +DECLARE_SOA_COLUMN(PdgCodeMFT, pdgCodeMFT, int); //! +DECLARE_SOA_COLUMN(PdgCodeMCHMID, pdgCodeMCHMID, int); //! DECLARE_SOA_COLUMN(MultMFT, multMFT, uint16_t); //! number of MFTsa tracks per collision } // namespace emmlfwdtrack DECLARE_SOA_TABLE(EMFwdTracksForML, "AOD", "EMFWDTRKML", //! o2::soa::Index<>, collision::PosZ, /*collision::NumContrib,*/ mult::MultFT0C, /*evsel::NumTracksInTimeRange,*/ evsel::SumAmpFT0CInTimeRange, emmltrack::HadronicRate, emmlfwdtrack::MultMFT, - // fwdtrack::TrackType, emmlfwdtrack::Signed1PtMFTatMP, emmlfwdtrack::TglMFTatMP, emmlfwdtrack::PhiMFTatMP, emmlfwdtrack::XMFTatMP, emmlfwdtrack::YMFTatMP, @@ -173,19 +175,18 @@ DECLARE_SOA_TABLE(EMFwdTracksForML, "AOD", "EMFWDTRKML", //! emmlfwdtrack::Signed1PtMCHMIDatMP, emmlfwdtrack::TglMCHMIDatMP, emmlfwdtrack::PhiMCHMIDatMP, emmlfwdtrack::XMCHMIDatMP, emmlfwdtrack::YMCHMIDatMP, - // fwdtrack::NClusters, fwdtrack::PDca, fwdtrack::RAtAbsorberEnd, - // fwdtrack::Chi2MatchMCHMID, fwdtrack::Chi2MatchMCHMFT, - // fwdtrack::MFTClusterSizesAndTrackFlags, emmlfwdtrack::Chi2MFTMCHMID, emmlfwdtrack::Chi2MCHMID, emmlfwdtrack::Chi2MFT, emmlfwdtrack::NClustersMFT, - mcparticle::PdgCode, emmlfwdtrack::IsPrimary, emmlfwdtrack::IsCorrectMatch, mcparticle::Pt, mcparticle::Eta, mcparticle::Phi); + emmlfwdtrack::PdgCodeMFT, emmlfwdtrack::IsPrimaryMFT, + emmlfwdtrack::PdgCodeMCHMID, emmlfwdtrack::IsPrimaryMCHMID, + emmlfwdtrack::IsCorrectMatch); // iterators using EMFwdTrackForML = EMFwdTracksForML::iterator; DECLARE_SOA_TABLE(EMFwdTrackErrsForML, "AOD", "EMFWDTRKERRML", //! Joinable with EMFwdTracksForML - emmlfwdtrack::Signed1PtErrMFTatMP, emmlfwdtrack::TglErrMFTatMP, emmlfwdtrack::PhiErrMFTatMP, + /*emmlfwdtrack::Signed1PtErrMFTatMP,*/ emmlfwdtrack::TglErrMFTatMP, emmlfwdtrack::PhiErrMFTatMP, emmlfwdtrack::XErrMFTatMP, emmlfwdtrack::YErrMFTatMP, - emmlfwdtrack::Signed1PtErrMCHMIDatMP, emmlfwdtrack::TglErrMCHMIDatMP, emmlfwdtrack::PhiErrMCHMIDatMP, + /*emmlfwdtrack::Signed1PtErrMCHMIDatMP,*/ emmlfwdtrack::TglErrMCHMIDatMP, emmlfwdtrack::PhiErrMCHMIDatMP, emmlfwdtrack::XErrMCHMIDatMP, emmlfwdtrack::YErrMCHMIDatMP); // iterators diff --git a/PWGEM/Dilepton/TableProducer/treeCreatorMuonML.cxx b/PWGEM/Dilepton/TableProducer/treeCreatorMuonML.cxx index b4c099ea89e..ece00db378a 100644 --- a/PWGEM/Dilepton/TableProducer/treeCreatorMuonML.cxx +++ b/PWGEM/Dilepton/TableProducer/treeCreatorMuonML.cxx @@ -119,10 +119,16 @@ struct TreeCreatorMuonML { struct : ConfigurableGroup { std::string prefix = "MFTCutGroup"; Configurable minPt{"minPt", 0.04f, "min. pT for MFTsa to reject crazy tracks"}; - Configurable minEta{"minEta", -4.1f, "min. eta acceptance for MFTsa to reject crazy tracks"}; - Configurable maxEta{"maxEta", -2.0f, "max. eta acceptance for MFTsa to reject crazy tracks"}; + Configurable minEta{"minEta", -4.6f, "min. eta acceptance for MFTsa to reject crazy tracks"}; + Configurable maxEta{"maxEta", -1.5f, "max. eta acceptance for MFTsa to reject crazy tracks"}; } MFTCutGroup; + struct : ConfigurableGroup { + std::string prefix = "MCHMIDCutGroup"; + Configurable minEta{"minEta", -4.6f, "min. eta acceptance for MCHMID to reject crazy tracks"}; + Configurable maxEta{"maxEta", -1.5f, "max. eta acceptance for MCHMID to reject crazy tracks"}; + } MCHMIDCutGroup; + o2::ccdb::CcdbApi ccdbApi; Service ccdb; std::mt19937 engine; @@ -215,44 +221,16 @@ struct TreeCreatorMuonML { hMuonType->GetXaxis()->SetBinLabel(4, "MCH-MID"); hMuonType->GetXaxis()->SetBinLabel(5, "MCH standalone"); - fRegistry.add("MFTMCHMID/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{100, 0.0f, 10}}, false); - fRegistry.add("MFTMCHMID/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {60, -5.f, -2.f}}, false); - fRegistry.add("MFTMCHMID/hEtaPhi_MatchedMCHMID", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {60, -5.f, -2.f}}, false); - fRegistry.add("MFTMCHMID/hDeltaPt_Pt", "#Deltap_{T}/p_{T} vs. p_{T};p_{T}^{gl} (GeV/c);(p_{T}^{sa} - p_{T}^{gl})/p_{T}^{gl}", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); - fRegistry.add("MFTMCHMID/hDeltaEta_Pt", "#Delta#eta vs. p_{T};p_{T}^{gl} (GeV/c);#Delta#eta", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); - fRegistry.add("MFTMCHMID/hDeltaPhi_Pt", "#Delta#varphi vs. p_{T};p_{T}^{gl} (GeV/c);#Delta#varphi (rad.)", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); - fRegistry.add("MFTMCHMID/hDeltaEtaAtMP_Pt", "#Delta#eta vs. p_{T} at MP;p_{T}^{gl} (GeV/c);#Delta#eta", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); - fRegistry.add("MFTMCHMID/hDeltaPhiAtMP_Pt", "#Delta#varphi vs. p_{T} at MP;p_{T}^{gl} (GeV/c);#Delta#varphi (rad.)", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); - fRegistry.add("MFTMCHMID/hSign", "sign;sign", kTH1F, {{3, -1.5, +1.5}}, false); - fRegistry.add("MFTMCHMID/hNclusters", "Nclusters;Nclusters", kTH1F, {{21, -0.5f, 20.5}}, false); - fRegistry.add("MFTMCHMID/hNclustersMFT", "NclustersMFT;Nclusters MFT", kTH1F, {{11, -0.5f, 10.5}}, false); - fRegistry.add("MFTMCHMID/hRatAbsorberEnd", "R at absorber end;R at absorber end (cm)", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("MFTMCHMID/hPDCA_Rabs", "pDCA vs. Rabs;R at absorber end (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 1000}}, false); - fRegistry.add("MFTMCHMID/hChi2", "chi2;chi2/ndf", kTH1F, {{100, 0.0f, 10}}, false); - fRegistry.add("MFTMCHMID/hChi2MFT", "chi2 MFT;chi2 MFT/ndf", kTH1F, {{100, 0.0f, 10}}, false); - fRegistry.add("MFTMCHMID/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("MFTMCHMID/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("MFTMCHMID/hDCAxy2D", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -1, 1}, {200, -1, +1}}, false); - fRegistry.add("MFTMCHMID/hDCAxy2DinSigma", "DCA x vs. y in sigma;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10, 10}, {200, -10, +10}}, false); - fRegistry.add("MFTMCHMID/hDCAxy", "DCAxy;DCA_{xy} (cm);", kTH1F, {{100, 0, 1}}, false); - fRegistry.add("MFTMCHMID/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{100, 0, 1}, {200, -0.1, 0.1}}, false); - fRegistry.add("MFTMCHMID/hDCAxyinSigma", "DCAxy in sigma;DCA_{xy} (#sigma);", kTH1F, {{100, 0, 10}}, false); - fRegistry.add("MFTMCHMID/hDCAx_PosZ", "DCAx vs. posZ;Z_{vtx} (cm);DCA_{x} (cm)", kTH2F, {{200, -10, +10}, {400, -0.2, +0.2}}, false); - fRegistry.add("MFTMCHMID/hDCAy_PosZ", "DCAy vs. posZ;Z_{vtx} (cm);DCA_{y} (cm)", kTH2F, {{200, -10, +10}, {400, -0.2, +0.2}}, false); - fRegistry.add("MFTMCHMID/hDCAx_Phi", "DCAx vs. #varphi;#varphi (rad.);DCA_{x} (cm)", kTH2F, {{90, 0, 2 * M_PI}, {400, -0.2, +0.2}}, false); - fRegistry.add("MFTMCHMID/hDCAy_Phi", "DCAy vs. #varphi;#varphi (rad.);DCA_{y} (cm)", kTH2F, {{90, 0, 2 * M_PI}, {400, -0.2, +0.2}}, false); - fRegistry.add("MFTMCHMID/hNmu", "#mu multiplicity;N_{#mu} per collision", kTH1F, {{21, -0.5, 20.5}}, false); - fRegistry.add("MFTMCHMID/hdR_Chi2MatchMCHMFT", "dr vs. matching chi2 MCH-MFT;chi2 match MCH-MFT;#DeltaR", kTH2F, {{200, 0, 50}, {200, 0, 0.5}}, false); - fRegistry.add("MFTMCHMID/hdR_Chi2", "dr vs. chi2;global chi2/ndf;#DeltaR", kTH2F, {{100, 0, 10}, {200, 0, 0.5}}, false); - fRegistry.add("MFTMCHMID/hChi2_Chi2MatchMCHMFT", "chi2 vs. matching chi2 MCH-MFT;chi2 match MCH-MFT;global chi2/ndf", kTH2F, {{200, 0, 50}, {100, 0, 10}}, false); - - fRegistry.addClone("MFTMCHMID/", "MCHMID/"); - fRegistry.add("MFTMCHMID/hDCAxResolutionvsPt", "DCA_{x} vs. p_{T};p_{T} (GeV/c);DCA_{x} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 500}}, false); - fRegistry.add("MFTMCHMID/hDCAyResolutionvsPt", "DCA_{y} vs. p_{T};p_{T} (GeV/c);DCA_{y} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 500}}, false); - fRegistry.add("MFTMCHMID/hDCAxyResolutionvsPt", "DCA_{xy} vs. p_{T};p_{T} (GeV/c);DCA_{y} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 500}}, false); - fRegistry.add("MCHMID/hDCAxResolutionvsPt", "DCA_{x} vs. p_{T};p_{T} (GeV/c);DCA_{x} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 5e+5}}, false); - fRegistry.add("MCHMID/hDCAyResolutionvsPt", "DCA_{y} vs. p_{T};p_{T} (GeV/c);DCA_{y} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 5e+5}}, false); - fRegistry.add("MCHMID/hDCAxyResolutionvsPt", "DCA_{xy} vs. p_{T};p_{T} (GeV/c);DCA_{y} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 5e+5}}, false); + // information at matching plane + fRegistry.add("MCHMID/correct/hP", "p;p (GeV/c)", kTH1F, {{1000, 0.0f, 100}}, false); + fRegistry.add("MCHMID/correct/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); + fRegistry.add("MCHMID/correct/hPz", "pz;p_{z} (GeV/c)", kTH1F, {{1000, -100, 0}}, false); + fRegistry.add("MCHMID/correct/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {80, -5.f, -1.f}}, false); + fRegistry.add("MCHMID/correct/hPEta", "p vs. #eta;#eta;p (GeV/c)", kTH2F, {{80, -5.f, -1.f}, {1000, 0, 100}}, false); + fRegistry.add("MCHMID/correct/hPtEta", "pT vs. #eta;#eta;p_{T} (GeV/c)", kTH2F, {{80, -5.f, -1.f}, {1000, 0, 10}}, false); + fRegistry.add("MCHMID/correct/hPzEta", "pz vs. #eta;#eta;p_{z} (GeV/c)", kTH2F, {{80, -5.f, -1.f}, {1000, -100, 0}}, false); + fRegistry.addClone("MCHMID/correct/", "MCHMID/wrong/"); + fRegistry.addClone("MCHMID/", "MFT/"); } template @@ -279,180 +257,19 @@ struct TreeCreatorMuonML { return true; } - template - bool fillFwdTrackTable(TCollision const& collision, TFwdTrack const& fwdtrack, TFwdTracks const&, TMFTTracks const&, TMFTTracksCov const& mftCovs, const float hadronicRate, const uint16_t nmft) - { - if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - return false; - } - - auto mchtrack = fwdtrack.template matchMCHTrack_as(); // MCH-MID - auto mfttrack = fwdtrack.template matchMFTTrack_as(); // MFTsa - if (!mfttrack.has_mcParticle() || !mchtrack.has_mcParticle() || !fwdtrack.has_mcParticle()) { - return false; - } - float chi2 = fwdtrack.chi2() / (2.f * (mchtrack.nClusters() + mfttrack.nClusters()) - 5.f); - float chi2mft = mfttrack.chi2() / (2.f * mfttrack.nClusters() - 5.f); - - if (mfttrack.eta() < MFTCutGroup.minEta || MFTCutGroup.maxEta < mfttrack.eta() || mfttrack.pt() < MFTCutGroup.minPt) { - return false; - } - - // auto mcParticle_MFTMCHMID = fwdtrack.template mcParticle_as(); // this is identical to mcParticle_MCHMID - auto mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID - auto mcParticle_MFT = mfttrack.template mcParticle_as(); - - int pdgCode = mcParticle_MCHMID.pdgCode(); - bool isPrimary = mcParticle_MCHMID.isPhysicalPrimary() || mcParticle_MCHMID.producedByGenerator(); - bool isMatched = (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) && (mcParticle_MFT.mcCollisionId() == mcParticle_MCHMID.mcCollisionId()); - - if (!isMatched && dist01(engine) > cfgDownSampling) { - return false; - } - - if (fwdtrack.chi2MatchMCHMID() < 0.f) { // this should never happen. only for protection. - return false; - } - if (fwdtrack.chi2MatchMCHMFT() < 0.f || maxMatchingChi2MCHMFT < fwdtrack.chi2MatchMCHMFT()) { - return false; - } - // if (fwdtrack.chi2() < 0.f || glMuonCutGroup.maxChi2 < chi2) { - // return false; - // } - if (mfttrack.chi2() < 0.f) { // this should never happen. only for protection. - return false; - } - - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); - float pt = propmuonAtPV.getPt(); - float eta = propmuonAtPV.getEta(); - float phi = propmuonAtPV.getPhi(); - // o2::math_utils::bringTo02Pi(phi); - phi = RecoDecay::constrainAngle(phi, 0, 1U); - - // if (eta < glMuonCutGroup.minEta || glMuonCutGroup.maxEta < eta) { - // return false; - // } - - float dcaX = propmuonAtPV.getX() - collision.posX(); - float dcaY = propmuonAtPV.getY() - collision.posY(); - // float dcaZ = propmuonAtPV.getZ() - collision.posZ(); - float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); - // float ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); - float etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); - float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); - // o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - phiMatchedMCHMID = RecoDecay::constrainAngle(phiMatchedMCHMID, 0, 1U); - if (glMuonCutGroup.refitGlobalMuon) { - pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); - } - - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift); - float dcaX_Matched = propmuonAtDCA_Matched.getX() - collision.posX(); - float dcaY_Matched = propmuonAtDCA_Matched.getY() - collision.posY(); - float dcaXY_Matched = std::sqrt(dcaX_Matched * dcaX_Matched + dcaY_Matched * dcaY_Matched); - float pDCA = mchtrack.p() * dcaXY_Matched; - float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack - - float xMFTatMP = 999.f, yMFTatMP = 999.f; - float xMCHMIDatMP = 999.f, yMCHMIDatMP = 999.f; - - float xErrMFTatMP = 999.f, yErrMFTatMP = 999.f; - float xErrMCHMIDatMP = 999.f, yErrMCHMIDatMP = 999.f; - float signed1PtMFTatMP = 999.f, tglMFTatMP = 999.f, phiMFTatMP = 999.f; - float signed1PtMCHMIDatMP = 999.f, tglMCHMIDatMP = 999.f, phiMCHMIDatMP = 999.f; - - float signed1PtErrMFTatMP = 999.f, tglErrMFTatMP = 999.f, phiErrMFTatMP = 999.f; - float signed1PtErrMCHMIDatMP = 999.f, tglErrMCHMIDatMP = 999.f, phiErrMCHMIDatMP = 999.f; - - if constexpr (withMFTCov) { - auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwdShift(mfttrack, mZShift, mfttrackcov); // values at innermost update - mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane - xMFTatMP = mftsaAtMP.getX(); - yMFTatMP = mftsaAtMP.getY(); - xErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2X()); - yErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Y()); - - signed1PtMFTatMP = mftsaAtMP.getInvQPt(); - tglMFTatMP = mftsaAtMP.getTanl(); - phiMFTatMP = RecoDecay::constrainAngle(mftsaAtMP.getPhi(), 0, 1U); - signed1PtErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2InvQPt()); - tglErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Tanl()); - phiErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Phi()); - - auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz, mZShift); // propagated to matching plane - xMCHMIDatMP = muonAtMP.getX(); - yMCHMIDatMP = muonAtMP.getY(); - xErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2X()); - yErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2Y()); - - signed1PtMCHMIDatMP = muonAtMP.getInvQPt(); - tglMCHMIDatMP = muonAtMP.getTanl(); - phiMCHMIDatMP = RecoDecay::constrainAngle(muonAtMP.getPhi(), 0, 1U); - signed1PtErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2InvQPt()); - tglErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2Tanl()); - phiErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2Phi()); - } - - float deta = etaMatchedMCHMID - eta; - float dphi = phiMatchedMCHMID - phi; - o2::math_utils::bringToPMPi(dphi); - - trackTable(collision.posZ(), /*collision.numContrib(),*/ collision.multFT0C(), /*collision.trackOccupancyInTimeRange(),*/ collision.ft0cOccupancyInTimeRange(), hadronicRate, nmft, - // fwdtrack.trackType(), - - signed1PtMFTatMP, tglMFTatMP, phiMFTatMP, - xMFTatMP, yMFTatMP, - - signed1PtMCHMIDatMP, tglMCHMIDatMP, phiMCHMIDatMP, - xMCHMIDatMP, yMCHMIDatMP, - - // fwdtrack.chi2MatchMCHMID(), - fwdtrack.chi2MatchMCHMFT(), - pdgCode, isPrimary, isMatched, - mcParticle_MCHMID.pt(), mcParticle_MCHMID.eta(), mcParticle_MCHMID.phi()); - - trackErrTable(signed1PtErrMFTatMP, tglErrMFTatMP, phiErrMFTatMP, - xErrMFTatMP, yErrMFTatMP, - signed1PtErrMCHMIDatMP, tglErrMCHMIDatMP, phiErrMCHMIDatMP, - xErrMCHMIDatMP, yErrMCHMIDatMP); - - fRegistry.fill(HIST("hMuonType"), fwdtrack.trackType()); - fRegistry.fill(HIST("MFTMCHMID/hPt"), pt); - fRegistry.fill(HIST("MFTMCHMID/hEtaPhi"), phi, eta); - fRegistry.fill(HIST("MFTMCHMID/hEtaPhi_MatchedMCHMID"), phiMatchedMCHMID, etaMatchedMCHMID); - fRegistry.fill(HIST("MFTMCHMID/hDeltaEta_Pt"), pt, deta); - fRegistry.fill(HIST("MFTMCHMID/hDeltaPhi_Pt"), pt, dphi); - fRegistry.fill(HIST("MFTMCHMID/hSign"), fwdtrack.sign()); - fRegistry.fill(HIST("MFTMCHMID/hNclusters"), fwdtrack.nClusters()); - fRegistry.fill(HIST("MFTMCHMID/hNclustersMFT"), mfttrack.nClusters()); - fRegistry.fill(HIST("MFTMCHMID/hPDCA_Rabs"), rAtAbsorberEnd, pDCA); - fRegistry.fill(HIST("MFTMCHMID/hRatAbsorberEnd"), rAtAbsorberEnd); - fRegistry.fill(HIST("MFTMCHMID/hChi2"), chi2); - fRegistry.fill(HIST("MFTMCHMID/hChi2MFT"), chi2mft); - fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); - fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); - fRegistry.fill(HIST("MFTMCHMID/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/hDCAxy"), dcaXY); - fRegistry.fill(HIST("MFTMCHMID/hDCAx_PosZ"), collision.posZ(), dcaX); - fRegistry.fill(HIST("MFTMCHMID/hDCAy_PosZ"), collision.posZ(), dcaY); - fRegistry.fill(HIST("MFTMCHMID/hDCAx_Phi"), phi, dcaX); - fRegistry.fill(HIST("MFTMCHMID/hDCAy_Phi"), phi, dcaY); - return true; - } - SliceCache cache; Preslice perCollision = o2::aod::fwdtrack::collisionId; Preslice perCollision_MFT = o2::aod::fwdtrack::collisionId; + PresliceUnsorted fwdtracksPerMCHTrack = aod::fwdtrack::matchMCHTrackId; + Partition glMuons = o2::aod::fwdtrack::trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack); + Partition saMuons = o2::aod::fwdtrack::trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack); std::unordered_map map_mfttrackcovs; - void processWithMFTCov(MyCollisionsMC const& collisions, aod::BCsWithTimestamps const&, MyFwdTracksMC const& fwdtracks, MyMFTTracksMC const& mfttracks, aod::MFTTracksCov const& mftCovs, aod::McParticles const&, aod::McCollisions const&) + + void processWithMFTCov(MyCollisionsMC const& collisions, aod::BCsWithTimestamps const&, MyFwdTracksMC const&, MyMFTTracksMC const& mfttracks, aod::MFTTracksCov const& mftCovs, aod::McParticles const&, aod::McCollisions const&) { - for (const auto& mfttrackConv : mftCovs) { - map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + for (const auto& mfttrackCov : mftCovs) { + map_mfttrackcovs[mfttrackCov.matchMFTTrackId()] = mfttrackCov.globalIndex(); } for (const auto& collision : collisions) { @@ -474,66 +291,139 @@ struct TreeCreatorMuonML { float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSourceForCptFetcher) * 1.e-3; // kHz auto mfttracks_per_collision = mfttracks.sliceBy(perCollision_MFT, collision.globalIndex()); - uint16_t nmft = mfttracks_per_collision.size(); + uint16_t multMFT = mfttracks_per_collision.size(); fRegistry.fill(HIST("Event/hCorrFT0CvsMFT"), collision.multFT0C(), mfttracks_per_collision.size()); - auto fwdtracks_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); - for (const auto& fwdtrack : fwdtracks_coll) { - if (!fwdtrack.has_mcParticle()) { + auto saMuons_per_coll = saMuons->sliceByCached(o2::aod::fwdtrack::collisionId, collision.globalIndex(), cache); + for (const auto& mchtrack : saMuons_per_coll) { + if (!mchtrack.has_mcParticle()) { continue; } - if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + if (mchtrack.chi2() < 0.f) { + continue; + } + if (mchtrack.chi2MatchMCHMID() < 0.f) { continue; } - fillFwdTrackTable(collision, fwdtrack, fwdtracks, mfttracks, mftCovs, hadronicRate, nmft); - - } // end of fwdtrack loop - } // end of collision loop - - map_mfttrackcovs.clear(); - } - PROCESS_SWITCH(TreeCreatorMuonML, processWithMFTCov, "produce ML input for single track level", true); - - void processWithoutMFTCov(MyCollisionsMC const& collisions, aod::BCsWithTimestamps const&, MyFwdTracksMC const& fwdtracks, MyMFTTracksMC const& mfttracks, aod::McParticles const&, aod::McCollisions const&) - { - for (const auto& collision : collisions) { - auto bc = collision.template foundBC_as(); - initCCDB(bc); - - if (!collision.has_mcCollision()) { - continue; - } - - if (!isSelectedCollision(collision)) { - continue; - } - - if (eventCutGroup.cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { - continue; - } - - float hadronicRate = mRateFetcher.fetch(ccdb.service, bc.timestamp(), bc.runNumber(), irSourceForCptFetcher) * 1.e-3; // kHz - - auto mfttracks_per_collision = mfttracks.sliceBy(perCollision_MFT, collision.globalIndex()); - uint16_t nmft = mfttracks_per_collision.size(); - fRegistry.fill(HIST("Event/hCorrFT0CvsMFT"), collision.multFT0C(), mfttracks_per_collision.size()); - - auto fwdtracks_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); - for (const auto& fwdtrack : fwdtracks_coll) { - if (!fwdtrack.has_mcParticle()) { + if (dist01(engine) > cfgDownSampling) { // random sampling, if necessary continue; } - if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + + auto mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID + bool isPrimary_MCHMID = mcParticle_MCHMID.isPhysicalPrimary() || mcParticle_MCHMID.producedByGenerator(); + + auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz, mZShift); // propagated to matching plane + float signed1PtMCHMIDatMP = muonAtMP.getInvQPt(); + float tglMCHMIDatMP = muonAtMP.getTanl(); + float phiMCHMIDatMP = RecoDecay::constrainAngle(muonAtMP.getPhi(), 0, 1U); + float xMCHMIDatMP = muonAtMP.getX(); + float yMCHMIDatMP = muonAtMP.getY(); + + float tglErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2Tanl()); + float phiErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2Phi()); + float xErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2X()); + float yErrMCHMIDatMP = std::sqrt(muonAtMP.getSigma2Y()); + if (muonAtMP.getEta() < MCHMIDCutGroup.minEta || MCHMIDCutGroup.maxEta < muonAtMP.getEta()) { continue; } - fillFwdTrackTable(collision, fwdtrack, fwdtracks, mfttracks, nullptr, hadronicRate, nmft); + // auto glMuons_per_MCHMID = fwdtracks.sliceBy(fwdtracksPerMCHTrack, mchtrack.globalIndex()); + auto glMuons_per_MCHMID = glMuons->sliceByCachedUnsorted(o2::aod::fwdtrack::matchMCHTrackId, mchtrack.globalIndex(), cache); + // LOGF(info, "glMuons_per_MCHMID.size() = %d", glMuons_per_MCHMID.size()); + + for (const auto& fwdtrack : glMuons_per_MCHMID) { + // LOGF(info, "mchtrack.globalIndex() = %d, fwdtrack.globalIndex() = %d, mchtrack.collisionId() = %d, fwdtrack.collisionId() = %d, fwdtrack.trackType() = %d", mchtrack.globalIndex(), fwdtrack.globalIndex(), mchtrack.collisionId(), fwdtrack.collisionId(), fwdtrack.trackType()); + if (!fwdtrack.has_mcParticle()) { + continue; + } + // auto mcParticle_MFTMCHMID = fwdtrack.template mcParticle_as(); // this is identical to mcParticle_MCHMID + if (fwdtrack.chi2MatchMCHMFT() < 0.f || maxMatchingChi2MCHMFT < fwdtrack.chi2MatchMCHMFT()) { + continue; + } + if (fwdtrack.chi2MatchMCHMID() < 0.f) { + continue; + } + + auto mfttrack = fwdtrack.template matchMFTTrack_as(); // MFTsa + if (mfttrack.eta() < MFTCutGroup.minEta || MFTCutGroup.maxEta < mfttrack.eta() || mfttrack.pt() < MFTCutGroup.minPt || mfttrack.chi2() < 0.f) { + continue; + } + if (!mfttrack.has_mcParticle()) { + continue; + } + auto mcParticle_MFT = mfttrack.template mcParticle_as(); + bool isPrimary_MFT = mcParticle_MFT.isPhysicalPrimary() || mcParticle_MFT.producedByGenerator(); + + bool isMatched = (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) && (mcParticle_MFT.mcCollisionId() == mcParticle_MCHMID.mcCollisionId()); + + auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwdShift(mfttrack, mZShift, mfttrackcov); // values at innermost update + mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane + + float signed1PtMFTatMP = mftsaAtMP.getInvQPt(); + float tglMFTatMP = mftsaAtMP.getTanl(); + float phiMFTatMP = RecoDecay::constrainAngle(mftsaAtMP.getPhi(), 0, 1U); + float xMFTatMP = mftsaAtMP.getX(); + float yMFTatMP = mftsaAtMP.getY(); + + float tglErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Tanl()); + float phiErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Phi()); + float xErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2X()); + float yErrMFTatMP = std::sqrt(mftsaAtMP.getSigma2Y()); + + trackTable(collision.posZ(), collision.multFT0C(), collision.ft0cOccupancyInTimeRange(), hadronicRate, multMFT, + signed1PtMFTatMP, tglMFTatMP, phiMFTatMP, xMFTatMP, yMFTatMP, + signed1PtMCHMIDatMP, tglMCHMIDatMP, phiMCHMIDatMP, xMCHMIDatMP, yMCHMIDatMP, + fwdtrack.chi2MatchMCHMFT(), + mcParticle_MFT.pdgCode(), isPrimary_MFT, + mcParticle_MCHMID.pdgCode(), isPrimary_MCHMID, + isMatched); + + trackErrTable(tglErrMFTatMP, phiErrMFTatMP, xErrMFTatMP, yErrMFTatMP, + tglErrMCHMIDatMP, phiErrMCHMIDatMP, xErrMCHMIDatMP, yErrMCHMIDatMP); + + if (isMatched) { + fRegistry.fill(HIST("MCHMID/correct/hPt"), muonAtMP.getPt()); + fRegistry.fill(HIST("MCHMID/correct/hPz"), muonAtMP.getPz()); + fRegistry.fill(HIST("MCHMID/correct/hP"), muonAtMP.getP()); + fRegistry.fill(HIST("MCHMID/correct/hEtaPhi"), RecoDecay::constrainAngle(muonAtMP.getPhi(), 0, 1U), muonAtMP.getEta()); + fRegistry.fill(HIST("MCHMID/correct/hPEta"), muonAtMP.getEta(), muonAtMP.getP()); + fRegistry.fill(HIST("MCHMID/correct/hPtEta"), muonAtMP.getEta(), muonAtMP.getPt()); + fRegistry.fill(HIST("MCHMID/correct/hPzEta"), muonAtMP.getEta(), muonAtMP.getPz()); + + fRegistry.fill(HIST("MFT/correct/hPt"), mftsaAtMP.getPt()); + fRegistry.fill(HIST("MFT/correct/hPz"), mftsaAtMP.getPz()); + fRegistry.fill(HIST("MFT/correct/hP"), mftsaAtMP.getP()); + fRegistry.fill(HIST("MFT/correct/hEtaPhi"), RecoDecay::constrainAngle(mftsaAtMP.getPhi(), 0, 1U), mftsaAtMP.getEta()); + fRegistry.fill(HIST("MFT/correct/hPEta"), mftsaAtMP.getEta(), mftsaAtMP.getP()); + fRegistry.fill(HIST("MFT/correct/hPtEta"), mftsaAtMP.getEta(), mftsaAtMP.getPt()); + fRegistry.fill(HIST("MFT/correct/hPzEta"), mftsaAtMP.getEta(), mftsaAtMP.getPz()); + } else { + fRegistry.fill(HIST("MCHMID/wrong/hPt"), muonAtMP.getPt()); + fRegistry.fill(HIST("MCHMID/wrong/hPz"), muonAtMP.getPz()); + fRegistry.fill(HIST("MCHMID/wrong/hP"), muonAtMP.getP()); + fRegistry.fill(HIST("MCHMID/wrong/hEtaPhi"), RecoDecay::constrainAngle(muonAtMP.getPhi(), 0, 1U), muonAtMP.getEta()); + fRegistry.fill(HIST("MCHMID/wrong/hPEta"), muonAtMP.getEta(), muonAtMP.getP()); + fRegistry.fill(HIST("MCHMID/wrong/hPtEta"), muonAtMP.getEta(), muonAtMP.getPt()); + fRegistry.fill(HIST("MCHMID/wrong/hPzEta"), muonAtMP.getEta(), muonAtMP.getPz()); + + fRegistry.fill(HIST("MFT/wrong/hPt"), mftsaAtMP.getPt()); + fRegistry.fill(HIST("MFT/wrong/hPz"), mftsaAtMP.getPz()); + fRegistry.fill(HIST("MFT/wrong/hP"), mftsaAtMP.getP()); + fRegistry.fill(HIST("MFT/wrong/hEtaPhi"), RecoDecay::constrainAngle(mftsaAtMP.getPhi(), 0, 1U), mftsaAtMP.getEta()); + fRegistry.fill(HIST("MFT/wrong/hPEta"), mftsaAtMP.getEta(), mftsaAtMP.getP()); + fRegistry.fill(HIST("MFT/wrong/hPtEta"), mftsaAtMP.getEta(), mftsaAtMP.getPt()); + fRegistry.fill(HIST("MFT/wrong/hPzEta"), mftsaAtMP.getEta(), mftsaAtMP.getPz()); + } + } + } // end of MCHMID track - } // end of fwdtrack loop } // end of collision loop + + map_mfttrackcovs.clear(); } - PROCESS_SWITCH(TreeCreatorMuonML, processWithoutMFTCov, "produce ML input for single track level", false); + PROCESS_SWITCH(TreeCreatorMuonML, processWithMFTCov, "produce ML input for single track level", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)