diff --git a/Common/Core/fwdtrackUtilities.h b/Common/Core/fwdtrackUtilities.h index a07f461394d..14bc744e634 100644 --- a/Common/Core/fwdtrackUtilities.h +++ b/Common/Core/fwdtrackUtilities.h @@ -18,17 +18,17 @@ #ifndef COMMON_CORE_FWDTRACKUTILITIES_H_ #define COMMON_CORE_FWDTRACKUTILITIES_H_ +#include "DetectorsBase/GeometryManager.h" +#include "Field/MagneticField.h" #include "Framework/AnalysisDataModel.h" -#include -#include -#include -#include -#include -#include +#include "GlobalTracking/MatchGlobalFwd.h" +#include "MCHTracking/TrackExtrap.h" +#include "ReconstructionDataFormats/GlobalFwdTrack.h" +#include "ReconstructionDataFormats/TrackFwd.h" -#include -#include -#include +#include "Math/MatrixRepresentationsStatic.h" +#include "Math/SMatrix.h" +#include "TGeoGlobalMagField.h" #include #include @@ -51,8 +51,7 @@ using SMatrix5 = ROOT::Math::SVector; template o2::track::TrackParCovFwd getTrackParCovFwd(TFwdTrack const& track, TFwdTrackCov const& cov) { - // This function works for both (saMuon, saMuon) and (MFTTrack, MFTTrackCov). - // Don't use covariant matrix of global muons stored in AO2D.root. + // This function works for (glMuon, glMuon), (saMuon, saMuon) and (MFTTrack, MFTTrackCov). double chi2 = track.chi2(); if constexpr (std::is_same_v, aod::MFTTracksCov::iterator>) { @@ -128,11 +127,12 @@ o2::dataformats::GlobalFwdTrack propagateTrackParCovFwd(TFwdTrackParCov const& f // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); // auto Bz = field->getBz(centerMFT); // Get field at centre of MFT in kG. - auto geoMan = o2::base::GeometryManager::meanMaterialBudget(fwdtrack.getX(), fwdtrack.getY(), fwdtrack.getZ(), collision.posX(), collision.posY(), collision.posZ()); - auto x2x0 = static_cast(geoMan.meanX2X0); - if (endPoint == propagationPoint::kToVertex) { - fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, bzkG, x2x0); + // auto geoMan = o2::base::GeometryManager::meanMaterialBudget(fwdtrack.getX(), fwdtrack.getY(), fwdtrack.getZ(), collision.posX(), collision.posY(), collision.posZ()); + // auto x2x0 = static_cast(geoMan.meanX2X0); + // fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, bzkG, x2x0); + std::array dcaInfOrig{999.f, 999.f, 999.f}; + fwdtrack.propagateToDCAhelix(bzkG, {collision.posX(), collision.posY(), collision.posZ()}, dcaInfOrig); } else if (endPoint == propagationPoint::kToDCA) { fwdtrack.propagateToZhelix(collision.posZ(), bzkG); } else if (endPoint == propagationPoint::kToMatchingPlane) { diff --git a/PWGEM/Dilepton/Core/Dilepton.h b/PWGEM/Dilepton/Core/Dilepton.h index 1cfe6bf0efb..fd01163f794 100644 --- a/PWGEM/Dilepton/Core/Dilepton.h +++ b/PWGEM/Dilepton/Core/Dilepton.h @@ -418,111 +418,84 @@ struct Dilepton { emh_pos = new TEMH(ndepth); emh_neg = new TEMH(ndepth); - emh_pair_uls = new MyEMH_pair(ndepth); - emh_pair_lspp = new MyEMH_pair(ndepth); - emh_pair_lsmm = new MyEMH_pair(ndepth); - if (accBins.ConfMllAccBins.value[0] == VARIABLE_WIDTH) { - mll_bin_edges = std::vector(accBins.ConfMllAccBins.value.begin(), accBins.ConfMllAccBins.value.end()); - mll_bin_edges.erase(mll_bin_edges.begin()); - for (const auto& edge : mll_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: mll_bin_edges = %f", edge); - } - } else { - int nbins = static_cast(accBins.ConfMllAccBins.value[0]); - float xmin = static_cast(accBins.ConfMllAccBins.value[1]); - float xmax = static_cast(accBins.ConfMllAccBins.value[2]); - mll_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - mll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: mll_bin_edges[%d] = %f", i, mll_bin_edges[i]); - } - } + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kPolarization)) { // only for polarization - if (accBins.ConfPtllAccBins.value[0] == VARIABLE_WIDTH) { - ptll_bin_edges = std::vector(accBins.ConfPtllAccBins.value.begin(), accBins.ConfPtllAccBins.value.end()); - ptll_bin_edges.erase(ptll_bin_edges.begin()); - for (const auto& edge : ptll_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: ptll_bin_edges = %f", edge); - } - } else { - int nbins = static_cast(accBins.ConfPtllAccBins.value[0]); - float xmin = static_cast(accBins.ConfPtllAccBins.value[1]); - float xmax = static_cast(accBins.ConfPtllAccBins.value[2]); - ptll_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - ptll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: ptll_bin_edges[%d] = %f", i, ptll_bin_edges[i]); - } - } + emh_pair_uls = new MyEMH_pair(ndepth); + emh_pair_lspp = new MyEMH_pair(ndepth); + emh_pair_lsmm = new MyEMH_pair(ndepth); - if (accBins.ConfEtallAccBins.value[0] == VARIABLE_WIDTH) { - etall_bin_edges = std::vector(accBins.ConfEtallAccBins.value.begin(), accBins.ConfEtallAccBins.value.end()); - etall_bin_edges.erase(etall_bin_edges.begin()); - for (const auto& edge : etall_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: etall_bin_edges = %f", edge); + if (accBins.ConfMllAccBins.value[0] == VARIABLE_WIDTH) { + mll_bin_edges = std::vector(accBins.ConfMllAccBins.value.begin(), accBins.ConfMllAccBins.value.end()); + mll_bin_edges.erase(mll_bin_edges.begin()); + for (const auto& edge : mll_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: mll_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfMllAccBins.value[0]); + float xmin = static_cast(accBins.ConfMllAccBins.value[1]); + float xmax = static_cast(accBins.ConfMllAccBins.value[2]); + mll_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + mll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: mll_bin_edges[%d] = %f", i, mll_bin_edges[i]); + } } - } else { - int nbins = static_cast(accBins.ConfEtallAccBins.value[0]); - float xmin = static_cast(accBins.ConfEtallAccBins.value[1]); - float xmax = static_cast(accBins.ConfEtallAccBins.value[2]); - etall_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - etall_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: etall_bin_edges[%d] = %f", i, etall_bin_edges[i]); + + if (accBins.ConfPtllAccBins.value[0] == VARIABLE_WIDTH) { + ptll_bin_edges = std::vector(accBins.ConfPtllAccBins.value.begin(), accBins.ConfPtllAccBins.value.end()); + ptll_bin_edges.erase(ptll_bin_edges.begin()); + for (const auto& edge : ptll_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: ptll_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfPtllAccBins.value[0]); + float xmin = static_cast(accBins.ConfPtllAccBins.value[1]); + float xmax = static_cast(accBins.ConfPtllAccBins.value[2]); + ptll_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + ptll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: ptll_bin_edges[%d] = %f", i, ptll_bin_edges[i]); + } } - } - if (accBins.ConfPhillAccBins.value[0] == VARIABLE_WIDTH) { - phill_bin_edges = std::vector(accBins.ConfPhillAccBins.value.begin(), accBins.ConfPhillAccBins.value.end()); - phill_bin_edges.erase(phill_bin_edges.begin()); - for (const auto& edge : phill_bin_edges) { - LOGF(info, "VARIABLE_WIDTH: phill_bin_edges = %f", edge); + if (accBins.ConfEtallAccBins.value[0] == VARIABLE_WIDTH) { + etall_bin_edges = std::vector(accBins.ConfEtallAccBins.value.begin(), accBins.ConfEtallAccBins.value.end()); + etall_bin_edges.erase(etall_bin_edges.begin()); + for (const auto& edge : etall_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: etall_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfEtallAccBins.value[0]); + float xmin = static_cast(accBins.ConfEtallAccBins.value[1]); + float xmax = static_cast(accBins.ConfEtallAccBins.value[2]); + etall_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + etall_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: etall_bin_edges[%d] = %f", i, etall_bin_edges[i]); + } } - } else { - int nbins = static_cast(accBins.ConfPhillAccBins.value[0]); - float xmin = static_cast(accBins.ConfPhillAccBins.value[1]); - float xmax = static_cast(accBins.ConfPhillAccBins.value[2]); - phill_bin_edges.resize(nbins + 1); - for (int i = 0; i < nbins + 1; i++) { - phill_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; - LOGF(info, "FIXED_WIDTH: phill_bin_edges[%d] = %f", i, phill_bin_edges[i]); + + if (accBins.ConfPhillAccBins.value[0] == VARIABLE_WIDTH) { + phill_bin_edges = std::vector(accBins.ConfPhillAccBins.value.begin(), accBins.ConfPhillAccBins.value.end()); + phill_bin_edges.erase(phill_bin_edges.begin()); + for (const auto& edge : phill_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: phill_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfPhillAccBins.value[0]); + float xmin = static_cast(accBins.ConfPhillAccBins.value[1]); + float xmax = static_cast(accBins.ConfPhillAccBins.value[2]); + phill_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + phill_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: phill_bin_edges[%d] = %f", i, phill_bin_edges[i]); + } } - } - int nM = mll_bin_edges.size(); - int nPt = ptll_bin_edges.size(); - int nEta = etall_bin_edges.size(); - int nPhi = phill_bin_edges.size(); - - // emhs_pair_uls.resize(nM); - // emhs_pair_lspp.resize(nM); - // emhs_pair_lsmm.resize(nM); - // for (int im = 0;im>>> emhs_pair_uls; // 4D{m, pt, eta, phi} - // std::vector>>> emhs_pair_lspp; // 4D{m, pt, eta, phi} - // std::vector>>> emhs_pair_lsmm; // 4D{m, pt, eta, phi} std::map, uint64_t> map_mixed_eventId_to_globalBC; std::vector used_trackIds_per_col; @@ -1730,19 +1676,6 @@ struct Dilepton { emh_pair_uls->AddCollisionIdAtLast(key_bin, key_df_collision); emh_pair_lspp->AddCollisionIdAtLast(key_bin, key_df_collision); emh_pair_lsmm->AddCollisionIdAtLast(key_bin, key_df_collision); - - // for (int im = 0;imAddCollisionIdAtLast(key_bin, key_df_collision); - // emhs_pair_lspp[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision); - // emhs_pair_lsmm[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision); - // } - // } - // } - // } - } // end of if pair exist } // end of collision loop diff --git a/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx b/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx index c5a23f31ffe..1e4e923be1e 100644 --- a/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx +++ b/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx @@ -181,8 +181,8 @@ struct CreateEMEventDilepton { q3xft0m = collision.qvecFT0MReVec()[1], q3xft0a = collision.qvecFT0AReVec()[1], q3xft0c = collision.qvecFT0CReVec()[1], q3xfv0a = collision.qvecFV0AReVec()[1], q3xbpos = collision.qvecBPosReVec()[1], q3xbneg = collision.qvecBNegReVec()[1], q3xbtot = collision.qvecBTotReVec()[1]; q3yft0m = collision.qvecFT0MImVec()[1], q3yft0a = collision.qvecFT0AImVec()[1], q3yft0c = collision.qvecFT0CImVec()[1], q3yfv0a = collision.qvecFV0AImVec()[1], q3ybpos = collision.qvecBPosImVec()[1], q3ybneg = collision.qvecBNegImVec()[1], q3ybtot = collision.qvecBTotImVec()[1]; } else if (collision.qvecFT0CReVec().size() >= 1) { // harmonics 2 - q2xft0m = collision.qvecFT0MReVec()[0], q2xft0a = collision.qvecFT0AReVec()[0], q2xft0c = collision.qvecFT0CReVec()[0], q2xbpos = collision.qvecBPosReVec()[0], q2xbneg = collision.qvecBNegReVec()[0], q2xbtot = collision.qvecBTotReVec()[0]; - q2yft0m = collision.qvecFT0MImVec()[0], q2yft0a = collision.qvecFT0AImVec()[0], q2yft0c = collision.qvecFT0CImVec()[0], q2ybpos = collision.qvecBPosImVec()[0], q2ybneg = collision.qvecBNegImVec()[0], q2ybtot = collision.qvecBTotImVec()[0]; + q2xft0m = collision.qvecFT0MReVec()[0], q2xft0a = collision.qvecFT0AReVec()[0], q2xft0c = collision.qvecFT0CReVec()[0], q2xfv0a = collision.qvecFV0AReVec()[0], q2xbpos = collision.qvecBPosReVec()[0], q2xbneg = collision.qvecBNegReVec()[0], q2xbtot = collision.qvecBTotReVec()[0]; + q2yft0m = collision.qvecFT0MImVec()[0], q2yft0a = collision.qvecFT0AImVec()[0], q2yft0c = collision.qvecFT0CImVec()[0], q2yfv0a = collision.qvecFV0AImVec()[0], q2ybpos = collision.qvecBPosImVec()[0], q2ybneg = collision.qvecBNegImVec()[0], q2ybtot = collision.qvecBTotImVec()[0]; } event_qvec( q2xft0m, q2yft0m, q2xft0a, q2yft0a, q2xft0c, q2yft0c, q2xfv0a, q2yfv0a, q2xbpos, q2ybpos, q2xbneg, q2ybneg, q2xbtot, q2ybtot, diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index 353fbc37fdd..d03e0907811 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -164,6 +164,7 @@ struct skimmerPrimaryMuon { 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.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); @@ -234,26 +235,25 @@ struct skimmerPrimaryMuon { float phi = propmuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phi); - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); - float cXXatDCA = propmuonAtDCA.getSigma2X(); - float cYYatDCA = propmuonAtDCA.getSigma2Y(); - float cXYatDCA = propmuonAtDCA.getSigmaXY(); - - float dcaX = propmuonAtDCA.getX() - collision.posX(); - float dcaY = propmuonAtDCA.getY() - collision.posY(); + 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); float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack + float cXX = propmuonAtPV.getSigma2X(); + float cYY = propmuonAtPV.getSigma2Y(); + float cXY = propmuonAtPV.getSigmaXY(); - float det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat + float det = cXX * cYY - cXY * cXY; // determinanat 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 + dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYY + dcaY * dcaY * cXX - 2.f * dcaX * dcaY * cXY) / det / 2.f)); // dca xy in sigma } float sigma_dcaXY = dcaXY / dcaXYinSigma; - float pDCA = fwdtrack.p() * dcaXY; + float pDCA = propmuonAtPV.getP() * dcaXY; int nClustersMFT = 0; float ptMatchedMCHMID = propmuonAtPV.getPt(); float etaMatchedMCHMID = propmuonAtPV.getEta(); @@ -317,51 +317,70 @@ struct skimmerPrimaryMuon { float dcaXY_Matched = std::sqrt(dcaX_Matched * dcaX_Matched + dcaY_Matched * dcaY_Matched); pDCA = mchtrack.p() * dcaXY_Matched; - if (refitGlobalMuon) { - // eta = mfttrack.eta(); - // phi = mfttrack.phi(); - // o2::math_utils::bringTo02Pi(phi); - eta = propmuonAtDCA.getEta(); - phi = propmuonAtDCA.getPhi(); + 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 + 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 = getTrackParCovFwd(mfttrack, 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(); + phi = globalMuon.getPhi(); o2::math_utils::bringTo02Pi(phi); - pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); - // x = mfttrack.x(); - // y = mfttrack.y(); - // z = mfttrack.z(); - // tgl = mfttrack.tgl(); - - 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 - 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 = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. - // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); - // pt = globalMuonAtDCA.getPt(); - // eta = globalMuonAtDCA.getEta(); - // phi = globalMuonAtDCA.getPhi(); - // o2::math_utils::bringTo02Pi(phi); - // cXXatDCA = globalMuonAtDCA.getSigma2X(); - // cYYatDCA = globalMuonAtDCA.getSigma2Y(); - // cXYatDCA = globalMuonAtDCA.getSigmaXY(); - // dcaX = globalMuonAtDCA.getX() - collision.posX(); - // dcaY = globalMuonAtDCA.getY() - collision.posY(); + cXX = globalMuon.getSigma2X(); + cYY = globalMuon.getSigma2Y(); + cXY = globalMuon.getSigmaXY(); + dcaX = globalMuon.getX() - collision.posX(); + dcaY = globalMuon.getY() - collision.posY(); + dcaZ = globalMuon.getZ() - collision.posZ(); + dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + det = cXX * cYY - cXY * cXY; // determinanat + dcaXYinSigma = 999.f; + if (det < 0) { + dcaXYinSigma = 999.f; + } else { + dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYY + dcaY * dcaY * cXX - 2.f * dcaX * dcaY * cXY) / det / 2.f)); // dca xy in sigma } + sigma_dcaXY = dcaXY / dcaXYinSigma; + } + + 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); // 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); + cXX = propmuonAtDCA.getSigma2X(); + cYY = propmuonAtDCA.getSigma2Y(); + cXY = propmuonAtDCA.getSigmaXY(); + dcaX = propmuonAtDCA.getX() - collision.posX(); + dcaY = propmuonAtDCA.getY() - collision.posY(); + dcaZ = propmuonAtDCA.getZ() - collision.posZ(); + dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + pDCA = fwdtrack.p() * dcaXY; + + det = cXX * cYY - cXY * cXY; // determinanat + dcaXYinSigma = 999.f; + if (det < 0) { + dcaXYinSigma = 999.f; + } else { + dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYY + dcaY * dcaY * cXX - 2.f * dcaX * dcaY * cXY) / det / 2.f)); // dca xy in sigma + } + sigma_dcaXY = dcaXY / dcaXYinSigma; } else { return; } @@ -383,12 +402,12 @@ struct skimmerPrimaryMuon { // LOGF(info, "isAmbiguous = %d, isAssociatedToMPC = %d, fwdtrack.globalIndex() = %d, fwdtrack.collisionId() = %d, collision.globalIndex() = %d", isAmbiguous, isAssociatedToMPC, fwdtrack.globalIndex(), fwdtrack.collisionId(), collision.globalIndex()); emprimarymuons(collision.globalIndex(), fwdtrack.globalIndex(), fwdtrack.matchMFTTrackId(), fwdtrack.matchMCHTrackId(), fwdtrack.trackType(), - pt, eta, phi, fwdtrack.sign(), dcaX, dcaY, cXXatDCA, cYYatDCA, cXYatDCA, ptMatchedMCHMID, etaMatchedMCHMID, phiMatchedMCHMID, + pt, eta, phi, fwdtrack.sign(), dcaX, dcaY, cXX, cYY, cXY, ptMatchedMCHMID, etaMatchedMCHMID, phiMatchedMCHMID, etaMatchedMCHMIDatMP, phiMatchedMCHMIDatMP, etaMatchedMFTatMP, phiMatchedMFTatMP, fwdtrack.nClusters(), pDCA, rAtAbsorberEnd, fwdtrack.chi2(), fwdtrack.chi2MatchMCHMID(), fwdtrack.chi2MatchMCHMFT(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mftClusterSizesAndTrackFlags, chi2mft, isAssociatedToMPC, isAmbiguous); - const auto& fwdcov = propmuonAtPV.getCovariances(); // covatiant matrix at PV + const auto& fwdcov = propmuonAtPV.getCovariances(); // covatiance matrix at PV emprimarymuonscov( fwdcov(0, 0), fwdcov(0, 1), fwdcov(1, 1), @@ -425,12 +444,13 @@ struct skimmerPrimaryMuon { 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(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); + 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(cXXatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um + 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 } else if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { fRegistry.fill(HIST("MCHMID/hPt"), pt); fRegistry.fill(HIST("MCHMID/hEtaPhi"), phi, eta); @@ -450,12 +470,13 @@ struct skimmerPrimaryMuon { 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(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); + 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(cXXatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MCHMID/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MCHMID/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um + 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 } } } diff --git a/PWGEM/Dilepton/Tasks/createResolutionMap.cxx b/PWGEM/Dilepton/Tasks/createResolutionMap.cxx index a887a713263..5df5855d171 100644 --- a/PWGEM/Dilepton/Tasks/createResolutionMap.cxx +++ b/PWGEM/Dilepton/Tasks/createResolutionMap.cxx @@ -557,19 +557,18 @@ struct CreateResolutionMap { return; } o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(muon, muon, collision, propagationPoint::kToVertex, muoncuts.matchingZ, mBzMFT); - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(muon, muon, collision, propagationPoint::kToDCA, muoncuts.matchingZ, mBzMFT); float pt = propmuonAtPV.getPt(); float eta = propmuonAtPV.getEta(); float phi = propmuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phi); - float dcaX = propmuonAtDCA.getX() - collision.posX(); - float dcaY = propmuonAtDCA.getY() - collision.posY(); + float dcaX = propmuonAtPV.getX() - collision.posX(); + float dcaY = propmuonAtPV.getY() - collision.posY(); float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); float rAtAbsorberEnd = muon.rAtAbsorberEnd(); // this works only for GlobalMuonTrack - float pDCA = muon.p() * dcaXY; + float pDCA = propmuonAtPV.getP() * dcaXY; int nClustersMFT = 0; float ptMatchedMCHMID = propmuonAtPV.getPt(); float etaMatchedMCHMID = propmuonAtPV.getEta(); @@ -602,11 +601,14 @@ struct CreateResolutionMap { etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, muoncuts.matchingZ, mBzMFT); - 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, muoncuts.matchingZ, mBzMFT); + // 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; + pDCA = propmuonAtPV.getP() * dcaXY; + nClustersMFT = mfttrack.nClusters(); float chi2mft = mfttrack.chi2() / (2.f * nClustersMFT - 5.f); @@ -622,40 +624,32 @@ struct CreateResolutionMap { return; } - if (muoncuts.refitGlobalMuon) { - // eta = mfttrack.eta(); - // phi = mfttrack.phi(); - // o2::math_utils::bringTo02Pi(phi); - eta = propmuonAtDCA.getEta(); - phi = propmuonAtDCA.getPhi(); + if constexpr (withMFTCov) { + auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, muoncuts.matchingZ, mBzMFT); // propagated to matching plane + o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + mftsaAtMP.propagateToZhelix(muoncuts.matchingZ, mBzMFT); // 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 = getTrackParCovFwd(mfttrack, 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, muon.trackType(), collision, propagationPoint::kToVertex, muoncuts.matchingZ, mBzMFT); + pt = globalMuonAtPV.getPt(); + eta = globalMuonAtPV.getEta(); + phi = globalMuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phi); - pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); + dcaX = globalMuonAtPV.getX() - collision.posX(); + dcaY = globalMuonAtPV.getY() - collision.posY(); + dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + } - if constexpr (withMFTCov) { - auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, muoncuts.matchingZ, mBzMFT); // propagated to matching plane - o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - mftsaAtMP.propagateToZhelix(muoncuts.matchingZ, mBzMFT); // 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 = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. - // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); - // pt = globalMuonAtDCA.getPt(); - // eta = globalMuonAtDCA.getEta(); - // phi = globalMuonAtDCA.getPhi(); - // o2::math_utils::bringTo02Pi(phi); - // cXXatDCA = globalMuonAtDCA.getSigma2X(); - // cYYatDCA = globalMuonAtDCA.getSigma2Y(); - // cXYatDCA = globalMuonAtDCA.getSigmaXY(); - // dcaX = globalMuonAtDCA.getX() - collision.posX(); - // dcaY = globalMuonAtDCA.getY() - collision.posY(); - } + if (muoncuts.refitGlobalMuon) { + pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); } float dpt = (ptMatchedMCHMID - pt) / pt; @@ -690,6 +684,12 @@ struct CreateResolutionMap { 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(muon, muon, collision, propagationPoint::kToDCA, muoncuts.matchingZ, mBzMFT); + dcaX = propmuonAtDCA.getX() - collision.posX(); + dcaY = propmuonAtDCA.getY() - collision.posY(); + dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + pDCA = muon.p() * dcaXY; } else { return; } diff --git a/PWGEM/Dilepton/Tasks/matchingMFT.cxx b/PWGEM/Dilepton/Tasks/matchingMFT.cxx index 154d26aded1..9b0d1d6d2d2 100644 --- a/PWGEM/Dilepton/Tasks/matchingMFT.cxx +++ b/PWGEM/Dilepton/Tasks/matchingMFT.cxx @@ -77,7 +77,6 @@ struct matchingMFT { Configurable maxChi2MFT{"maxChi2MFT", 1e+6f, "max. chi2/ndf for MFT track in global muon"}; Configurable minNclustersMFT{"minNclustersMFT", 5, "min nclusters MFT"}; Configurable refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"}; - Configurable propagateToDCAhelix{"propagateToDCAhelix", false, "flag to use propagateToDCAhelix"}; 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"}; @@ -113,10 +112,6 @@ struct matchingMFT { LOGF(fatal, "Cannot enable doprocessWithoutFTTCA and doprocessWithFTTCA at the same time. Please choose one."); } - if (refitGlobalMuon && propagateToDCAhelix) { - LOGF(fatal, "Cannot enable doprocessWithoutFTTCA and doprocessWithFTTCA at the same time. Please choose one."); - } - ccdb->setURL(ccdburl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -191,7 +186,7 @@ struct matchingMFT { fRegistry.add("MFTMCHMID/primary/correct/hNclustersMFT", "NclustersMFT;Nclusters MFT", kTH1F, {{11, -0.5f, 10.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hMFTClusterMap", "MFT cluster map", kTH1F, {{1024, -0.5, 1023.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hRatAbsorberEnd", "R at absorber end;R at absorber end (cm)", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("MFTMCHMID/primary/correct/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/primary/correct/hPDCA_Rabs", "pDCA vs. Rabs;R at absorber end (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hChi2", "chi2;chi2/ndf", kTH1F, {{100, 0.0f, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hChi2MFT", "chi2 MFT/ndf;chi2 MFT/ndf", kTH1F, {{100, 0.0f, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); @@ -216,7 +211,6 @@ struct matchingMFT { fRegistry.add("MFTMCHMID/primary/correct/hDeltaPhi_Neg", "#varphi resolution;p_{T}^{gen} (GeV/c);#varphi^{rec} - #varphi^{gen} (rad.)", kTH2F, {{100, 0, 10}, {400, -0.2, +0.2}}, false); fRegistry.addClone("MFTMCHMID/primary/correct/", "MFTMCHMID/primary/wrong/"); fRegistry.addClone("MFTMCHMID/primary/", "MFTMCHMID/secondary/"); - // fRegistry.addClone("MFTMCHMID/", "MCHMID/"); } bool isSelected(const float pt, const float eta, const float rAtAbsorberEnd, const float pDCA, const float chi2_per_ndf, const uint8_t trackType, const float dcaXY) @@ -331,15 +325,6 @@ struct matchingMFT { deta = muonAtMP.getEta() - mftsaAtMP.getEta(); dphi = muonAtMP.getPhi() - mftsaAtMP.getPhi(); o2::math_utils::bringToPMPi(dphi); - // reldpt = (muonAtMP.getPt() - mftsaAtMP.getPt()) / muonAtMP.getPt(); - - // o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV, mftsa); // this is track at IU. - // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToMatchingPlane, matchingZ); - // deta = propmuonAtPV.getEta() - globalMuonAtDCA.getEta(); - // dphi = propmuonAtPV.getPhi() - globalMuonAtDCA.getPhi(); - // o2::math_utils::bringToPMPi(dphi); - // reldpt = (globalMuonAtDCA.getPt() - propmuonAtPV.getPt()) / propmuonAtPV.getPt(); } template @@ -396,31 +381,30 @@ struct matchingMFT { 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 propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, 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); float pt = propmuonAtPV.getPt(); float eta = propmuonAtPV.getEta(); float phi = propmuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phi); - float cXXatDCA = propmuonAtDCA.getSigma2X(); - float cYYatDCA = propmuonAtDCA.getSigma2Y(); - float cXYatDCA = propmuonAtDCA.getSigmaXY(); + float cXX = propmuonAtPV.getSigma2X(); + float cYY = propmuonAtPV.getSigma2Y(); + float cXY = propmuonAtPV.getSigmaXY(); - float dcaX = propmuonAtDCA.getX() - collision.posX(); - float dcaY = propmuonAtDCA.getY() - collision.posY(); - float dcaZ = propmuonAtDCA.getZ() - collision.posZ(); // 0 at this point. - float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack + float dcaX = propmuonAtPV.getX() - collision.posX(); + float dcaY = propmuonAtPV.getY() - collision.posY(); + float dcaZ = propmuonAtPV.getZ() - collision.posZ(); // 0 at this point. + float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - float det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat + float det = cXX * cYY - cXY * cXY; // determinanat float dcaXYinSigma = 999.f; if (det < 0) { dcaXYinSigma = 999.f; } else { - dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2. * dcaX * dcaY * cXYatDCA) / det / 2.)); // dca xy in sigma + dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYY + dcaY * dcaY * cXX - 2. * dcaX * dcaY * cXY) / det / 2.)); // dca xy in sigma } float sigma_dcaXY = dcaXY / dcaXYinSigma; @@ -428,94 +412,46 @@ struct matchingMFT { 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 pDCA = propmuonAtPV.getP() * dcaXY; float dphiMP = 999.f, detaMP = 999.f; - if (refitGlobalMuon) { - // eta = mfttrack.eta(); - // phi = mfttrack.phi(); - // o2::math_utils::bringTo02Pi(phi); - eta = propmuonAtDCA.getEta(); - phi = propmuonAtDCA.getPhi(); - o2::math_utils::bringTo02Pi(phi); - pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); - - if constexpr (withMFTCov) { - // auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - // o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update - // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. - // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); - - // eta = globalMuonAtDCA.getEta(); - // phi = globalMuonAtDCA.getPhi(); - // o2::math_utils::bringTo02Pi(phi); - // pt = globalMuonAtDCA.getPt(); - // p = globalMuonAtDCA.getP(); - - // eta = globalMuonRefit.getEta(); - // phi = globalMuonRefit.getPhi(); - // o2::math_utils::bringTo02Pi(phi); - // pt = globalMuonRefit.getPt(); - // p = globalMuonRefit.getP(); - - // cXXatDCA = globalMuonAtDCA.getSigma2X(); - // cYYatDCA = globalMuonAtDCA.getSigma2Y(); - // cXYatDCA = globalMuonAtDCA.getSigmaXY(); - // dcaX = globalMuonAtDCA.getX() - collision.posX(); - // dcaY = globalMuonAtDCA.getY() - collision.posY(); - // dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - // det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat - // if (det < 0) { - // dcaXYinSigma = 999.f; - // } else { - // dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2. * dcaX * dcaY * cXYatDCA) / det / 2.)); // dca xy in sigma - // } - // sigma_dcaXY = dcaXY / dcaXYinSigma; - - // cXXatDCA = mftsaAtDCA.getSigma2X(); - // cYYatDCA = mftsaAtDCA.getSigma2Y(); - // cXYatDCA = mftsaAtDCA.getSigmaXY(); - // dcaX = mftsaAtDCA.getX() - collision.posX(); - // dcaY = mftsaAtDCA.getY() - collision.posY(); - // dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - // det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat - // if (det < 0) { - // dcaXYinSigma = 999.f; - // } else { - // dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2. * dcaX * dcaY * cXYatDCA) / det / 2.)); // dca xy in sigma - // } - // sigma_dcaXY = dcaXY / dcaXYinSigma; + pt = propmuonAtPV.getPt(); + eta = propmuonAtPV.getEta(); + phi = propmuonAtPV.getPhi(); + o2::math_utils::bringTo02Pi(phi); - getDeltaEtaDeltaPhiAtMatchingPlane(collision, fwdtrack, mftCovs, detaMP, dphiMP); - o2::math_utils::bringToPMPi(dphiMP); - } - } else if (propagateToDCAhelix) { - auto fwdTrackParCov = getTrackParCovFwd(fwdtrack, fwdtrack); // values at innermost update - std::array dcaInfOrig{999.f, 999.f, 999.f}; - fwdTrackParCov.propagateToDCAhelix(mBz, {collision.posX(), collision.posY(), collision.posZ()}, dcaInfOrig); + if constexpr (withMFTCov) { + auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, 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); - eta = fwdTrackParCov.getEta(); - phi = fwdTrackParCov.getPhi(); + pt = globalMuonAtPV.getPt(); + eta = globalMuonAtPV.getEta(); + phi = globalMuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phi); - // pt = fwdTrackParCov.getPt(); - pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); - cXXatDCA = fwdTrackParCov.getSigma2X(); - cYYatDCA = fwdTrackParCov.getSigma2Y(); - cXYatDCA = fwdTrackParCov.getSigmaXY(); - - dcaX = fwdTrackParCov.getX() - collision.posX(); - dcaY = fwdTrackParCov.getY() - collision.posY(); - dcaZ = fwdTrackParCov.getZ() - collision.posZ(); + cXX = globalMuonAtPV.getSigma2X(); + cYY = globalMuonAtPV.getSigma2Y(); + cXY = globalMuonAtPV.getSigmaXY(); + dcaX = globalMuonAtPV.getX() - collision.posX(); + dcaY = globalMuonAtPV.getY() - collision.posY(); + dcaZ = globalMuonAtPV.getZ() - collision.posZ(); dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat - - dcaXYinSigma = 999.f; + det = cXX * cYY - cXY * cXY; // determinanat if (det < 0) { dcaXYinSigma = 999.f; } else { - dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2. * dcaX * dcaY * cXYatDCA) / det / 2.)); // dca xy in sigma + dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYY + dcaY * dcaY * cXX - 2. * dcaX * dcaY * cXY) / det / 2.)); // dca xy in sigma } sigma_dcaXY = dcaXY / dcaXYinSigma; + + getDeltaEtaDeltaPhiAtMatchingPlane(collision, fwdtrack, mftCovs, detaMP, dphiMP); + o2::math_utils::bringToPMPi(dphiMP); + } + + if (refitGlobalMuon) { + pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); } float ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); @@ -569,16 +505,16 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/correct/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxyz"), dcaXY, dcaZ); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxyinSigma"), dcaXYinSigma); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hMCHBitMap"), fwdtrack.mchBitMap()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hMIDBitMap"), fwdtrack.midBitMap()); - fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxResolutionvsPt"), pt, std::sqrt(cXXatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxResolutionvsPt"), pt, std::sqrt(cXX) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAyResolutionvsPt"), pt, std::sqrt(cYY) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/primary/correct/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt()); @@ -607,16 +543,16 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxyz"), dcaXY, dcaZ); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxyinSigma"), dcaXYinSigma); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hMCHBitMap"), fwdtrack.mchBitMap()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hMIDBitMap"), fwdtrack.midBitMap()); - fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxResolutionvsPt"), pt, std::sqrt(cXXatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxResolutionvsPt"), pt, std::sqrt(cXX) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAyResolutionvsPt"), pt, std::sqrt(cYY) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt()); @@ -647,16 +583,16 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxyz"), dcaXY, dcaZ); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxyinSigma"), dcaXYinSigma); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hMCHBitMap"), fwdtrack.mchBitMap()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hMIDBitMap"), fwdtrack.midBitMap()); - fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxResolutionvsPt"), pt, std::sqrt(cXXatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxResolutionvsPt"), pt, std::sqrt(cXX) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAyResolutionvsPt"), pt, std::sqrt(cYY) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt()); @@ -685,16 +621,16 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxyz"), dcaXY, dcaZ); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxyinSigma"), dcaXYinSigma); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hMCHBitMap"), fwdtrack.mchBitMap()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hMIDBitMap"), fwdtrack.midBitMap()); - fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxResolutionvsPt"), pt, std::sqrt(cXXatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um - fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxResolutionvsPt"), pt, std::sqrt(cXX) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAyResolutionvsPt"), pt, std::sqrt(cYY) * 1e+4); // convert cm to um + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt());