diff --git a/Common/Core/fwdtrackUtilities.h b/Common/Core/fwdtrackUtilities.h index 29a9db2c194..932a77d9ae3 100644 --- a/Common/Core/fwdtrackUtilities.h +++ b/Common/Core/fwdtrackUtilities.h @@ -108,16 +108,16 @@ o2::track::TrackParCovFwd getTrackParCovFwdShift(TFwdTrack const& track, float z return getTrackParCovFwd3DShift(track, 0.f, 0.f, zshift, covOpt...); } -inline o2::track::TrackParCovFwd getTrackParCovFwdShiftManual( +inline o2::track::TrackParCovFwd getTrackParCovFwd3DShiftManual( const double x, const double y, const double phi, const double tgl, const double signed1Pt, const double cXX, const double cXY, const double cYY, const double cPhiX, const double cPhiY, const double cPhiPhi, const double cTglX, const double cTglY, const double cTglPhi, const double cTglTgl, const double c1PtX, const double c1PtY, const double c1PtPhi, const double c1PtTgl, const double c1Pt21Pt2, - const float z, const float zshift, const float chi2) + const float z, const float xshift, const float yshift, const float zshift, const float chi2) { - SMatrix5 tpars(x, y, phi, tgl, signed1Pt); + SMatrix5 tpars(x + xshift, y + yshift, phi, tgl, signed1Pt); SMatrix55 tcovs; std::vector v1{ @@ -164,15 +164,15 @@ o2::track::TrackParCovFwd getTrackParCovFwd(TFwdTrack const& track, TFwdTrackCov /// propagate fwdtrack to a certain point. template -o2::dataformats::GlobalFwdTrack propagateMuon(TFwdTrack const& muon, TFwdTrackCov const& cov, TCollision const& collision, const propagationPoint endPoint, const float matchingZ, const float bzkG, const float zshift = 0.f) +o2::dataformats::GlobalFwdTrack propagateMuon(TFwdTrack const& muon, TFwdTrackCov const& cov, TCollision const& collision, const propagationPoint endPoint, const float matchingZ, const float bzkG, const float xshift = 0.f, const float yshift = 0.f, const float zshift = 0.f) { o2::track::TrackParCovFwd trackParCovFwd; if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - trackParCovFwd = getTrackParCovFwdShift(muon, zshift, cov); + trackParCovFwd = getTrackParCovFwd3DShift(muon, xshift, yshift, zshift, cov); } else if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - trackParCovFwd = getTrackParCovFwdShift(muon, zshift, muon); + trackParCovFwd = getTrackParCovFwd3DShift(muon, xshift, yshift, zshift, muon); } else { - trackParCovFwd = getTrackParCovFwdShift(muon, zshift, muon); + trackParCovFwd = getTrackParCovFwd3DShift(muon, xshift, yshift, zshift, muon); } o2::dataformats::GlobalFwdTrack propmuon = propagateTrackParCovFwd(trackParCovFwd, muon.trackType(), collision, endPoint, matchingZ, bzkG); @@ -340,7 +340,7 @@ float getFwdChi2IP(TTrackParCovFwd const& inputTrk, TCollision const& collision, } template -float getFwdChi2IP(const TFullFwdTrack& fwdtrack, const TCollision& collision, const float bz, const float zShift) +float getFwdChi2IP(const TFullFwdTrack& fwdtrack, const TCollision& collision, const float bz, const float xShift, const float yShift, const float zShift) { // this function returns imcompatibility of fwdtrack respect to a given PV. // fwdtracks are never PV contributors in ALICE. @@ -348,7 +348,7 @@ float getFwdChi2IP(const TFullFwdTrack& fwdtrack, const TCollision& collision, c // chi2IP cannot be used to decide the best fwdtrack-to-collision match or MFT-MCH match, because it gives biases toward small muon impact parameter. // chi2IP should be used only after the best fwdtrack-to-collision association and the best MFT-MCH match are defined. - o2::track::TrackParCovFwd trk = getTrackParCovFwdShift(fwdtrack, zShift, fwdtrack); + o2::track::TrackParCovFwd trk = getTrackParCovFwd3DShift(fwdtrack, xShift, yShift, zShift, fwdtrack); if (std::abs(bz) < 1e-12) { trk.propagateToZlinear(collision.posZ()); diff --git a/EventFiltering/PWGEM/globalDimuonFilter.cxx b/EventFiltering/PWGEM/globalDimuonFilter.cxx index 48cf5515f5c..8d5e8f2ae8d 100644 --- a/EventFiltering/PWGEM/globalDimuonFilter.cxx +++ b/EventFiltering/PWGEM/globalDimuonFilter.cxx @@ -147,13 +147,17 @@ struct globalDimuonFilter { // for z shift for propagation o2::framework::Configurable cfgApplyZShiftFromCCDB{"cfgApplyZShiftFromCCDB", false, "flag to apply z shift"}; o2::framework::Configurable cfgZShiftPath{"cfgZShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"}; - o2::framework::Configurable cfgManualZShift{"cfgManualZShift", 0, "manual z-shift for propagation of global muon to PV"}; + o2::framework::Configurable cfgManualXShift{"cfgManualXShift", 0, "manual x shift for propagation of global muon to PV"}; + o2::framework::Configurable cfgManualYShift{"cfgManualYShift", 0, "manual y shift for propagation of global muon to PV"}; + o2::framework::Configurable cfgManualZShift{"cfgManualZShift", 0, "manual z shift for propagation of global muon to PV"}; o2::framework::HistogramRegistry fRegistry{"output", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, false}; o2::ccdb::CcdbApi ccdbApi; o2::framework::Service ccdb; int mRunNumber = 0; float mBz = 0; + float mXShift = 0; + float mYShift = 0; float mZShift = 0; void init(o2::framework::InitContext&) @@ -165,6 +169,8 @@ struct globalDimuonFilter { ccdbApi.init(ccdburl); mRunNumber = 0; mBz = 0; + mXShift = 0; + mYShift = 0; mZShift = 0; addHistograms(); @@ -203,6 +209,8 @@ struct globalDimuonFilter { } } else { LOGF(info, "z shift is manually set to %f cm", cfgManualZShift.value); + mXShift = cfgManualXShift; + mYShift = cfgManualYShift; mZShift = cfgManualZShift; } } @@ -418,12 +426,12 @@ struct globalDimuonFilter { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, glMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, glMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); float etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); phiMatchedMCHMID = RecoDecay::constrainAngle(phiMatchedMCHMID, 0, 1U); - o2::dataformats::GlobalFwdTrack propmuonAtPV = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, glMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, glMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); pt = propmuonAtPV.getPt(); eta = propmuonAtPV.getEta(); phi = propmuonAtPV.getPhi(); @@ -440,7 +448,7 @@ struct globalDimuonFilter { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtDCA = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, glMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtDCA = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, glMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); float dcaX = propmuonAtDCA.getX() - collision.posX(); float dcaY = propmuonAtDCA.getY() - collision.posY(); float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); @@ -474,7 +482,7 @@ struct globalDimuonFilter { } float sigma_dcaXY = dcaXY / dcaXYinSigma / std::sqrt(2.f); - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, glMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, glMuonCutGroup.matchingZ, mBz, mXShift, mYShift, 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); @@ -488,7 +496,7 @@ struct globalDimuonFilter { return false; } - float chi2IP = o2::aod::fwdtrackutils::getFwdChi2IP(fwdtrack, collision, mBz, mZShift); + float chi2IP = o2::aod::fwdtrackutils::getFwdChi2IP(fwdtrack, collision, mBz, mXShift, mYShift, mZShift); if constexpr (fillHistograms) { if (fwdtrack.sign() > 0) { @@ -605,12 +613,12 @@ struct globalDimuonFilter { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, tagMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, tagMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); float etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); phiMatchedMCHMID = RecoDecay::constrainAngle(phiMatchedMCHMID, 0, 1U); - o2::dataformats::GlobalFwdTrack propmuonAtPV = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, tagMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, tagMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); pt = propmuonAtPV.getPt(); eta = propmuonAtPV.getEta(); phi = propmuonAtPV.getPhi(); @@ -693,12 +701,12 @@ struct globalDimuonFilter { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, probeMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, probeMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); float etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); phiMatchedMCHMID = RecoDecay::constrainAngle(phiMatchedMCHMID, 0, 1U); - o2::dataformats::GlobalFwdTrack propmuonAtPV = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, probeMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV = o2::aod::fwdtrackutils::propagateMuon(fwdtrack, fwdtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToVertex, probeMuonCutGroup.matchingZ, mBz, mXShift, mYShift, mZShift); pt = propmuonAtPV.getPt(); eta = propmuonAtPV.getEta(); phi = propmuonAtPV.getPhi(); @@ -722,7 +730,7 @@ struct globalDimuonFilter { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, probeMuonCutGroup.matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = o2::aod::fwdtrackutils::propagateMuon(mchtrack, mchtrack, collision, o2::aod::fwdtrackutils::propagationPoint::kToDCA, probeMuonCutGroup.matchingZ, mBz, mXShift, mYShift, 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); diff --git a/PWGEM/Dilepton/TableProducer/CMakeLists.txt b/PWGEM/Dilepton/TableProducer/CMakeLists.txt index 5d2889041ae..8a00046766a 100644 --- a/PWGEM/Dilepton/TableProducer/CMakeLists.txt +++ b/PWGEM/Dilepton/TableProducer/CMakeLists.txt @@ -46,11 +46,6 @@ o2physics_add_dpl_workflow(skimmer-primary-muon PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking COMPONENT_NAME Analysis) -o2physics_add_dpl_workflow(skimmer-primary-muon-qc - SOURCES skimmerPrimaryMuonQC.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::GlobalTracking - COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(skimmer-primary-track SOURCES skimmerPrimaryTrack.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index 8706432f69d..27ad0708be0 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -112,16 +112,26 @@ struct skimmerPrimaryMuon { // for z shift for propagation Configurable cfgApplyZShiftFromCCDB{"cfgApplyZShiftFromCCDB", false, "flag to apply z shift"}; Configurable cfgZShiftPath{"cfgZShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"}; - Configurable cfgManualZShift{"cfgManualZShift", 0, "manual z-shift for propagation of global muon to PV"}; + Configurable cfgManualXShiftMFTtop{"cfgManualXShiftMFTtop", 0, "manual x shift on MFT top when propagating global muon to PV"}; + Configurable cfgManualYShiftMFTtop{"cfgManualYShiftMFTtop", 0, "manual y shift on MFT top when propagating global muon to PV"}; + Configurable cfgManualZShiftMFTtop{"cfgManualZShiftMFTtop", 0, "manual z shift on MFT top when propagating global muon to PV"}; + Configurable cfgManualXShiftMFTbottom{"cfgManualXShiftMFTbottom", 0, "manual x shift on MFT bottom when propagating global muon to PV"}; + Configurable cfgManualYShiftMFTbottom{"cfgManualYShiftMFTbottom", 0, "manual y shift on MFT bottom when propagating global muon to PV"}; + Configurable cfgManualZShiftMFTbottom{"cfgManualZShiftMFTbottom", 0, "manual z shift on MFT bottom when propagating global muon to PV"}; o2::ccdb::CcdbApi ccdbApi; Service ccdb; int mRunNumber = 0; float mBz = 0; - float mZShift = 0; + float mXShiftMFTtop = 0; + float mYShiftMFTtop = 0; + float mZShiftMFTtop = 0; + float mXShiftMFTbottom = 0; + float mYShiftMFTbottom = 0; + float mZShiftMFTbottom = 0; HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - static constexpr std::string_view muon_types[5] = {"MFTMCHMID/", "MFTMCHMIDOtherMatch/", "MFTMCH/", "MCHMID/", "MCH/"}; + // static constexpr std::string_view muon_types[5] = {"MFTMCHMID/", "MFTMCHMIDOtherMatch/", "MFTMCH/", "MCHMID/", "MCH/"}; void init(InitContext&) { @@ -140,7 +150,12 @@ struct skimmerPrimaryMuon { } mRunNumber = 0; mBz = 0; - mZShift = 0; + mXShiftMFTtop = 0; + mYShiftMFTtop = 0; + mZShiftMFTtop = 0; + mXShiftMFTbottom = 0; + mYShiftMFTbottom = 0; + mZShiftMFTbottom = 0; } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -168,14 +183,22 @@ struct skimmerPrimaryMuon { auto* zShift = ccdb->getForTimeStamp>(cfgZShiftPath, bc.timestamp()); if (zShift != nullptr && !zShift->empty()) { LOGF(info, "reading z shift %f from %s", (*zShift)[0], cfgZShiftPath.value); - mZShift = (*zShift)[0]; + mZShiftMFTtop = (*zShift)[0]; + mZShiftMFTbottom = (*zShift)[0]; } else { LOGF(info, "z shift is not found in ccdb path %s. set to 0 cm", cfgZShiftPath.value); - mZShift = 0; + mZShiftMFTtop = 0; + mZShiftMFTbottom = 0; } } else { - LOGF(info, "z shift is manually set to %f cm", cfgManualZShift.value); - mZShift = cfgManualZShift; + LOGF(info, "z shift on MFT top is manually set to %f cm", cfgManualZShiftMFTtop.value); + LOGF(info, "z shift on MFT bottom is manually set to %f cm", cfgManualZShiftMFTbottom.value); + mXShiftMFTtop = cfgManualXShiftMFTtop; + mYShiftMFTtop = cfgManualYShiftMFTtop; + mZShiftMFTtop = cfgManualZShiftMFTtop; + mXShiftMFTbottom = cfgManualXShiftMFTbottom; + mYShiftMFTbottom = cfgManualYShiftMFTbottom; + mZShiftMFTbottom = cfgManualZShiftMFTbottom; } } @@ -209,8 +232,10 @@ struct skimmerPrimaryMuon { fRegistry.add("MFTMCHMID/hDCAxyinSigma", "DCAxy in sigma;DCA_{xy} (#sigma);", kTH1F, {{100, 0, 10}}, false); fRegistry.add("MFTMCHMID/hLog10Chi2IP", "chi2IP;log_{10}(#chi^{2}_{IP})", kTH1F, {{100, -5, 5}}, false); fRegistry.add("MFTMCHMID/hSqrtChi2IP", "chi2IP;#sqrt{#chi^{2}_{IP}}", 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_PosZ_MFTtop", "DCAx vs. posZ;Z_{vtx} (cm);DCA_{x} (cm)", kTH2F, {{200, -10, +10}, {400, -0.2, +0.2}}, false); + fRegistry.add("MFTMCHMID/hDCAy_PosZ_MFTtop", "DCAy vs. posZ;Z_{vtx} (cm);DCA_{y} (cm)", kTH2F, {{200, -10, +10}, {400, -0.2, +0.2}}, false); + fRegistry.add("MFTMCHMID/hDCAx_PosZ_MFTbottom", "DCAx vs. posZ;Z_{vtx} (cm);DCA_{x} (cm)", kTH2F, {{200, -10, +10}, {400, -0.2, +0.2}}, false); + fRegistry.add("MFTMCHMID/hDCAy_PosZ_MFTbottom", "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, {{180, -M_PI, M_PI}, {400, -0.2, +0.2}}, false); fRegistry.add("MFTMCHMID/hDCAy_Phi", "DCAy vs. #varphi;#varphi (rad.);DCA_{y} (cm)", kTH2F, {{180, -M_PI, M_PI}, {400, -0.2, +0.2}}, false); fRegistry.add("MFTMCHMID/hMeanDCAx", ";X_{IU} (cm);Y_{IU} (cm); (cm)", kTProfile2D, {{240, -12, +12}, {240, -12, +12}}, false); @@ -279,51 +304,24 @@ struct skimmerPrimaryMuon { 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); + o2::dataformats::GlobalFwdTrack propmuonAtPV; + float pt = 0.f, eta = 0.f, phi = 0.f; - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift); - - float dcaX = propmuonAtDCA.getX() - collision.posX(); - float dcaY = propmuonAtDCA.getY() - collision.posY(); - float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack - float cXXatDCA = propmuonAtDCA.getSigma2X(); - float cYYatDCA = propmuonAtDCA.getSigma2Y(); - float cXYatDCA = propmuonAtDCA.getSigmaXY(); - - float det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat + o2::dataformats::GlobalFwdTrack propmuonAtDCA; + float dcaX = 999.f, dcaY = 999.f, dcaXY = 999.f; + float rAtAbsorberEnd = 1e+10; + float cXXatDCA = 1e+10, cYYatDCA = 1e+10, cXYatDCA = 1e+10; float dcaXYinSigma = 999.f; - if (det < 0) { - dcaXYinSigma = 999.f; - } else { - dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2.f * dcaX * dcaY * cXYatDCA) / det / 2.f)); // dca xy in sigma - } - float sigma_dcaXY = dcaXY / dcaXYinSigma; + float sigma_dcaXY = 1e+10; - float pDCA = propmuonAtPV.getP() * dcaXY; + float pDCA = 1e+10; int nClustersMFT = 0; - float ptMatchedMCHMID = propmuonAtPV.getPt(); - float etaMatchedMCHMID = propmuonAtPV.getEta(); - float phiMatchedMCHMID = propmuonAtPV.getPhi(); - o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - // float x = fwdtrack.x(); - // float y = fwdtrack.y(); - // float z = fwdtrack.z(); - // float tgl = fwdtrack.tgl(); + float ptMatchedMCHMID = 1e+10, etaMatchedMCHMID = 1e+10, phiMatchedMCHMID = 1e+10; float chi2mft = 0.f; uint64_t mftClusterSizesAndTrackFlags = 0; int ndf_mchmft = 1; int ndf_mft = 1; - // float etaMatchedMCHMIDatMP = 999.f; - float phiMatchedMCHMIDatMP = 999.f; - // float etaMatchedMFTatMP = 999.f; - float phiMatchedMFTatMP = 999.f; - float deta = 999.f; float dphi = 999.f; @@ -341,17 +339,40 @@ struct skimmerPrimaryMuon { return false; } - // apply dca cut here to minimize the number of calling propagateMuon. - if (maxDCAxy < dcaXY) { - return false; - } - auto mchtrack = fwdtrack.template matchMCHTrack_as(); // MCH-MID auto mfttrack = fwdtrack.template matchMFTTrack_as(); // MFTsa if (mfttrack.chi2() < 0.f) { return false; } + propmuonAtPV = std::atan2(mfttrack.y(), mfttrack.x()) > 0.f ? propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mXShiftMFTtop, mYShiftMFTtop, mZShiftMFTtop) : propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mXShiftMFTbottom, mYShiftMFTbottom, mZShiftMFTbottom); + pt = propmuonAtPV.getPt(); + eta = propmuonAtPV.getEta(); + phi = propmuonAtPV.getPhi(); + o2::math_utils::bringTo02Pi(phi); + + propmuonAtDCA = std::atan2(mfttrack.y(), mfttrack.x()) > 0.f ? propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mXShiftMFTtop, mYShiftMFTtop, mZShiftMFTtop) : propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mXShiftMFTbottom, mYShiftMFTbottom, mZShiftMFTbottom); + dcaX = propmuonAtDCA.getX() - collision.posX(); + dcaY = propmuonAtDCA.getY() - collision.posY(); + dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + cXXatDCA = propmuonAtDCA.getSigma2X(); + cYYatDCA = propmuonAtDCA.getSigma2Y(); + cXYatDCA = propmuonAtDCA.getSigmaXY(); + rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack + + float det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinant + if (det < 0) { + dcaXYinSigma = 999.f; + } else { + dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2.f * dcaX * dcaY * cXYatDCA) / det / 2.f)); // dca xy in sigma + } + sigma_dcaXY = dcaXY / dcaXYinSigma / std::sqrt(2); + + // apply dca cut here to minimize the number of calling propagateMuon. + if (maxDCAxy < dcaXY) { + return false; + } + xMFT = mfttrack.x(); yMFT = mfttrack.y(); @@ -369,39 +390,30 @@ struct skimmerPrimaryMuon { ndf_mchmft = 2.f * (mchtrack.nClusters() + nClustersMFT) - 5.f; ndf_mft = 2.f * nClustersMFT - 5.f; chi2mft = mfttrack.chi2(); - // chi2mft = mfttrack.chi2() / (2.f * nClustersMFT - 5.f); // apply chi2/ndf cut here to minimize the number of calling propagateMuon. if (maxChi2GL < fwdtrack.chi2() / ndf_mchmft) { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, 0.f, 0.f, 0.f); ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - 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); - pDCA = mchtrack.p() * dcaXY_Matched; + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, 0.f, 0.f, 0.f); + // float dcaXY_Matched = std::hypot(propmuonAtDCA_Matched.getX() - collision.posX(), propmuonAtDCA_Matched.getY() - collision.posY()); + pDCA = mchtrack.p() * std::hypot(propmuonAtDCA_Matched.getX() - collision.posX(), propmuonAtDCA_Matched.getY() - collision.posY()); if constexpr (withMFTCov) { auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz, mZShift); // propagated to matching plane - o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwdShift(mfttrack, mZShift, mfttrackcov); // values at innermost update - mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane - // etaMatchedMFTatMP = mftsaAtMP.getEta(); - phiMatchedMFTatMP = mftsaAtMP.getPhi(); - // etaMatchedMCHMIDatMP = muonAtMP.getEta(); - phiMatchedMCHMIDatMP = muonAtMP.getPhi(); - o2::math_utils::bringTo02Pi(phiMatchedMCHMIDatMP); - o2::math_utils::bringTo02Pi(phiMatchedMFTatMP); - - o2::track::TrackParCovFwd mftsa = getTrackParCovFwdShift(mfttrack, mZShift, mfttrackcov); // values at innermost update - o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. + // auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz, 0.f, 0.f, 0.f); // propagated to matching plane + // o2::track::TrackParCovFwd mftsaAtMP = std::atan2(mfttrack.y(), mfttrack.x()) > 0.f ? getTrackParCovFwd3DShift(mfttrack, mXShiftMFTtop, mYShiftMFTtop, mZShiftMFTtop, mfttrackcov) : getTrackParCovFwd3DShift(mfttrack, mXShiftMFTbottom, mYShiftMFTbottom, mZShiftMFTbottom, mfttrackcov); // values at innermost update + // mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane + + o2::track::TrackParCovFwd mftsa = std::atan2(mfttrack.y(), mfttrack.x()) > 0.f ? getTrackParCovFwd3DShift(mfttrack, mXShiftMFTtop, mYShiftMFTtop, mZShiftMFTtop, mfttrackcov) : getTrackParCovFwd3DShift(mfttrack, mXShiftMFTbottom, mYShiftMFTbottom, mZShiftMFTbottom, mfttrackcov); // values at innermost update + o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. auto globalMuon = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToVertex, matchingZ, mBz); pt = globalMuon.getPt(); eta = globalMuon.getEta(); @@ -414,8 +426,7 @@ struct skimmerPrimaryMuon { dcaX = globalMuon.getX() - collision.posX(); dcaY = globalMuon.getY() - collision.posY(); dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat - dcaXYinSigma = 999.f; + det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinant if (det < 0) { dcaXYinSigma = 999.f; } else { @@ -428,7 +439,7 @@ struct skimmerPrimaryMuon { dphi = phiMatchedMCHMID - phi; o2::math_utils::bringToPMPi(dphi); - chi2IP = getFwdChi2IP(fwdtrack, collision, mBz, mZShift); + chi2IP = std::atan2(mfttrack.y(), mfttrack.x()) > 0.f ? getFwdChi2IP(fwdtrack, collision, mBz, mXShiftMFTtop, mYShiftMFTtop, mZShiftMFTtop) : getFwdChi2IP(fwdtrack, collision, mBz, mXShiftMFTbottom, mYShiftMFTbottom, mZShiftMFTbottom); if (std::sqrt(std::pow(deta / maxDEta, 2) + std::pow(dphi / maxDPhi, 2)) > 1.f) { return false; @@ -440,14 +451,14 @@ struct skimmerPrimaryMuon { const auto& fwdcov = propmuonAtPV.getCovariances(); // covatiance matrix at PV - auto globalMuonManual = getTrackParCovFwdShiftManual( + auto globalMuonManual = getTrackParCovFwd3DShiftManual( propmuonAtPV.getX(), propmuonAtPV.getY(), propmuonAtPV.getPhi(), propmuonAtPV.getTgl(), propmuonAtPV.getInvQPt() / sfPt, fwdcov(0, 0), fwdcov(1, 0), fwdcov(1, 1), fwdcov(2, 0), fwdcov(2, 1), fwdcov(2, 2), fwdcov(3, 0), fwdcov(3, 1), fwdcov(3, 2), fwdcov(3, 3), fwdcov(4, 0) / sfPt, fwdcov(4, 1) / sfPt, fwdcov(4, 2) / sfPt, fwdcov(4, 3) / sfPt, fwdcov(4, 4) / sfPt / sfPt, - propmuonAtPV.getZ(), 0.0, fwdtrack.chi2()); + propmuonAtPV.getZ(), 0.0, 0.0, 0.0, fwdtrack.chi2()); auto globalMuonManualAtZPV = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonManual, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); @@ -471,24 +482,29 @@ struct skimmerPrimaryMuon { cXYatDCA = cXYatDCA + resPVXY; cYYatDCA = cYYatDCA + resPVYY; - det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat - dcaXYinSigma = 999.f; + det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinant if (det < 0) { dcaXYinSigma = 999.f; } else { dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2.f * dcaX * dcaY * cXYatDCA) / det / 2.f)); // dca xy in sigma } - sigma_dcaXY = dcaXY / dcaXYinSigma; + sigma_dcaXY = dcaXY / dcaXYinSigma / std::sqrt(2); chi2IP = getFwdChi2IP(globalMuonManual, collision, mBz); } } else if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToRabs, matchingZ, mBz, mZShift); // this is necessary only for MuonStandaloneTrack + propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, 0.f, 0.f, 0.f); + pt = propmuonAtPV.getPt(); + eta = propmuonAtPV.getEta(); + phi = propmuonAtPV.getPhi(); + o2::math_utils::bringTo02Pi(phi); + + o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToRabs, matchingZ, mBz, 0.f, 0.f, 0.f); // this is necessary only for MuonStandaloneTrack float xAbs = propmuonAtRabs.getX(); float yAbs = propmuonAtRabs.getY(); rAtAbsorberEnd = std::sqrt(xAbs * xAbs + yAbs * yAbs); // Redo propagation only for muon tracks // propagation of MFT tracks alredy done in reconstruction - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift); + propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, 0.f, 0.f, 0.f); cXXatDCA = propmuonAtDCA.getSigma2X(); cYYatDCA = propmuonAtDCA.getSigma2Y(); cXYatDCA = propmuonAtDCA.getSigmaXY(); @@ -497,14 +513,13 @@ struct skimmerPrimaryMuon { dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); pDCA = fwdtrack.p() * dcaXY; - det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat - dcaXYinSigma = 999.f; + float det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinant if (det < 0) { dcaXYinSigma = 999.f; } else { dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2.f * dcaX * dcaY * cXYatDCA) / det / 2.f)); // dca xy in sigma } - sigma_dcaXY = dcaXY / dcaXYinSigma; + sigma_dcaXY = dcaXY / dcaXYinSigma / std::sqrt(2); } else { return false; } @@ -516,10 +531,6 @@ struct skimmerPrimaryMuon { if constexpr (fillTable) { float dpt = (ptMatchedMCHMID - pt) / pt; - // float detaMP = etaMatchedMCHMIDatMP - etaMatchedMFTatMP; - // float dphiMP = phiMatchedMCHMIDatMP - phiMatchedMFTatMP; - // o2::math_utils::bringToPMPi(dphiMP); - bool isAssociatedToMPC = fwdtrack.collisionId() == collision.globalIndex(); // LOGF(info, "isAmbiguous = %d, isAssociatedToMPC = %d, fwdtrack.globalIndex() = %d, fwdtrack.collisionId() = %d, collision.globalIndex() = %d", isAmbiguous, isAssociatedToMPC, fwdtrack.globalIndex(), fwdtrack.collisionId(), collision.globalIndex()); @@ -572,8 +583,13 @@ struct skimmerPrimaryMuon { fRegistry.fill(HIST("MFTMCHMID/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/hLog10Chi2IP"), std::log10(chi2IP)); fRegistry.fill(HIST("MFTMCHMID/hSqrtChi2IP"), std::sqrt(chi2IP)); - fRegistry.fill(HIST("MFTMCHMID/hDCAx_PosZ"), collision.posZ(), dcaX); - fRegistry.fill(HIST("MFTMCHMID/hDCAy_PosZ"), collision.posZ(), dcaY); + if (std::atan2(yMFT, xMFT) > 0.f) { + fRegistry.fill(HIST("MFTMCHMID/hDCAx_PosZ_MFTtop"), collision.posZ(), dcaX); + fRegistry.fill(HIST("MFTMCHMID/hDCAy_PosZ_MFTtop"), collision.posZ(), dcaY); + } else { + fRegistry.fill(HIST("MFTMCHMID/hDCAx_PosZ_MFTbottom"), collision.posZ(), dcaX); + fRegistry.fill(HIST("MFTMCHMID/hDCAy_PosZ_MFTbottom"), collision.posZ(), dcaY); + } fRegistry.fill(HIST("MFTMCHMID/hDCAx_Phi"), std::atan2(yMFT, xMFT), dcaX); fRegistry.fill(HIST("MFTMCHMID/hDCAy_Phi"), std::atan2(yMFT, xMFT), dcaY); fRegistry.fill(HIST("MFTMCHMID/hMeanDCAx"), fwdtrack.x(), fwdtrack.y(), dcaX); @@ -627,17 +643,6 @@ struct skimmerPrimaryMuon { // LOGF(info, "stanadalone: muon.globalIndex() = %d, muon.chi2MatchMCHMFT() = %f", muon.globalIndex(), muon.chi2MatchMCHMFT()); // LOGF(info, "muons_per_MCHMID.size() = %d", muons_per_MCHMID.size()); - // o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); - // float etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); - // float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); - // o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - - // o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(fwdtrack, fwdtrack, 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 = fwdtrack.p() * dcaXY_Matched; - float min_chi2MatchMCHMFT = 1e+10; std::tuple tupleId_at_min_chi2mftmch; std::vector vec_chi2tmp; @@ -653,41 +658,6 @@ struct skimmerPrimaryMuon { continue; } - // if (muon_tmp.chi2() < 0.f || muon_tmp.chi2MatchMCHMFT() < 0.f || muon_tmp.chi2MatchMCHMID() < 0.f || mfttrack.chi2() < 0.f) { // reject negative chi2, i.e. wrong. - // continue; - // } - - // o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(muon_tmp, muon_tmp, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); - // float pt = propmuonAtPV.getPt(); - // float eta = propmuonAtPV.getEta(); - // float phi = propmuonAtPV.getPhi(); - // o2::math_utils::bringTo02Pi(phi); - - // if (refitGlobalMuon) { - // pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); - // } - - // float deta = etaMatchedMCHMID - eta; - // float dphi = phiMatchedMCHMID - phi; - // o2::math_utils::bringToPMPi(dphi); - // int ndf = 2 * (mchtrack.nClusters() + mfttrack.nClusters()) - 5; - - // float dcaX = propmuonAtPV.getX() - collision.posX(); - // float dcaY = propmuonAtPV.getY() - collision.posY(); - // float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - - // if (cfgApplyPreselectionInBestMatch) { - // if (!isSelected(pt, eta, muon_tmp.rAtAbsorberEnd(), pDCA, muon_tmp.chi2() / ndf, muon_tmp.trackType(), dcaXY)) { - // continue; - // } - // if (std::sqrt(std::pow(deta / maxDEta, 2) + std::pow(dphi / maxDPhi, 2)) > 1.f) { - // continue; - // } - // if (muon_tmp.chi2MatchMCHMFT() > maxMatchingChi2MCHMFT) { - // continue; - // } - // } - vec_chi2tmp.emplace_back(muon_tmp.chi2MatchMCHMFT()); if (0.f < muon_tmp.chi2MatchMCHMFT() && muon_tmp.chi2MatchMCHMFT() < min_chi2MatchMCHMFT) { min_chi2MatchMCHMFT = muon_tmp.chi2MatchMCHMFT(); diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuonQC.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuonQC.cxx deleted file mode 100644 index 246634a975d..00000000000 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuonQC.cxx +++ /dev/null @@ -1,797 +0,0 @@ -// 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 write relevant information for muons. -/// \author daiki.sekihata@cern.ch - -#include "PWGEM/Dilepton/DataModel/dileptonTables.h" - -#include "Common/Core/fwdtrackUtilities.h" -#include "Common/DataModel/CollisionAssociationTables.h" -#include "Common/DataModel/EventSelection.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) -#include -#include - -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::soa; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::constants::physics; -using namespace o2::aod::fwdtrackutils; - -struct skimmerPrimaryMuonQC { - using MyCollisions = soa::Join; - using MyCollisionsWithSWT = soa::Join; - - using MyFwdTracks = soa::Join; // muon tracks are repeated. i.e. not exclusive. - using MyFwdTrack = MyFwdTracks::iterator; - - using MyFwdTracksMC = soa::Join; - using MyFwdTrackMC = MyFwdTracksMC::iterator; - - using MFTTracksMC = soa::Join; - using MFTTrackMC = MFTTracksMC::iterator; - - Produces emprimarymuons; - Produces emprimarymuonscov; - - // Configurables - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - Configurable fillQAHistograms{"fillQAHistograms", false, "flag to fill QA histograms"}; - - // for z shift for propagation - Configurable cfgApplyZShiftFromCCDB{"cfgApplyZShiftFromCCDB", false, "flag to apply z shift"}; - Configurable cfgZShiftPath{"cfgZShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"}; - Configurable cfgManualZShift{"cfgManualZShift", 0, "manual z-shift for propagation of global muon to PV"}; - Configurable matchingZ{"matchingZ", -77.5, "z position where matching is performed"}; - Configurable refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"}; - - struct : ConfigurableGroup { // tight cut - std::string prefix = "tagMuonCut"; - Configurable minPt{"minPt", 0.8, "min pt for muon"}; - Configurable maxPt{"maxPt", 1e+10, "max pt for muon"}; - Configurable minEta{"minEta", -3.6, "min. eta acceptance for MFT-MCH-MID"}; - Configurable maxEta{"maxEta", -2.5, "max. eta acceptance for MFT-MCH-MID"}; - Configurable minRabs{"minRabs", 27.6, "min. R at absorber end for global muon (min. eta = -3.6)"}; // std::tan(2.f * std::atan(std::exp(- -3.6)) ) * -505. = 27.6 - Configurable maxRabs{"maxRabs", 89.5, "max. R at absorber end"}; - Configurable maxDCAxy{"maxDCAxy", 0.06, "max. DCAxy for global muons"}; - Configurable maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 40.f, "max. chi2 for MCH-MFT matching"}; - Configurable maxChi2{"maxChi2", 4.f, "max. chi2/ndf for global muon"}; - // Configurable minNclsMFT{"minNclsMFT", 5, "min ncluster of MFT"}; - // Configurable minNclsMCH{"minNclsMCH", 5, "min ncluster of MCH"}; - Configurable maxDEta{"maxDEta", 0.08, "max. deta between MFT-MCH-MID and MCH-MID"}; - Configurable maxDPhi{"maxDPhi", 0.08, "max. dphi between MFT-MCH-MID and MCH-MID"}; - } tagMuonCut; - - struct : ConfigurableGroup { // loose cut - std::string prefix = "probeMuonCut"; - Configurable minPt{"minPt", 0.01, "min pt for muon"}; - Configurable maxPt{"maxPt", 1e+10, "max pt for muon"}; - Configurable minEtaSA{"minEtaSA", -4.0, "min. eta acceptance for MFT-MCH"}; - Configurable maxEtaSA{"maxEtaSA", -2.5, "max. eta acceptance for MFT-MCH"}; - Configurable minEtaGL{"minEtaGL", -3.6, "min. eta acceptance for MFT-MCH-MID"}; - Configurable maxEtaGL{"maxEtaGL", -2.5, "max. eta acceptance for MFT-MCH-MID"}; - Configurable minRabs{"minRabs", 17.6, "min. R at absorber end for global muon (min. eta = -3.6)"}; // std::tan(2.f * std::atan(std::exp(- -3.6)) ) * -505. = 27.6 - Configurable midRabs{"midRabs", 26.5, "middle R at absorber end for pDCA cut"}; - Configurable maxRabs{"maxRabs", 89.5, "max. R at absorber end"}; - Configurable maxDCAxy{"maxDCAxy", 1.f, "max. DCAxy for global muons"}; - Configurable maxPDCAforLargeR{"maxPDCAforLargeR", 324.f, "max. pDCA for large R at absorber end"}; - Configurable maxPDCAforSmallR{"maxPDCAforSmallR", 594.f, "max. pDCA for small R at absorber end"}; - Configurable maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 100, "max. chi2 for MCH-MFT matching"}; - Configurable maxChi2{"maxChi2", 1e+10, "max. chi2/ndf for global muon"}; - // Configurable minNclsMFT{"minNclsMFT", 5, "min ncluster of MFT"}; - // Configurable minNclsMCH{"minNclsMCH", 5, "min ncluster of MCH"}; - Configurable maxDEta{"maxDEta", 1e+10, "max. deta between MFT-MCH-MID and MCH-MID"}; - Configurable maxDPhi{"maxDPhi", 1e+10, "max. dphi between MFT-MCH-MID and MCH-MID"}; - } probeMuonCut; - - struct : ConfigurableGroup { - std::string prefix = "pairCuts"; - Configurable minMass{"minMass", 0.21, "min mass"}; - Configurable maxMass{"maxMass", 0.30, "max mass"}; - } pairCuts; - - o2::ccdb::CcdbApi ccdbApi; - Service ccdb; - int mRunNumber = 0; - float mBz = 0; - float mZShift = 0; - - HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - // static constexpr std::string_view muon_types[5] = {"MFTMCHMID/", "MFTMCHMIDOtherMatch/", "MFTMCH/", "MCHMID/", "MCH/"}; - - void init(InitContext&) - { - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - ccdbApi.init(ccdburl); - - if (fillQAHistograms) { - addHistograms(); - } - mRunNumber = 0; - mBz = 0; - mZShift = 0; - } - - void initCCDB(aod::BCsWithTimestamps::iterator const& bc) - { - if (mRunNumber == bc.runNumber()) { - return; - } - mRunNumber = bc.runNumber(); - - std::map metadata; - auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, mRunNumber); - auto ts = soreor.first; - auto grpmag = ccdbApi.retrieveFromTFileAny(grpmagPath, metadata, ts); - o2::base::Propagator::initFieldFromGRP(grpmag); - if (!o2::base::GeometryManager::isGeometryLoaded()) { - ccdb->get(geoPath); - } - o2::mch::TrackExtrap::setField(); - const double centerMFT[3] = {0, 0, -61.4}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - mBz = field->getBz(centerMFT); // Get field at centre of MFT - LOGF(info, "Bz at center of MFT = %f kZG", mBz); - - if (cfgApplyZShiftFromCCDB) { - auto* zShift = ccdb->getForTimeStamp>(cfgZShiftPath, bc.timestamp()); - if (zShift != nullptr && !zShift->empty()) { - LOGF(info, "reading z shift %f from %s", (*zShift)[0], cfgZShiftPath.value); - mZShift = (*zShift)[0]; - } else { - LOGF(info, "z shift is not found in ccdb path %s. set to 0 cm", cfgZShiftPath.value); - mZShift = 0; - } - } else { - LOGF(info, "z shift is manually set to %f cm", cfgManualZShift.value); - mZShift = cfgManualZShift; - } - } - - void addHistograms() - { - // auto hMuonType = fRegistry.add("hMuonType", "muon type", kTH1F, {{5, -0.5f, 4.5f}}, false); - // hMuonType->GetXaxis()->SetBinLabel(1, "MFT-MCH-MID (global muon)"); - // hMuonType->GetXaxis()->SetBinLabel(2, "MFT-MCH-MID (global muon other match)"); - // hMuonType->GetXaxis()->SetBinLabel(3, "MFT-MCH"); - // hMuonType->GetXaxis()->SetBinLabel(4, "MCH-MID"); - // hMuonType->GetXaxis()->SetBinLabel(5, "MCH standalone"); - - // fRegistry.add("MFTMCHMID/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{200, 0.0f, 10}}, false); - // fRegistry.add("MFTMCHMID/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {80, -4.f, -2.f}}, false); - // fRegistry.add("MFTMCHMID/hEtaPhi_MatchedMCHMID", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {80, -4.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/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, {{200, 0.0f, 20}}, false); - // fRegistry.add("MFTMCHMID/hChi2MFT", "chi2 MFT;chi2 MFT/ndf", kTH1F, {{200, 0.0f, 20}}, false); - // fRegistry.add("MFTMCHMID/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{200, 0.0f, 20}}, false); - // fRegistry.add("MFTMCHMID/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{200, 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.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); - - fRegistry.add("Pair/uls/gl_gl/hMvsPt", "dimuon;m_{#mu#mu} (GeV/c^{2});p_{T,#mu} (GeV/c);", kTH2F, {{380, 0.2, 4.f}, {100, 0, 10}}, false); - fRegistry.add("Pair/uls/gl_sa/hMvsPt", "dimuon;m_{#mu#mu} (GeV/c^{2});p_{T,#mu} (GeV/c);", kTH2F, {{380, 0.2, 4.f}, {100, 0, 10}}, false); - } - - struct Muon { - int globalIndex = -1; - int collisionId = -1; - int matchMCHTrackId = -1; - int matchMFTTrackId = -1; - uint8_t trackType = 99; - int8_t sign = 0; - float pt = 0; - float eta = 0; - float phi = 0; - float dcaX = 0; // in cm - float dcaY = 0; // in cm - float dcaXY = 0; // in cm - float cXX = 0; - float cYY = 0; - float cXY = 0; - float rAtAbsorberEnd = 0; - float pDCA = 0; - float chi2ndf = 0; - float chi2MatchMCHMID = 0; - - // only for global muons - float ptMatchedMCHMID = 0; - float etaMatchedMCHMID = 0; - float phiMatchedMCHMID = 0; - float chi2MatchMCHMFT = 0; - float chi2mft = -999.f; - uint64_t mftClusterSizesAndTrackFlags = 0; - }; - - bool isSelected(Muon const& muon) - { - if (muon.pt < probeMuonCut.minPt || probeMuonCut.maxPt < muon.pt) { - return false; - } - - if (muon.rAtAbsorberEnd < probeMuonCut.minRabs || probeMuonCut.maxRabs < muon.rAtAbsorberEnd) { - return false; - } - - if (muon.trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack)) { - if (muon.eta < probeMuonCut.minEtaGL || probeMuonCut.maxEtaGL < muon.eta) { - return false; - } - if (probeMuonCut.maxDCAxy < muon.dcaXY) { - return false; - } - if (probeMuonCut.maxMatchingChi2MCHMFT < muon.chi2MatchMCHMFT) { - return false; - } - } else if (muon.trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack)) { - if (muon.eta < probeMuonCut.minEtaSA || probeMuonCut.maxEtaSA < muon.eta) { - return false; - } - if (muon.rAtAbsorberEnd < probeMuonCut.midRabs ? muon.pDCA > probeMuonCut.maxPDCAforSmallR : muon.pDCA > probeMuonCut.maxPDCAforLargeR) { - return false; - } - } else { - return false; - } - - return true; - } - - bool isSelectedTight(Muon const& muon) - { - if (muon.trackType != static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack)) { // tag muon should be tight. - return false; - } - - if (muon.pt < tagMuonCut.minPt || tagMuonCut.maxPt < muon.pt) { - return false; - } - - if (muon.rAtAbsorberEnd < tagMuonCut.minRabs || tagMuonCut.maxRabs < muon.rAtAbsorberEnd) { - return false; - } - - if (muon.chi2ndf < 0.f || tagMuonCut.maxChi2 < muon.chi2ndf) { - return false; - } - - if (muon.eta < tagMuonCut.minEta || tagMuonCut.maxEta < muon.eta) { - return false; - } - - if (tagMuonCut.maxDCAxy < muon.dcaXY) { - return false; - } - - if (tagMuonCut.maxMatchingChi2MCHMFT < muon.chi2MatchMCHMFT) { - return false; - } - - float deta = muon.etaMatchedMCHMID - muon.eta; - float dphi = muon.phiMatchedMCHMID - muon.phi; - // LOGF(info, "muon.trackType = %d, deta = %f, dphi = %f", muon.trackType, deta, dphi); - if (std::sqrt(std::pow(deta / tagMuonCut.maxDEta, 2) + std::pow(dphi / tagMuonCut.maxDPhi, 2)) > 1.f) { - return false; - } - - return true; - } - - template - bool fillMuonInfo(TCollision const& collision, TFwdTrack fwdtrack) - { - if (fwdtrack.trackType() != static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) && fwdtrack.trackType() != static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack)) { - return false; - } - - if (fwdtrack.chi2MatchMCHMID() < 0.f) { // this should never happen. only for protection. - return false; - } - - if (fwdtrack.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); - - float dcaX = propmuonAtPV.getX() - collision.posX(); - float dcaY = propmuonAtPV.getY() - collision.posY(); - float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack - float cXX = propmuonAtPV.getSigma2X(); - float cYY = propmuonAtPV.getSigma2Y(); - float cXY = propmuonAtPV.getSigmaXY(); - - float pDCA = propmuonAtPV.getP() * dcaXY; - int nClustersMFT = 0; - float ptMatchedMCHMID = propmuonAtPV.getPt(); - float etaMatchedMCHMID = propmuonAtPV.getEta(); - float phiMatchedMCHMID = propmuonAtPV.getPhi(); - o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - float chi2mft = -999.f; - uint64_t mftClusterSizesAndTrackFlags = 0; - int ndf_mchmft = 1; - - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - if (fwdtrack.chi2MatchMCHMFT() < 0.f) { - return false; - } // Users have to decide the best match between MFT and MCH-MID at analysis level. The same global muon is repeatedly stored. - - auto mchtrack = fwdtrack.template matchMCHTrack_as(); // MCH-MID - auto mfttrack = fwdtrack.template matchMFTTrack_as(); // MFTsa - if (mfttrack.chi2() < 0.f) { - return false; - } - - if constexpr (isMC) { - if (!mfttrack.has_mcParticle() || !mchtrack.has_mcParticle() || !fwdtrack.has_mcParticle()) { - 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(); - } - - nClustersMFT = mfttrack.nClusters(); - mftClusterSizesAndTrackFlags = mfttrack.mftClusterSizesAndTrackFlags(); - ndf_mchmft = 2.f * (mchtrack.nClusters() + nClustersMFT) - 5.f; - chi2mft = mfttrack.chi2(); - - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); - ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); - etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); - phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); - o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - - 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); - pDCA = mchtrack.p() * dcaXY_Matched; - - if (refitGlobalMuon) { - pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); - } - } else if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToRabs, matchingZ, mBz, mZShift); // this is necessary only for MuonStandaloneTrack - float xAbs = propmuonAtRabs.getX(); - float yAbs = propmuonAtRabs.getY(); - rAtAbsorberEnd = std::sqrt(xAbs * xAbs + yAbs * yAbs); // Redo propagation only for muon tracks // propagation of MFT tracks alredy done in reconstruction - ndf_mchmft = 1; // chi2 is already normalized by ndf for MCH-MID tracks. - - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift); - cXX = propmuonAtDCA.getSigma2X(); - cYY = propmuonAtDCA.getSigma2Y(); - cXY = propmuonAtDCA.getSigmaXY(); - dcaX = propmuonAtDCA.getX() - collision.posX(); - dcaY = propmuonAtDCA.getY() - collision.posY(); - dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - pDCA = fwdtrack.p() * dcaXY; - } else { - return false; - } - - Muon muon; - muon.globalIndex = fwdtrack.globalIndex(); - muon.collisionId = collision.globalIndex(); - muon.trackType = fwdtrack.trackType(); - muon.sign = fwdtrack.sign(); - muon.pt = pt; - muon.eta = eta; - muon.phi = phi; - muon.dcaX = dcaX; - muon.dcaY = dcaY; - muon.dcaXY = dcaXY; - muon.cXX = cXX; - muon.cYY = cYY; - muon.cXY = cXY; - muon.rAtAbsorberEnd = rAtAbsorberEnd; - muon.pDCA = pDCA; - muon.chi2ndf = fwdtrack.chi2() / ndf_mchmft; - muon.ptMatchedMCHMID = ptMatchedMCHMID; - muon.etaMatchedMCHMID = etaMatchedMCHMID; - muon.phiMatchedMCHMID = phiMatchedMCHMID; - muon.chi2mft = chi2mft; - muon.matchMCHTrackId = fwdtrack.matchMCHTrackId(); - muon.matchMFTTrackId = fwdtrack.matchMFTTrackId(); - muon.mftClusterSizesAndTrackFlags = mftClusterSizesAndTrackFlags; - - vecMuons.emplace_back(muon); - return true; - } - - template - bool fillFwdTrackTable(TCollision const& collision, TMuon const& muon, TFwdTrack const& fwdtrack) - { - emprimarymuons(collision.globalIndex(), fwdtrack.globalIndex(), fwdtrack.matchMFTTrackId(), fwdtrack.matchMCHTrackId(), fwdtrack.trackType(), - muon.pt, muon.eta, muon.phi, fwdtrack.sign(), muon.dcaX, muon.dcaY, muon.cXX, muon.cYY, muon.cXY, muon.ptMatchedMCHMID, muon.etaMatchedMCHMID, muon.phiMatchedMCHMID, - fwdtrack.nClusters(), muon.pDCA, muon.rAtAbsorberEnd, fwdtrack.chi2(), fwdtrack.chi2MatchMCHMID(), fwdtrack.chi2MatchMCHMFT(), 9999.f, - fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), muon.mftClusterSizesAndTrackFlags, muon.chi2mft, true, false); - - // if (fillQAHistograms) { - // fRegistry.fill(HIST("hMuonType"), fwdtrack.trackType()); - // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - // 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/hDeltaPt_Pt"), pt, dpt); - // 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"), nClustersMFT); - // fRegistry.fill(HIST("MFTMCHMID/hPDCA_Rabs"), rAtAbsorberEnd, pDCA); - // fRegistry.fill(HIST("MFTMCHMID/hRatAbsorberEnd"), rAtAbsorberEnd); - // fRegistry.fill(HIST("MFTMCHMID/hChi2"), fwdtrack.chi2() / ndf_mchmft); - // fRegistry.fill(HIST("MFTMCHMID/hChi2MFT"), chi2mft / ndf_mft); - // 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/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); - // fRegistry.fill(HIST("MFTMCHMID/hDCAxy"), dcaXY); - // fRegistry.fill(HIST("MFTMCHMID/hDCAxyz"), dcaXY, dcaZ); - // fRegistry.fill(HIST("MFTMCHMID/hDCAxyinSigma"), dcaXYinSigma); - // fRegistry.fill(HIST("MFTMCHMID/hDCAxResolutionvsPt"), pt, std::sqrt(cXX) * 1e+4); // convert cm to um - // fRegistry.fill(HIST("MFTMCHMID/hDCAyResolutionvsPt"), pt, std::sqrt(cYY) * 1e+4); // convert cm to um - // fRegistry.fill(HIST("MFTMCHMID/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um - // 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); - // } else if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - // fRegistry.fill(HIST("MCHMID/hPt"), pt); - // fRegistry.fill(HIST("MCHMID/hEtaPhi"), phi, eta); - // fRegistry.fill(HIST("MCHMID/hEtaPhi_MatchedMCHMID"), phiMatchedMCHMID, etaMatchedMCHMID); - // fRegistry.fill(HIST("MCHMID/hDeltaPt_Pt"), pt, dpt); - // fRegistry.fill(HIST("MCHMID/hDeltaEta_Pt"), pt, deta); - // fRegistry.fill(HIST("MCHMID/hDeltaPhi_Pt"), pt, dphi); - // fRegistry.fill(HIST("MCHMID/hSign"), fwdtrack.sign()); - // fRegistry.fill(HIST("MCHMID/hNclusters"), fwdtrack.nClusters()); - // fRegistry.fill(HIST("MCHMID/hNclustersMFT"), nClustersMFT); - // fRegistry.fill(HIST("MCHMID/hPDCA_Rabs"), rAtAbsorberEnd, pDCA); - // fRegistry.fill(HIST("MCHMID/hRatAbsorberEnd"), rAtAbsorberEnd); - // fRegistry.fill(HIST("MCHMID/hChi2"), fwdtrack.chi2()); - // fRegistry.fill(HIST("MCHMID/hChi2MFT"), chi2mft / ndf_mft); - // fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); - // fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); - // fRegistry.fill(HIST("MCHMID/hDCAxy2D"), dcaX, dcaY); - // fRegistry.fill(HIST("MCHMID/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); - // fRegistry.fill(HIST("MCHMID/hDCAxy"), dcaXY); - // fRegistry.fill(HIST("MCHMID/hDCAxyz"), dcaXY, dcaZ); - // fRegistry.fill(HIST("MCHMID/hDCAxyinSigma"), dcaXYinSigma); - // fRegistry.fill(HIST("MCHMID/hDCAxResolutionvsPt"), pt, std::sqrt(cXX) * 1e+4); // convert cm to um - // fRegistry.fill(HIST("MCHMID/hDCAyResolutionvsPt"), pt, std::sqrt(cYY) * 1e+4); // convert cm to um - // fRegistry.fill(HIST("MCHMID/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um - // } - // } - return true; - } - - SliceCache cache; - Preslice perCollision = o2::aod::fwdtrack::collisionId; - Preslice fwdtrackIndicesPerCollision = aod::track_association::collisionId; - PresliceUnsorted fwdtrackIndicesPerFwdTrack = aod::track_association::fwdtrackId; - PresliceUnsorted fwdtracksPerMCHTrack = aod::fwdtrack::matchMCHTrackId; - - // Filter trackFilter = o2::aod::fwdtrack::trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) || o2::aod::fwdtrack::trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack); - // using filteredMyFwdTracks = soa::Filtered; - // Partition posTracks = o2::aod::fwdtrack::signed1Pt > 0.f; - // Partition negTracks = o2::aod::fwdtrack::signed1Pt < 0.f; - - std::vector vecMuons; - - void processRec(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&) - { - vecMuons.reserve(fwdtracks.size()); - - for (const auto& collision : collisions) { - auto bc = collision.template bc_as(); - initCCDB(bc); - - if (!collision.isSelected()) { - continue; - } - - auto fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); - for (const auto& fwdtrack : fwdtracks_per_coll) { - fillMuonInfo(collision, fwdtrack); - } - - auto pos_muons_per_col = std::views::filter(vecMuons, [](Muon muon) { return muon.sign > 0; }); - auto neg_muons_per_col = std::views::filter(vecMuons, [](Muon muon) { return muon.sign < 0; }); - - // ULS - for (const auto& pos : pos_muons_per_col) { - if (!isSelectedTight(pos)) { // pos is tag, neg is probe - continue; - } - for (const auto& neg : neg_muons_per_col) { - if (!isSelected(neg)) { - continue; - } - ROOT::Math::PtEtaPhiMVector v1(pos.pt, pos.eta, pos.phi, o2::constants::physics::MassMuon); // tag - ROOT::Math::PtEtaPhiMVector v2(neg.pt, neg.eta, neg.phi, o2::constants::physics::MassMuon); // probe - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - if (neg.trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack)) { - if (pos.matchMCHTrackId == neg.matchMCHTrackId || pos.matchMFTTrackId == neg.matchMFTTrackId) { // this should not happen in ULS. only for protection. - continue; - } - fRegistry.fill(HIST("Pair/uls/gl_gl/hMvsPt"), v12.M(), v2.Pt()); - } else if (neg.trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack)) { - if (pos.matchMCHTrackId == neg.globalIndex) { // this should not happen in ULS. only for protection. - continue; - } - fRegistry.fill(HIST("Pair/uls/gl_sa/hMvsPt"), v12.M(), v2.Pt()); - } - if (pairCuts.minMass < v12.M() && v12.M() < pairCuts.maxMass) { - fillFwdTrackTable(collision, neg, fwdtracks.rawIteratorAt(neg.globalIndex)); - } - } // end of neg - } // end of pos - - // ULS - for (const auto& neg : neg_muons_per_col) { - if (!isSelectedTight(neg)) { // neg is tag, pos is probe - continue; - } - for (const auto& pos : pos_muons_per_col) { - if (!isSelected(pos)) { - continue; - } - ROOT::Math::PtEtaPhiMVector v1(neg.pt, neg.eta, neg.phi, o2::constants::physics::MassMuon); // tag - ROOT::Math::PtEtaPhiMVector v2(pos.pt, pos.eta, pos.phi, o2::constants::physics::MassMuon); // probe - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - if (pos.trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack)) { - if (pos.matchMCHTrackId == neg.matchMCHTrackId || pos.matchMFTTrackId == neg.matchMFTTrackId) { // this should not happen in ULS. only for protection. - continue; - } - fRegistry.fill(HIST("Pair/uls/gl_gl/hMvsPt"), v12.M(), v2.Pt()); - } else if (pos.trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack)) { - if (neg.matchMCHTrackId == pos.globalIndex) { // this should not happen in ULS. only for protection. - continue; - } - fRegistry.fill(HIST("Pair/uls/gl_sa/hMvsPt"), v12.M(), v2.Pt()); - } - if (pairCuts.minMass < v12.M() && v12.M() < pairCuts.maxMass) { - fillFwdTrackTable(collision, pos, fwdtracks.rawIteratorAt(pos.globalIndex)); - } - } // end of pos - } // end of neg - - } // end of collision loop - - vecMuons.clear(); - vecMuons.shrink_to_fit(); - } - PROCESS_SWITCH(skimmerPrimaryMuonQC, processRec, "process reconstructed info", false); - - void processRec_SWT(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&) - { - vecMuons.reserve(fwdtracks.size()); - - for (const auto& collision : collisions) { - const auto& bc = collision.template bc_as(); - initCCDB(bc); - - if (!collision.isSelected()) { - continue; - } - - if (collision.triggerMask_raw() == 0) { - continue; - } - - // const auto& fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); - // for (const auto& fwdtrack : fwdtracks_per_coll) { - // if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - // continue; - // } - - // if (!fillFwdTrackTable(collision, fwdtrack, false)) { - // continue; - // } - - // } // end of fwdtrack loop - } // end of collision loop - - vecMuons.clear(); - vecMuons.shrink_to_fit(); - } - PROCESS_SWITCH(skimmerPrimaryMuonQC, processRec_SWT, "process reconstructed info only with standalone", false); - - using filteredMyFwdTracksMC = soa::Filtered; - void processMC(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const&, aod::BCsWithTimestamps const&, aod::McParticles const&) - { - vecMuons.reserve(fwdtracks.size()); - - for (const auto& collision : collisions) { - auto bc = collision.template bc_as(); - initCCDB(bc); - if (!collision.isSelected()) { - continue; - } - if (!collision.has_mcCollision()) { - continue; - } - - // auto fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); - // for (const auto& fwdtrack : fwdtracks_per_coll) { - // if (!fwdtrack.has_mcParticle()) { - // continue; - // } - // if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - // continue; - // } - - // if (!fillFwdTrackTable(collision, fwdtrack, false)) { - // continue; - // } - - // } // end of fwdtrack loop - } // end of collision loop - - vecMuons.clear(); - vecMuons.shrink_to_fit(); - } - PROCESS_SWITCH(skimmerPrimaryMuonQC, processMC, "process reconstructed and MC info", false); - - void processDummy(aod::Collisions const&) {} - PROCESS_SWITCH(skimmerPrimaryMuonQC, processDummy, "process dummy", true); -}; -struct associateAmbiguousMuon { - Produces em_amb_muon_ids; - - SliceCache cache; - PresliceUnsorted perTrack = o2::aod::emprimarymuon::fwdtrackId; - std::vector ambmuon_self_Ids; - - void process(aod::EMPrimaryMuons const& muons) - { - for (const auto& muon : muons) { - auto muons_with_same_trackId = muons.sliceBy(perTrack, muon.fwdtrackId()); - ambmuon_self_Ids.reserve(muons_with_same_trackId.size()); - for (const auto& amb_muon : muons_with_same_trackId) { - if (amb_muon.globalIndex() == muon.globalIndex()) { // don't store myself. - continue; - } - ambmuon_self_Ids.emplace_back(amb_muon.globalIndex()); - } - em_amb_muon_ids(ambmuon_self_Ids); - ambmuon_self_Ids.clear(); - ambmuon_self_Ids.shrink_to_fit(); - } - } -}; - -struct associateSameMuonElement { - Produces glmuon_same_ids; - - SliceCache cache; - PresliceUnsorted perMFTTrack = o2::aod::emprimarymuon::mfttrackId; - PresliceUnsorted perMCHTrack = o2::aod::emprimarymuon::mchtrackId; - std::vector selfIds_per_MFT; - std::vector selfIds_per_MCHMID; - - // Multiple MCH-MID tracks can match with the same MFTsa. This function is to reject such global muons. - void process(aod::EMPrimaryMuons const& muons) - { - for (const auto& muon : muons) { - if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - auto muons_with_same_mfttrackId = muons.sliceBy(perMFTTrack, muon.mfttrackId()); - auto muons_with_same_mchtrackId = muons.sliceBy(perMCHTrack, muon.mchtrackId()); - selfIds_per_MFT.reserve(muons_with_same_mfttrackId.size()); - selfIds_per_MCHMID.reserve(muons_with_same_mchtrackId.size()); - // LOGF(info, "muons_with_same_mchtrackId.size() = %d, muons_with_same_mfttrackId.size() = %d", muons_with_same_mchtrackId.size(), muons_with_same_mfttrackId.size()); - - for (const auto& global_muon : muons_with_same_mfttrackId) { - // LOGF(info, "same MFT: global_muon.globalIndex() = %d, global_muon.mchtrackId() = %d, global_muon.mfttrackId() = %d, global_muon.collisionId() = %d", global_muon.globalIndex(), global_muon.mchtrackId(), global_muon.mfttrackId(), global_muon.collisionId()); - if (global_muon.globalIndex() == muon.globalIndex()) { // don't store myself. - continue; - } - if (global_muon.collisionId() == muon.collisionId()) { // the same global muon is repeatedly stored and associated to different collisions if FTTCA is used. - selfIds_per_MFT.emplace_back(global_muon.globalIndex()); - } - } - - for (const auto& global_muon : muons_with_same_mchtrackId) { - // LOGF(info, "same MCH: global_muon.globalIndex() = %d, global_muon.mchtrackId() = %d, global_muon.mfttrackId() = %d, global_muon.collisionId() = %d", global_muon.globalIndex(), global_muon.mchtrackId(), global_muon.mfttrackId(), global_muon.collisionId()); - if (global_muon.globalIndex() == muon.globalIndex()) { // don't store myself. - continue; - } - if (global_muon.collisionId() == muon.collisionId()) { // the same global muon is repeatedly stored and associated to different collisions if FTTCA is used. - selfIds_per_MCHMID.emplace_back(global_muon.globalIndex()); - } - } - - glmuon_same_ids(selfIds_per_MCHMID, selfIds_per_MFT); - selfIds_per_MFT.clear(); - selfIds_per_MFT.shrink_to_fit(); - selfIds_per_MCHMID.clear(); - selfIds_per_MCHMID.shrink_to_fit(); - } else { - glmuon_same_ids(std::vector{}, std::vector{}); // empty for standalone muons - selfIds_per_MFT.clear(); - selfIds_per_MFT.shrink_to_fit(); - selfIds_per_MCHMID.clear(); - selfIds_per_MCHMID.shrink_to_fit(); - } - } // end of muon loop - } -}; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc, TaskName{"skimmer-primary-muon-qc"}), - adaptAnalysisTask(cfgc, TaskName{"associate-ambiguous-muon"}), - adaptAnalysisTask(cfgc, TaskName{"associate-same-muon-element"})}; -}