diff --git a/Common/Core/fwdtrackUtilities.h b/Common/Core/fwdtrackUtilities.h index b30a0316439..883000d6696 100644 --- a/Common/Core/fwdtrackUtilities.h +++ b/Common/Core/fwdtrackUtilities.h @@ -48,33 +48,6 @@ using SMatrix55 = ROOT::Math::SMatrix; using SMatrix5 = ROOT::Math::SVector; -template -o2::track::TrackParCovFwd getTrackParCovFwd(TFwdTrack const& track, TFwdTrackCov const& cov) -{ - // This function works for (glMuon, glMuon), (saMuon, saMuon) and (MFTTrack, MFTTrackCov). - - double chi2 = track.chi2(); - if constexpr (std::is_same_v, aod::MFTTracksCov::iterator>) { - chi2 = track.chi2(); - } else { - if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - chi2 = track.chi2(); - } else if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - chi2 = track.chi2() * (2.f * track.nClusters() - 5.f); - } - } - - SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); - std::vector v1{cov.cXX(), cov.cXY(), cov.cYY(), cov.cPhiX(), cov.cPhiY(), - cov.cPhiPhi(), cov.cTglX(), cov.cTglY(), cov.cTglPhi(), cov.cTglTgl(), - cov.c1PtX(), cov.c1PtY(), cov.c1PtPhi(), cov.c1PtTgl(), cov.c1Pt21Pt2()}; - SMatrix55 tcovs(v1.begin(), v1.end()); - o2::track::TrackParCovFwd trackparCov{track.z(), tpars, tcovs, chi2}; // this is chi2! Not chi2/ndf. - v1.clear(); - v1.shrink_to_fit(); - return trackparCov; -} - /// Produce TrackParCovFwds for MFT and FwdTracks, w/ or w/o cov, with z shift template o2::track::TrackParCovFwd getTrackParCovFwdShift(TFwdTrack const& track, float zshift, TCovariance const&... covOpt) @@ -107,6 +80,8 @@ o2::track::TrackParCovFwd getTrackParCovFwdShift(TFwdTrack const& track, float z cov.cPhiPhi(), cov.cTglX(), cov.cTglY(), cov.cTglPhi(), cov.cTglTgl(), cov.c1PtX(), cov.c1PtY(), cov.c1PtPhi(), cov.c1PtTgl(), cov.c1Pt21Pt2()}; tcovs = SMatrix55(v1.begin(), v1.end()); + v1.clear(); + v1.shrink_to_fit(); } else { tcovs = SMatrix55{}; } @@ -114,17 +89,46 @@ o2::track::TrackParCovFwd getTrackParCovFwdShift(TFwdTrack const& track, float z return o2::track::TrackParCovFwd(track.z() + zshift, tpars, tcovs, chi2); } +template +o2::track::TrackParCovFwd getTrackParCovFwd(TFwdTrack const& track, TFwdTrackCov const& cov) +{ + return getTrackParCovFwdShift(track, 0, cov); + + // // This function works for (glMuon, glMuon), (saMuon, saMuon) and (MFTTrack, MFTTrackCov). + + // double chi2 = track.chi2(); + // if constexpr (std::is_same_v, aod::MFTTracksCov::iterator>) { + // chi2 = track.chi2(); + // } else { + // if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + // chi2 = track.chi2(); + // } else if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + // chi2 = track.chi2() * (2.f * track.nClusters() - 5.f); + // } + // } + + // SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + // std::vector v1{cov.cXX(), cov.cXY(), cov.cYY(), cov.cPhiX(), cov.cPhiY(), + // cov.cPhiPhi(), cov.cTglX(), cov.cTglY(), cov.cTglPhi(), cov.cTglTgl(), + // cov.c1PtX(), cov.c1PtY(), cov.c1PtPhi(), cov.c1PtTgl(), cov.c1Pt21Pt2()}; + // SMatrix55 tcovs(v1.begin(), v1.end()); + // o2::track::TrackParCovFwd trackparCov{track.z(), tpars, tcovs, chi2}; // this is chi2! Not chi2/ndf. + // v1.clear(); + // v1.shrink_to_fit(); + // return trackparCov; +} + /// 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) +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::track::TrackParCovFwd trackParCovFwd; if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - trackParCovFwd = getTrackParCovFwd(muon, cov); + trackParCovFwd = getTrackParCovFwdShift(muon, zshift, cov); } else if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - trackParCovFwd = getTrackParCovFwd(muon, muon); + trackParCovFwd = getTrackParCovFwdShift(muon, zshift, muon); } else { - trackParCovFwd = getTrackParCovFwd(muon, muon); + trackParCovFwd = getTrackParCovFwdShift(muon, zshift, muon); } o2::dataformats::GlobalFwdTrack propmuon = propagateTrackParCovFwd(trackParCovFwd, muon.trackType(), collision, endPoint, matchingZ, bzkG); diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index 018abc465c4..b595c111543 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -71,13 +71,13 @@ struct skimmerPrimaryMuon { 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"}; - Configurable minPt{"minPt", 0.1, "min pt for muon"}; + 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 MCH-MID"}; Configurable maxEtaSA{"maxEtaSA", -2.5, "max. eta acceptance for MCH-MID"}; 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 minRabsGL{"minRabsGL", 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. + Configurable minRabsGL{"minRabsGL", 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. Configurable minRabs{"minRabs", 17.6, "min. R at absorber end"}; Configurable midRabs{"midRabs", 26.5, "middle R at absorber end for pDCA cut"}; Configurable maxRabs{"maxRabs", 89.5, "max. R at absorber end"}; @@ -93,10 +93,16 @@ struct skimmerPrimaryMuon { Configurable maxDEta{"maxDEta", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; Configurable maxDPhi{"maxDPhi", 1e+10f, "max. dphi between MFT-MCH-MID and MCH-MID"}; + // 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"}; + 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/"}; @@ -114,6 +120,7 @@ struct skimmerPrimaryMuon { } mRunNumber = 0; mBz = 0; + mZShift = 0; } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -136,6 +143,20 @@ struct skimmerPrimaryMuon { 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() @@ -238,7 +259,7 @@ struct skimmerPrimaryMuon { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); + 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(); @@ -314,13 +335,13 @@ struct skimmerPrimaryMuon { return false; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); + 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); + 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); @@ -328,9 +349,9 @@ struct skimmerPrimaryMuon { if constexpr (withMFTCov) { auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz); // propagated to matching plane - o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane + 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(); @@ -338,7 +359,7 @@ struct skimmerPrimaryMuon { o2::math_utils::bringTo02Pi(phiMatchedMCHMIDatMP); o2::math_utils::bringTo02Pi(phiMatchedMFTatMP); - o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + 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 globalMuon = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToVertex, matchingZ, mBz); pt = globalMuon.getPt(); @@ -367,12 +388,12 @@ struct skimmerPrimaryMuon { 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); // this is necessary only for 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 - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); + o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift); cXX = propmuonAtDCA.getSigma2X(); cYY = propmuonAtDCA.getSigma2Y(); cXY = propmuonAtDCA.getSigmaXY(); diff --git a/PWGEM/Dilepton/Tasks/matchingMFT.cxx b/PWGEM/Dilepton/Tasks/matchingMFT.cxx index 8dbf1b1b136..f5183fb0c12 100644 --- a/PWGEM/Dilepton/Tasks/matchingMFT.cxx +++ b/PWGEM/Dilepton/Tasks/matchingMFT.cxx @@ -82,6 +82,11 @@ struct matchingMFT { Configurable minNclustersMFT{"minNclustersMFT", 5, "min nclusters MFT"}; Configurable refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"}; + // for z shift for propagation + Configurable cfgApplyZShiftFromCCDB{"cfgApplyZShiftFromCCDB", false, "flag to apply z shift from CCDB"}; + 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 requireTrueAssociation{"requireTrueAssociation", false, "flag to require true mc collision association"}; Configurable maxRelDPt{"maxRelDPt", 1e+10f, "max. relative dpt between MFT-MCH-MID and MCH-MID"}; Configurable maxDEta{"maxDEta", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; @@ -130,6 +135,7 @@ struct matchingMFT { Service ccdb; int mRunNumber = -1; float mBz = 0; + float mZShift = 0; template void initCCDB(TBC const& bc) @@ -152,6 +158,20 @@ struct matchingMFT { 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() @@ -358,11 +378,11 @@ struct matchingMFT { return; // do nothing } - auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz); // propagated to matching plane + auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz, mZShift); // propagated to matching plane auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane + o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwdShift(mfttrack, mZShift, mfttrackcov); // values at innermost update + mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane dx = muonAtMP.getX() - mftsaAtMP.getX(); dy = muonAtMP.getY() - mftsaAtMP.getY(); // o2::math_utils::bringToPMPi(dphi); @@ -430,9 +450,9 @@ struct matchingMFT { bool isPrimary = mcParticle_MCHMID.isPhysicalPrimary() || mcParticle_MCHMID.producedByGenerator(); bool isMatched = (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) && (mcParticle_MFT.mcCollisionId() == mcParticle_MCHMID.mcCollisionId()); - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); + o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz, mZShift); float pt = propmuonAtPV.getPt(); float eta = propmuonAtPV.getEta(); @@ -472,7 +492,7 @@ struct matchingMFT { if constexpr (withMFTCov) { auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + 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 globalMuonAtPV = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToVertex, matchingZ, mBz);