From 9bde930e7e5c397ceb0fe624fca58fcca9cc58ec Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Tue, 11 Aug 2026 16:57:47 +0200 Subject: [PATCH 1/2] [PWGDQ] add dimuon analysis in global alignment task Added comparison of dimuon properties (invariant mass, angles, DCA) for different alignment configurations. --- PWGDQ/Tasks/muonGlobalAlignment.cxx | 1695 +++++++++++++++++---------- 1 file changed, 1081 insertions(+), 614 deletions(-) diff --git a/PWGDQ/Tasks/muonGlobalAlignment.cxx b/PWGDQ/Tasks/muonGlobalAlignment.cxx index 7014a0ea0c7..dd8cc3bee63 100644 --- a/PWGDQ/Tasks/muonGlobalAlignment.cxx +++ b/PWGDQ/Tasks/muonGlobalAlignment.cxx @@ -9,20 +9,23 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -/// \file muonDCA.cxx -/// \brief Task to compute and evaluate DCA quantities -/// \author Nicolas Bizé , SUBATECH -// +/// \file muonGlobalAlignment.cxx +/// \brief Analysis of global alignment between MFT and MCH-MID +/// \author Andrea Ferrero , CEA-Saclay #include "PWGDQ/Core/VarManager.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/Core/RecoDecay.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" +#include "Common/Core/fwdtrackUtilities.h" #include #include +#include +#include #include #include #include @@ -80,8 +83,6 @@ #include #include -#include - using namespace o2; using namespace o2::mch; using namespace o2::framework; @@ -110,7 +111,6 @@ using MyMFTCovariance = MyMFTCovariances::iterator; using SMatrix55 = ROOT::Math::SMatrix>; using SMatrix5 = ROOT::Math::SVector; -static o2::globaltracking::MatchGlobalFwd sExtrap; using o2::dataformats::GlobalFwdTrack; using o2::track::TrackParCovFwd; @@ -126,86 +126,100 @@ DECLARE_SOA_TABLE(CompactMFTTracks, "AOD", "COMPACTMFT", //! standalone table fo using CompactMFTTrack = CompactMFTTracks; } // namespace o2::aod -struct muonGlobalAlignment { +struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struct (exception) + + static constexpr int GlobalTrackTypeMax = 2; + static constexpr int NMchChambers = 10; + static constexpr int NMchDetElems = 156; + static constexpr int ThetaAbsBoundaryDeg = 3; + static constexpr double SlopeResolutionZ = 535.; + static constexpr double AbsorberBackZ = -505.f; + static constexpr double BransonPlaneZ = -466.f; Produces mftTable; Configurable cfgProduceMFTTable{"cfgProduceMFTTable", false, "flag to produce MFTsa table"}; //// Variables for selecting MCH and MFT tracks - Configurable fTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; - Configurable fPtMchLow{"cfgPtMchLow", 0.7f, ""}; - Configurable fEtaMftLow{"cfgEtaMftlow", -3.6f, ""}; - Configurable fEtaMftUp{"cfgEtaMftup", -2.5f, ""}; - Configurable fRabsLow{"cfgRabsLow", 17.6f, ""}; - Configurable fRabsUp{"cfgRabsUp", 89.5f, ""}; - Configurable fSigmaPdcaUp{"cfgPdcaUp", 6.f, ""}; + Configurable cfgTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; + Configurable cfgPtMchLow{"cfgPtMchLow", 0.7f, ""}; + Configurable cfgEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; + Configurable cfgEtaMchUp{"cfgEtaMchUp", -2.5f, ""}; + Configurable cfgEtaMftLow{"cfgEtaMftLow", -3.6f, ""}; + Configurable cfgEtaMftUp{"cfgEtaMftUp", -2.5f, ""}; + Configurable cfgRabsLow{"cfgRabsLow", 17.6f, ""}; + Configurable cfgRabsUp{"cfgRabsUp", 89.5f, ""}; + Configurable fSigmaPdcaUp{"fSigmaPdcaUp", 6.f, ""}; + + Configurable cfgTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; + Configurable cfgTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; - Configurable fTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; - Configurable fTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; + Configurable cfgMftDcaMatchChi2Up{"cfgMftDcaMatchChi2Up", 10.f, ""}; - Configurable fMftMchResidualsPLow{"cfgMftMchResidualsPLow", 30.f, ""}; - Configurable fMftMchResidualsPtLow{"cfgMftMchResidualsPtLow", 4.f, ""}; + Configurable cfgMftMchResidualsPLow{"cfgMftMchResidualsPLow", 30.f, ""}; + Configurable cfgMftMchResidualsPtLow{"cfgMftMchResidualsPtLow", 4.f, ""}; - Configurable fMftTracksMultiplicityMax{"cfgMftTracksMultiplicityMax", 0, "Maximum number of MFT tracks to be processed per event (zero means no limit)"}; + Configurable cfgMftTracksMultiplicityMax{"cfgMftTracksMultiplicityMax", 0, "Maximum number of MFT tracks to be processed per event (zero means no limit)"}; // Magnetic field position bias Configurable cfgFieldOriginBiasZ{"cfgFieldOriginBiasZ", 0.0f, "Bias applied to the magnetic field z position"}; - Configurable fVertexZshift{"cfgVertexZshift", 0.0f, "Correction to the vertex z position"}; - Configurable fDipoleZshift{"cfgDipoleZshift", 0.0f, "Correction to the dipole z position"}; + Configurable cfgDipoleZshift{"cfgDipoleZshift", 0.0f, "Correction to the dipole z position"}; + Configurable cfgVertexZshift{"cfgVertexZshift", 0.0f, "Correction to the vertex z position"}; //// Variables for MFT alignment corrections struct : ConfigurableGroup { - Configurable fEnableMFTAlignmentCorrections{"cfgEnableMFTAlignmentCorrections", false, ""}; + Configurable cfgEnableMFTAlignmentCorrections{"cfgEnableMFTAlignmentCorrections", false, ""}; // slope corrections - Configurable fMFTAlignmentCorrXSlopeTop{"cfgMFTAlignmentCorrXSlopeTop", (-0.0006696 - 0.0005621) / 2.f, "MFT X slope correction - top half"}; - Configurable fMFTAlignmentCorrXSlopeBottom{"cfgMFTAlignmentCorrXSlopeBottom", (0.00105 + 0.001007) / 2.f, "MFT X slope correction - bottom half"}; - Configurable fMFTAlignmentCorrYSlopeTop{"cfgMFTAlignmentCorrYSlopeTop", (-0.002299 - 0.002442) / 2.f, "MFT Y slope correction - top half"}; - Configurable fMFTAlignmentCorrYSlopeBottom{"cfgMFTAlignmentCorrYSlopeBottom", (-0.0005339 - 0.0006921) / 2.f, "MFT Y slope correction - bottom half"}; + Configurable cfgMFTAlignmentCorrXSlopeTop{"cfgMFTAlignmentCorrXSlopeTop", (-0.0006696 - 0.0005621) / 2.f, "MFT X slope correction - top half"}; + Configurable cfgMFTAlignmentCorrXSlopeBottom{"cfgMFTAlignmentCorrXSlopeBottom", (0.00105 + 0.001007) / 2.f, "MFT X slope correction - bottom half"}; + Configurable cfgMFTAlignmentCorrYSlopeTop{"cfgMFTAlignmentCorrYSlopeTop", (-0.002299 - 0.002442) / 2.f, "MFT Y slope correction - top half"}; + Configurable cfgMFTAlignmentCorrYSlopeBottom{"cfgMFTAlignmentCorrYSlopeBottom", (-0.0005339 - 0.0006921) / 2.f, "MFT Y slope correction - bottom half"}; // offset corrections - Configurable fMFTAlignmentCorrXOffsetTop{"cfgMFTAlignmentCorrXOffsetTop", 0.f, "MFT X offset correction - top half"}; - Configurable fMFTAlignmentCorrXOffsetBottom{"cfgMFTAlignmentCorrXOffsetBottom", 0.f, "MFT X offset correction - bottom half"}; - Configurable fMFTAlignmentCorrYOffsetTop{"cfgMFTAlignmentCorrYOffsetTop", 0.f, "MFT Y offset correction - top half"}; - Configurable fMFTAlignmentCorrYOffsetBottom{"cfgMFTAlignmentCorrYOffsetBottom", 0.f, "MFT Y offset correction - bottom half"}; + Configurable cfgMFTAlignmentCorrXOffsetTop{"cfgMFTAlignmentCorrXOffsetTop", 0.f, "MFT X offset correction - top half"}; + Configurable cfgMFTAlignmentCorrXOffsetBottom{"cfgMFTAlignmentCorrXOffsetBottom", 0.f, "MFT X offset correction - bottom half"}; + Configurable cfgMFTAlignmentCorrYOffsetTop{"cfgMFTAlignmentCorrYOffsetTop", 0.f, "MFT Y offset correction - top half"}; + Configurable cfgMFTAlignmentCorrYOffsetBottom{"cfgMFTAlignmentCorrYOffsetBottom", 0.f, "MFT Y offset correction - bottom half"}; } configMFTAlignmentCorrections; //// Variables for re-alignment setup struct : ConfigurableGroup { - Configurable fEnableMCHRealign{"cfgEnableMCHRealign", true, "Enable re-alignment of MCH clusters and tracks"}; - Configurable fChamberResolutionX{"cfgChamberResolutionX", 0.4, "Chamber resolution along X configuration for refit"}; // 0.4cm pp, 0.2cm PbPb - Configurable fChamberResolutionY{"cfgChamberResolutionY", 0.4, "Chamber resolution along Y configuration for refit"}; // 0.4cm pp, 0.2cm PbPb - Configurable fSigmaCutImprove{"cfgSigmaCutImprove", 6., "Sigma cut for track improvement"}; - Configurable fMCHRealignCorrections{"cfgMCHRealignCorrections", "", "MCH DE positions/angles corrections in JSON format"}; + Configurable cfgEnableMCHRefit{"cfgEnableMCHRefit", false, "Enable re-fitting of MCH tracks"}; + Configurable cfgEnableMCHRealign{"cfgEnableMCHRealign", false, "Enable re-alignment of MCH clusters and tracks"}; + Configurable cfgChamberResolutionX{"cfgChamberResolutionX", 0.4, "Chamber resolution along X configuration for refit"}; // 0.4cm pp, 0.2cm PbPb + Configurable cfgChamberResolutionY{"cfgChamberResolutionY", 0.4, "Chamber resolution along Y configuration for refit"}; // 0.4cm pp, 0.2cm PbPb + Configurable cfgSigmaCutImprove{"cfgSigmaCutImprove", 6., "Sigma cut for track improvement"}; + Configurable cfgMCHRealignCorrections{"cfgMCHRealignCorrections", "", "MCH DE positions/angles corrections in JSON format"}; } configRealign; //// Variables for ccdb struct : ConfigurableGroup { - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; // Configurable geoPathRealign{"geoPathRealign", "Users/j/jcastill/GeometryAlignedFix10Fix15ShiftCh1BNew2", "Path of the geometry file"}; Configurable geoPathRealign{"geoPathRealign", "Users/j/jcastill/GeometryAlignedLoczzm4pLHC24anap1sR5a", "Path of the geometry file"}; - Configurable nolaterthan{"ccdb-no-later-than-ref", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of reference basis"}; - Configurable nolaterthanRealign{"ccdb-no-later-than-new", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of new basis"}; + Configurable cfgCcdbNoLaterThanRef{"cfgCcdbNoLaterThanRef", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of reference basis"}; + Configurable cfgCcdbNoLaterThanNew{"cfgCcdbNoLaterThanNew", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of new basis"}; } configCCDB; - Configurable fRequireGoodRCT{"cfgRequireGoodRCT", true, "Require good detector flags in Run Condition Table"}; + Configurable cfgRequireGoodRCT{"cfgRequireGoodRCT", true, "Require good detector flags in Run Condition Table"}; - Configurable fEnableVertexShiftAnalysis{"cfgEnableVertexShiftAnalysis", true, "Enable the analysis of vertex shift"}; - Configurable fEnableMftDcaAnalysis{"cfgEnableMftDcaAnalysis", true, "Enable the analysis of DCA-based MFT alignment"}; - Configurable fEnableMftDcaExtraPlots{"cfgEnableMftDcaExtraPlots", false, "Enable additional plots for the analysis of DCA-based MFT alignment"}; - Configurable fEnableGlobalFwdDcaAnalysis{"cfgEnableGlobalFwdDcaAnalysis", true, "Enable the analysis of DCA-based MFT alignment using global forward tracks"}; - Configurable fEnableMftMchResidualsAnalysis{"cfgEnableMftMchResidualsAnalysis", true, "Enable the analysis of residuals between MFT tracks and MCH clusters"}; - Configurable fEnableMftMchResidualsExtraPlots{"cfgEnableMftMchResidualsExtraPlots", false, "Enable additional plots for the analysis of residuals between MFT tracks and MCH clusters"}; - Configurable fEnableMftMchMatchingAnalysis{"cfgEnableMftMchMatchingAnalysis", false, "Enable the analysis of residuals between MFT and MCH tracks at reference planes"}; + Configurable cfgEnableVertexShiftAnalysis{"cfgEnableVertexShiftAnalysis", false, "Enable the analysis of vertex shift"}; + Configurable cfgEnableMftDcaAnalysis{"cfgEnableMftDcaAnalysis", false, "Enable the analysis of DCA-based MFT alignment"}; + Configurable cfgEnableMftDcaExtraPlots{"cfgEnableMftDcaExtraPlots", false, "Enable additional plots for the analysis of DCA-based MFT alignment"}; + Configurable cfgEnableGlobalFwdDcaAnalysis{"cfgEnableGlobalFwdDcaAnalysis", false, "Enable the analysis of DCA-based MFT alignment using global forward tracks"}; + Configurable cfgEnableMftMchResidualsAnalysis{"cfgEnableMftMchResidualsAnalysis", true, "Enable the analysis of residuals between MFT tracks and MCH clusters"}; + Configurable cfgEnableMftMchResidualsExtraPlots{"cfgEnableMftMchResidualsExtraPlots", false, "Enable additional plots for the analysis of residuals between MFT tracks and MCH clusters"}; + Configurable cfgEnableMftMchMatchingAnalysis{"cfgEnableMftMchMatchingAnalysis", false, "Enable the analysis of residuals between MFT and MCH tracks at reference planes"}; + Configurable cfgEnableDimuonAnalysis{"cfgEnableDimuonAnalysis", false, "Enable the analysis of di-muon pairs"}; - Configurable fRefPlaneZMFT{"cfgRefPlaneZMFT", o2::mft::constants::mft::LayerZCoordinate()[0], "Reference plane on MFT side"}; - Configurable fRefPlaneZMCH{"cfgRefPlaneZMCH", -526.0, "Reference plane on MCH side"}; + Configurable cfgRefPlaneZMFT{"cfgRefPlaneZMFT", o2::mft::constants::mft::LayerZCoordinate()[0], "Reference plane on MFT side"}; + Configurable cfgRefPlaneZMCH{"cfgRefPlaneZMCH", -526.0, "Reference plane on MCH side"}; int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field - Service ccdbManager; - o2::field::MagneticField* fieldB; + Service ccdbManager{}; + o2::field::MagneticField* fieldB{nullptr}; o2::ccdb::CcdbApi ccdbApi; // Derived version of mch::Track class that handles the associated clusters as internal objects and deletes them in the destructor @@ -224,13 +238,59 @@ struct muonGlobalAlignment { } }; + class TrackParExt : public o2::track::TrackParCovFwd + { + public: + TrackParExt() = default; + TrackParExt(const TrackParExt& t) = default; + explicit TrackParExt(o2::track::TrackParCovFwd const& t, int nc = -1, bool r = false) + : TrackParCovFwd(t), nClusters(nc), removable(r) {} + ~TrackParExt() = default; + + TrackParExt& operator=(const TrackParCovFwd& tpf) + { + o2::track::TrackParCovFwd::operator=(tpf); + return *this; + } + TrackParExt& operator=(const TrackParExt& tpe) + { + o2::track::TrackParCovFwd::operator=(tpe); + nClusters = tpe.getNClusters(); + removable = tpe.isRemovable(); + return *this; + } + + void setNClusters(int n) { nClusters = n; } + [[nodiscard]] int getNClusters() const { return nClusters; } + + void setRemovable() { removable = true; } + [[nodiscard]] bool isRemovable() const { return removable; } + + [[nodiscard]] o2::track::TrackParCovFwd asTrackParCovFwd() const + { + return {static_cast(*this)}; + } + + private: + int nClusters{-1}; + bool removable{false}; + }; + + std::unordered_map mMchTrackPars; + std::unordered_map mMftTrackPars; + std::unordered_map mMchTrackParsNew; + std::unordered_map mMftTrackParsNew; + + using MuonPair = std::pair; + using GlobalMuonPair = std::pair, std::vector>; + geo::TransformationCreator transformation; std::map transformRef; // reference geometry w.r.t track data std::map transformNew; // new geometry TGeoManager* geoNew = nullptr; TGeoManager* geoRef = nullptr; - TrackFitter trackFitter; // Track fitter from MCH tracking library - double mImproveCutChi2; // Chi2 cut for track improvement. + TrackFitter trackFitter; // Track fitter from MCH tracking library + double mImproveCutChi2{0}; // Chi2 cut for track improvement. struct AlignmentCorrections { double x{0}; @@ -268,100 +328,12 @@ struct muonGlobalAlignment { std::map> globalMuonTracks; }; - void InitCollisions(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - std::map& collisionInfos) - { - // fill collision information for global muon tracks (MFT-MCH-MID matches) - for (auto muonTrack : muonTracks) { - if (!muonTrack.has_collision()) - continue; - - auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); - - if (fRequireGoodRCT && !rctChecker(collision)) - continue; - - uint64_t collisionIndex = collision.globalIndex(); - - auto bc = bcs.rawIteratorAt(collision.bcId()); - - auto& collisionInfo = collisionInfos[collisionIndex]; - collisionInfo.bc = bc.globalBC(); - collisionInfo.zVertex = collision.posZ(); - - if (static_cast(muonTrack.trackType()) > 2) { - // standalone MCH or MCH-MID tracks - uint64_t mchTrackIndex = muonTrack.globalIndex(); - collisionInfo.mchTracks.push_back(mchTrackIndex); - } else { - // global muon tracks (MFT-MCH or MFT-MCH-MID) - uint64_t muonTrackIndex = muonTrack.globalIndex(); - auto const& mchTrack = muonTrack.template matchMCHTrack_as(); - uint64_t mchTrackIndex = mchTrack.globalIndex(); - - // check if a vector of global muon candidates is already available for the current MCH index - // if not, initialize a new one and add the current global muon track - // bool globalMuonTrackFound = false; - auto matchingCandidateIterator = collisionInfo.globalMuonTracks.find(mchTrackIndex); - if (matchingCandidateIterator != collisionInfo.globalMuonTracks.end()) { - matchingCandidateIterator->second.push_back(muonTrackIndex); - // globalMuonTrackFound = true; - } else { - collisionInfo.globalMuonTracks[mchTrackIndex].push_back(muonTrackIndex); - } - } - } - - // sort the vectors of matching candidates in ascending order based on the matching chi2 value - auto compareChi2 = [&muonTracks](uint64_t trackIndex1, uint64_t trackIndex2) -> bool { - auto const& track1 = muonTracks.rawIteratorAt(trackIndex1); - auto const& track2 = muonTracks.rawIteratorAt(trackIndex2); - - return (track1.chi2MatchMCHMFT() < track2.chi2MatchMCHMFT()); - }; - - for (auto& [collisionIndex, collisionInfo] : collisionInfos) { - for (auto& [mchIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { - std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareChi2); - } - } - } - - void InitCollisions(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - MyMFTs const& mftTracks, - std::map& collisionInfos) - { - InitCollisions(collisions, bcs, muonTracks, collisionInfos); - - // fill collision information for MFT standalone tracks - for (auto mftTrack : mftTracks) { - if (!mftTrack.has_collision()) - continue; - - auto collision = collisions.rawIteratorAt(mftTrack.collisionId()); - uint64_t collisionIndex = collision.globalIndex(); - - auto bc = bcs.rawIteratorAt(collision.bcId()); - - uint64_t mftTrackIndex = mftTrack.globalIndex(); - - auto& collisionInfo = collisionInfos[collisionIndex]; - collisionInfo.bc = bc.globalBC(); - collisionInfo.zVertex = collision.posZ(); - - collisionInfo.mftTracks.push_back(mftTrackIndex); - } - } - template void initCCDB(BC const& bc) { - if (mRunNumber == bc.runNumber()) + if (mRunNumber == bc.runNumber()) { return; + } mRunNumber = bc.runNumber(); ccdbManager->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); @@ -372,31 +344,31 @@ struct muonGlobalAlignment { TrackExtrap::setField(); TrackExtrap::useExtrapV2(); fieldB = static_cast(TGeoGlobalMagField::Instance()->GetField()); // for MFT - double centerMFT[3] = {0, 0, -61.4}; // or use middle point between Vtx and MFT? - mBzAtMftCenter = fieldB->getBz(centerMFT); + std::array centerMFT{0, 0, -61.4}; // or use middle point between Vtx and MFT? + mBzAtMftCenter = fieldB->getBz(centerMFT.data()); } else { LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", bc.timestamp()); } // Load geometry information from CCDB/local - LOGF(info, "Loading reference aligned geometry from CCDB no later than %d", configCCDB.nolaterthan.value); - ccdbManager->setCreatedNotAfter(configCCDB.nolaterthan); // this timestamp has to be consistent with what has been used in reco + LOGF(info, "Loading reference aligned geometry from CCDB no later than %d", configCCDB.cfgCcdbNoLaterThanRef.value); + ccdbManager->setCreatedNotAfter(configCCDB.cfgCcdbNoLaterThanRef); // this timestamp has to be consistent with what has been used in reco geoRef = ccdbManager->getForTimeStamp(configCCDB.geoPath, bc.timestamp()); ccdbManager->clearCache(configCCDB.geoPath); - if (configRealign.fEnableMCHRealign && fEnableMftMchResidualsAnalysis) { + if (configRealign.cfgEnableMCHRealign && (cfgEnableMftMchResidualsAnalysis)) { if (geoRef != nullptr) { transformation = geo::transformationFromTGeoManager(*geoRef); } else { LOGF(fatal, "Reference aligned geometry object is not available in CCDB at timestamp=%llu", bc.timestamp()); } - for (int i = 0; i < 156; i++) { + for (int i = 0; i < NMchDetElems; i++) { int iDEN = GetDetElemId(i); transformRef[iDEN] = transformation(iDEN); } - LOGF(info, "Loading new aligned geometry from CCDB no later than %d", configCCDB.nolaterthanRealign.value); - ccdbManager->setCreatedNotAfter(configCCDB.nolaterthanRealign); // make sure this timestamp can be resolved regarding the reference one + LOGF(info, "Loading new aligned geometry from CCDB no later than %d", configCCDB.cfgCcdbNoLaterThanNew.value); + ccdbManager->setCreatedNotAfter(configCCDB.cfgCcdbNoLaterThanNew); // make sure this timestamp can be resolved regarding the reference one geoNew = ccdbManager->getForTimeStamp(configCCDB.geoPathRealign, bc.timestamp()); ccdbManager->clearCache(configCCDB.geoPathRealign); if (geoNew != nullptr) { @@ -404,7 +376,7 @@ struct muonGlobalAlignment { } else { LOGF(fatal, "New aligned geometry object is not available in CCDB at timestamp=%llu", bc.timestamp()); } - for (int i = 0; i < 156; i++) { + for (int i = 0; i < NMchDetElems; i++) { int iDEN = GetDetElemId(i); transformNew[iDEN] = transformation(iDEN); } @@ -414,10 +386,10 @@ struct muonGlobalAlignment { void init(o2::framework::InitContext&) { // Load geometry - ccdbManager->setURL(configCCDB.ccdburl); + ccdbManager->setURL(configCCDB.ccdbUrl); ccdbManager->setCaching(true); ccdbManager->setLocalObjectValidityChecking(); - ccdbApi.init(configCCDB.ccdburl); + ccdbApi.init(configCCDB.ccdbUrl); mRunNumber = 0; // configure magnetic field position bias @@ -426,10 +398,10 @@ struct muonGlobalAlignment { // Configuration for track fitter const auto& trackerParam = TrackerParam::Instance(); trackFitter.setBendingVertexDispersion(trackerParam.bendingVertexDispersion); - trackFitter.setChamberResolution(configRealign.fChamberResolutionX, configRealign.fChamberResolutionY); + trackFitter.setChamberResolution(configRealign.cfgChamberResolutionX, configRealign.cfgChamberResolutionY); trackFitter.smoothTracks(true); trackFitter.useChamberResolution(); - mImproveCutChi2 = 2. * configRealign.fSigmaCutImprove * configRealign.fSigmaCutImprove; + mImproveCutChi2 = 2. * configRealign.cfgSigmaCutImprove * configRealign.cfgSigmaCutImprove; // use the Runge-Kutta extrapolation v2 TrackExtrap::useExtrapV2(); @@ -437,7 +409,7 @@ struct muonGlobalAlignment { // Fill table of MCH alignment corrections rapidjson::Document document; // Check that the json is parsed correctly - rapidjson::ParseResult jsonOk = document.Parse(configRealign.fMCHRealignCorrections.value.c_str()); + rapidjson::ParseResult jsonOk = document.Parse(configRealign.cfgMCHRealignCorrections.value.c_str()); if (jsonOk) { for (rapidjson::Value::ConstMemberIterator it = document.MemberBegin(); it != document.MemberEnd(); it++) { LOG(info) << "DE" << it->name.GetString() << " alignment corrections:"; @@ -445,10 +417,10 @@ struct muonGlobalAlignment { LOG(info) << " y: " << it->value["y"].GetDouble(); LOG(info) << " z: " << it->value["z"].GetDouble(); - mMchAlignmentCorrections[std::stoi(it->name.GetString())] = { - it->value["x"].GetDouble(), - it->value["y"].GetDouble(), - it->value["z"].GetDouble()}; + mMchAlignmentCorrections[std::stoi(it->name.GetString())] = AlignmentCorrections{ + .x = it->value["x"].GetDouble(), + .y = it->value["y"].GetDouble(), + .z = it->value["z"].GetDouble()}; } } else { LOG(error) << "JSON parse error: " << rapidjson::GetParseErrorFunc(jsonOk.Code()) << " (" << jsonOk.Offset() << ")"; @@ -482,12 +454,12 @@ struct muonGlobalAlignment { registry.add("vertex_y_vs_x", std::format("Vertex y vs. x").c_str(), {HistType::kTH2F, {vxAxis, vyAxis}}); registry.add("vertex_z", std::format("Vertex z").c_str(), {HistType::kTH1F, {vzAxis}}); - if (fEnableVertexShiftAnalysis || fEnableMftDcaAnalysis) { + if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { registry.add("DCA/MFT/nTracksMFT", std::format("Number of MFT tracks per collision").c_str(), {HistType::kTH1F, {{100, 0, 1000, "# of MFT tracks"}}}); registry.add("DCA/MFT/DCA_y_vs_x", std::format("DCA y vs. x").c_str(), {HistType::kTH2F, {dcaxMFTAxis, dcayMFTAxis}}); } - if (fEnableVertexShiftAnalysis) { + if (cfgEnableVertexShiftAnalysis) { registry.add("DCA/MFT/DCA_x_vs_phi_vs_zshift", std::format("DCA(x) vs. #phi vs. z shift").c_str(), {HistType::kTH3F, {zshiftAxis, phiAxis, dcaxMFTAxis}}); registry.add("DCA/MFT/DCA_y_vs_phi_vs_zshift", std::format("DCA(y) vs. #phi vs. z shift").c_str(), {HistType::kTH3F, {zshiftAxis, phiAxis, dcayMFTAxis}}); @@ -497,13 +469,13 @@ struct muonGlobalAlignment { registry.add("DCA/MFT/DCA_y_vs_slopey_vs_zshift", std::format("DCA(y) vs. y slope vs. z shift").c_str(), {HistType::kTH3F, {zshiftAxis, syAxis, dcayMFTAxis}}); } - if (fEnableMftDcaAnalysis) { + if (cfgEnableMftDcaAnalysis) { registry.add("DCA/MFT/DCA_x", "DCA(x) vs. vz, tx, ty, nclus", HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); registry.add("DCA/MFT/DCA_y", "DCA(y) vs. vz, tx, ty, nclus", HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); - if (fEnableMftDcaExtraPlots) { + if (cfgEnableMftDcaExtraPlots) { registry.add("DCA/MFT/layers", "Layers vs. tx, ty, nclus", HistType::kTHnSparseF, {mftLayerAxis, txAxis, tyAxis, nMftClustersAxis}); registry.add("DCA/MFT/trackChi2", "Track #chi^{2} vs. tx, ty, nclus, layer", @@ -523,16 +495,16 @@ struct muonGlobalAlignment { } } - if (fEnableGlobalFwdDcaAnalysis) { + if (cfgEnableGlobalFwdDcaAnalysis) { registry.add("DCA/GlobalFwd/DCA_x", "DCA(x) vs. vz, tx, ty, nclus", HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); registry.add("DCA/GlobalFwd/DCA_y", "DCA(y) vs. vz, tx, ty, nclus", HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); } - if (fEnableMftMchResidualsAnalysis) { - AxisSpec dxAxis = {400, -20.0, 20.0, "#Delta x (cm)"}; - AxisSpec dyAxis = {400, -20.0, 20.0, "#Delta y (cm)"}; + if (cfgEnableMftMchResidualsAnalysis) { + AxisSpec dxAxis = {400, -10.0, 10.0, "#Delta x (cm)"}; + AxisSpec dyAxis = {400, -10.0, 10.0, "#Delta y (cm)"}; registry.add("DCA/MCH/DCA_y_vs_x", std::format("DCA y vs. x").c_str(), {HistType::kTH2F, {dcaxMCHAxis, dcayMCHAxis}}); registry.add("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom", std::format("DCA(x) vs. p, quadrant, chargeSign").c_str(), {HistType::kTHnSparseF, {{20, 0, 100.0, "p (GeV/c)"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dcaxMCHAxis}}); @@ -574,7 +546,7 @@ struct muonGlobalAlignment { registry.get(HIST("residuals/de_alignment_corrections_y"))->SetBinError(deIndex + 1, 0.1); } - if (fEnableMftMchResidualsExtraPlots) { + if (cfgEnableMftMchResidualsExtraPlots) { registry.add("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz", std::format("DCA(x) vs. vz, quadrant, chargeSign").c_str(), {HistType::kTHnSparseF, {dcazAxis, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dcaxMCHAxis}}); registry.add("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_vz", std::format("DCA(y) vs. vz, quadrant, chargeSign").c_str(), {HistType::kTHnSparseF, {dcazAxis, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dcayMCHAxis}}); @@ -583,7 +555,7 @@ struct muonGlobalAlignment { } } - if (fEnableMftMchMatchingAnalysis) { + if (cfgEnableMftMchMatchingAnalysis) { AxisSpec dxAxis = {200, -10.0, 10.0, "#Deltax (cm)"}; AxisSpec dyAxis = {200, -10.0, 10.0, "#Deltay (cm)"}; AxisSpec dsxAxis = {200, -0.1, 0.1, "#Deltaslope(x) (rad)"}; @@ -614,24 +586,66 @@ struct muonGlobalAlignment { registry.add("matching/dphiAtMCH", "Tracks #Delta#phi at MCH reference plane", {HistType::kTHnSparseF, {dphiAxis, {20, -100.0, 100.0, "track x (cm)"}, {20, -100.0, 100.0, "track y (cm)"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}}}); } + + if (cfgEnableDimuonAnalysis) { + AxisSpec invMassAxis = {1500, 0, 15, "M_{#mu^{+}#mu^{-}} (GeV/c^{2})"}; + AxisSpec pTAxis = {30, 0, 30, "#mu^{+}#mu^{-} p_{T} (GeV/c)"}; + AxisSpec pAxis = {50, 0, 200, "#mu^{+}#mu^{-} p (GeV/c)"}; + AxisSpec muPosQuadrantAxis = {4, 0, 4, "#mu^{+} quadrant"}; + AxisSpec muNegQuadrantAxis = {4, 0, 4, "#mu^{-} quadrant"}; + AxisSpec angleAxis = {100, 0, 0.5, "#mu^{+}#mu^{-} angle (rad)"}; + AxisSpec angleDiffAxis = {500, -0.05, 0.05, "#mu^{+}#mu^{-} angle difference (rad)"}; + AxisSpec dcaxAxis = {400, -10.0, 10.0, "#mu^{+}#mu^{-} DCA_{x} (cm)"}; + AxisSpec dcayAxis = {400, -10.0, 10.0, "#mu^{+}#mu^{-} DCA_{y} (cm)"}; + AxisSpec dcaMftxAxis = {400, -0.5, 0.5, "#mu^{+}#mu^{-} DCA_{x} (cm)"}; + AxisSpec dcaMftyAxis = {400, -0.5, 0.5, "#mu^{+}#mu^{-} DCA_{y} (cm)"}; + + // di-muon invariant mass distributions + registry.add("dimuon/invariantMass_MuonKine_MuonCuts", "#mu^{+}#mu^{-} invariant mass (muon cuts)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // difference in mu+mu- opening angle between MCH and global muon tracks + registry.add("dimuon/angle_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} opening angle difference (global muon cuts, good matches)", {HistType::kTHnSparseF, {angleDiffAxis, angleAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // mu+mu- DCA + registry.add("dimuon/dcax_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{x} (muon cuts)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcay_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{y} (muon cuts)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcay_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftyAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + + // di-muon invariant mass distributions (realigned/refitted tracks) + registry.add("dimuon/realign/invariantMass_MuonKine_MuonCuts", "#mu^{+}#mu^{-} invariant mass (muon cuts)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // difference in mu+mu- opening angle between MCH and global muon tracks (realigned/refitted tracks) + registry.add("dimuon/realign/angle_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} opening angle difference (global muon cuts, good matches)", {HistType::kTHnSparseF, {angleDiffAxis, angleAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // mu+mu- DCA + registry.add("dimuon/realign/dcax_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{x} (muon cuts)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcay_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{y} (muon cuts)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcay_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftyAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + } } int GetDetElemId(int iDetElemNumber) { - const int fgNCh = 10; - const int fgNDetElemCh[fgNCh] = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; - const int fgSNDetElemCh[fgNCh + 1] = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; + constexpr int fgNCh = 10; + const std::array fgNDetElemCh{4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; + const std::array fgSNDetElemCh{0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; // make sure detector number is valid - if (!(iDetElemNumber >= fgSNDetElemCh[0] && - iDetElemNumber < fgSNDetElemCh[10])) { + if (iDetElemNumber < fgSNDetElemCh[0] || + iDetElemNumber >= fgSNDetElemCh[10]) { LOGF(fatal, "Invalid detector element number: %d", iDetElemNumber); } /// get det element number from ID // get chamber and element number in chamber int iCh = 0; int iDet = 0; - for (int i = 1; i <= 10; i++) { + for (int i = 1; i <= NMchChambers; i++) { if (iDetElemNumber < fgSNDetElemCh[i]) { iCh = i; iDet = iDetElemNumber - fgSNDetElemCh[i - 1]; @@ -640,7 +654,7 @@ struct muonGlobalAlignment { } // make sure detector index is valid - if (!(iCh > 0 && iCh <= 10 && iDet < fgNDetElemCh[iCh - 1])) { + if (iCh <= 0 || iCh > NMchChambers || iDet >= fgNDetElemCh[iCh - 1]) { LOGF(fatal, "Invalid detector element id: %d", 100 * iCh + iDet); } @@ -683,7 +697,7 @@ struct muonGlobalAlignment { { static int nDE = 0; if (nDE <= 0) { - for (int c = 0; c < 10; c++) { + for (int c = 0; c < NMchChambers; c++) { nDE += getNumDEinChamber(c); } } @@ -713,18 +727,18 @@ struct muonGlobalAlignment { return idx + offset; } - int GetQuadrant(double phi) + int GetQuadrant(float phi) { - if (phi >= 0 && phi < 90) { + if (phi >= 0 && phi < o2::constants::math::PIHalf) { return 0; } - if (phi >= 90 && phi <= 180) { + if (phi >= o2::constants::math::PIHalf && phi <= o2::constants::math::PI) { return 1; } - if (phi >= -180 && phi < -90) { + if (phi >= -o2::constants::math::PI && phi < -o2::constants::math::PIHalf) { return 2; } - if (phi >= -90 && phi < 0) { + if (phi >= -o2::constants::math::PIHalf && phi < 0) { return 3; } return -1; @@ -733,8 +747,7 @@ struct muonGlobalAlignment { template int GetQuadrant(const T& track) { - double phi = track.phi() * 180 / TMath::Pi(); - return GetQuadrant(phi); + return GetQuadrant(static_cast(track.phi())); } template @@ -747,18 +760,16 @@ struct muonGlobalAlignment { // MCH track format. // Parameter conversion - double alpha1, alpha3, alpha4, x2, x3, x4; - - x2 = fwdtrack.getPhi(); - x3 = fwdtrack.getTanl(); - x4 = fwdtrack.getInvQPt(); + double x2 = fwdtrack.getPhi(); + double x3 = fwdtrack.getTanl(); + double x4 = fwdtrack.getInvQPt(); auto sinx2 = TMath::Sin(x2); auto cosx2 = TMath::Cos(x2); - alpha1 = cosx2 / x3; - alpha3 = sinx2 / x3; - alpha4 = x4 / TMath::Sqrt(x3 * x3 + sinx2 * sinx2); + double alpha1 = cosx2 / x3; + double alpha3 = sinx2 / x3; + double alpha4 = x4 / TMath::Sqrt(x3 * x3 + sinx2 * sinx2); auto K = TMath::Sqrt(x3 * x3 + sinx2 * sinx2); auto K3 = K * K * K; @@ -803,11 +814,11 @@ struct muonGlobalAlignment { // jacobian*covariances*jacobian^T covariances = ROOT::Math::Similarity(jacobian, covariances); - double cov[] = {covariances(0, 0), covariances(1, 0), covariances(1, 1), covariances(2, 0), covariances(2, 1), covariances(2, 2), covariances(3, 0), covariances(3, 1), covariances(3, 2), covariances(3, 3), covariances(4, 0), covariances(4, 1), covariances(4, 2), covariances(4, 3), covariances(4, 4)}; - double param[] = {fwdtrack.getX(), alpha1, fwdtrack.getY(), alpha3, alpha4}; + const std::array cov{covariances(0, 0), covariances(1, 0), covariances(1, 1), covariances(2, 0), covariances(2, 1), covariances(2, 2), covariances(3, 0), covariances(3, 1), covariances(3, 2), covariances(3, 3), covariances(4, 0), covariances(4, 1), covariances(4, 2), covariances(4, 3), covariances(4, 4)}; + const std::array param{fwdtrack.getX(), alpha1, fwdtrack.getY(), alpha3, alpha4}; - o2::mch::TrackParam convertedTrack(fwdtrack.getZ(), param, cov); - return o2::mch::TrackParam(convertedTrack); + o2::mch::TrackParam convertedTrack(fwdtrack.getZ(), param.data(), cov.data()); + return {convertedTrack}; } static o2::dataformats::GlobalFwdTrack MCHtoFwd(const o2::mch::TrackParam& mchParam) @@ -821,15 +832,13 @@ struct muonGlobalAlignment { o2::dataformats::GlobalFwdTrack convertedTrack; // Parameter conversion - double alpha1, alpha3, alpha4, x2, x3, x4; - - alpha1 = mchParam.getNonBendingSlope(); - alpha3 = mchParam.getBendingSlope(); - alpha4 = mchParam.getInverseBendingMomentum(); + double alpha1 = mchParam.getNonBendingSlope(); + double alpha3 = mchParam.getBendingSlope(); + double alpha4 = mchParam.getInverseBendingMomentum(); - x2 = TMath::ATan2(-alpha3, -alpha1); - x3 = -1. / TMath::Sqrt(alpha3 * alpha3 + alpha1 * alpha1); - x4 = alpha4 * -x3 * TMath::Sqrt(1 + alpha3 * alpha3); + double x2 = TMath::ATan2(-alpha3, -alpha1); + double x3 = -1. / TMath::Sqrt(alpha3 * alpha3 + alpha1 * alpha1); + double x4 = alpha4 * -x3 * TMath::Sqrt(1 + alpha3 * alpha3); auto K = alpha1 * alpha1 + alpha3 * alpha3; auto K32 = K * TMath::Sqrt(K); @@ -938,22 +947,22 @@ struct muonGlobalAlignment { auto itNextToNextParam = (itNextParam == track.end()) ? itNextParam : std::next(itNextParam); itStartingParam = track.rbegin(); - if (track.getNClusters() < 10) { + if (track.getNClusters() < NMchChambers) { removeTrack = true; break; - } else { - while (itNextToNextParam != track.end()) { - if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { - itStartingParam = std::make_reverse_iterator(++itNextParam); - break; - } - ++itNextToNextParam; + } + + while (itNextToNextParam != track.end()) { + if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { + itStartingParam = std::make_reverse_iterator(++itNextParam); + break; } + ++itNextToNextParam; } } if (!removeTrack) { - for (auto& param : track) { + for (auto& param : track) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) param.setParameters(param.getSmoothParameters()); param.setCovariances(param.getSmoothCovariances()); } @@ -968,12 +977,14 @@ struct muonGlobalAlignment { int nClustersCut) { // chi2 cut - if (mftTrack.chi2() > chi2Cut) + if (mftTrack.chi2() > chi2Cut) { return false; + } // number of clusters cut - if (mftTrack.nClusters() < nClustersCut) + if (mftTrack.nClusters() < nClustersCut) { return false; + } return true; } @@ -981,7 +992,7 @@ struct muonGlobalAlignment { template bool IsGoodMFT(const T& mftTrack) { - return IsGoodMFT(mftTrack, fTrackChi2MftUp, fTrackNClustMftLow); + return IsGoodMFT(mftTrack, cfgTrackChi2MftUp, cfgTrackNClustMftLow); } template @@ -1001,16 +1012,13 @@ struct muonGlobalAlignment { double p = mchTrackAtVertex.getP(); double pDCA = mchTrack.pDca(); - double sigmaPDCA = (thetaAbs < 3) ? sigmaPDCA23 : sigmaPDCA310; + double sigmaPDCA = (thetaAbs < ThetaAbsBoundaryDeg) ? sigmaPDCA23 : sigmaPDCA310; double nrp = nSigmaPDCA * relPRes * p; double pResEffect = sigmaPDCA / (1. - nrp / (1. + nrp)); - double slopeResEffect = 535. * slopeRes * p; + double slopeResEffect = SlopeResolutionZ * slopeRes * p; double sigmaPDCAWithRes = TMath::Sqrt(pResEffect * pResEffect + slopeResEffect * slopeResEffect); - if (pDCA > nSigmaPDCA * sigmaPDCAWithRes) { - return false; - } - return true; + return (pDCA <= nSigmaPDCA * sigmaPDCAWithRes); } template @@ -1022,11 +1030,12 @@ struct muonGlobalAlignment { std::array rAbsCut, double nSigmaPdcaCut) { - auto const& mchTrack = (static_cast(muonTrack.trackType()) <= 2) ? muonTrack.template matchMCHTrack_as() : muonTrack; + auto const& mchTrack = (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) ? muonTrack.template matchMCHTrack_as() : muonTrack; // chi2 cut - if (mchTrack.chi2() > chi2Cut) + if (mchTrack.chi2() > chi2Cut) { return false; + } // momentum cut if (mchTrack.p() < pCut) { @@ -1058,8 +1067,35 @@ struct muonGlobalAlignment { return true; } + template + bool isGoodGlobalMatching(const TMUON& muon, + double matchChi2Cut) + { + if (static_cast(muon.trackType()) > GlobalTrackTypeMax) + return false; + + // MFT-MCH match chi2 cut + if (muon.chi2MatchMCHMFT() > matchChi2Cut) + return false; + + return true; + } + + template + o2::track::TrackParCovFwd TrackToParCovFwd(const T& track) + { + double chi2 = track.chi2(); + SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + std::vector v1{0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0}; + SMatrix55 tcovs(v1.begin(), v1.end()); + o2::track::TrackParCovFwd trackparCov{track.z(), tpars, tcovs, chi2}; + return trackparCov; + } + template - o2::dataformats::GlobalFwdTrack FwdToTrackPar(const T& track) + o2::dataformats::GlobalFwdTrack TrackToGlobalFwd(const T& track) { double chi2 = track.chi2(); SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); @@ -1085,49 +1121,26 @@ struct muonGlobalAlignment { double xSlope = track.getNonBendingSlope(); double ySlope = track.getBendingSlope(); - double xSlopeCorrection = (y > 0) ? configMFTAlignmentCorrections.fMFTAlignmentCorrXSlopeTop : configMFTAlignmentCorrections.fMFTAlignmentCorrXSlopeBottom; + double xSlopeCorrection = (y > 0) ? configMFTAlignmentCorrections.cfgMFTAlignmentCorrXSlopeTop : configMFTAlignmentCorrections.cfgMFTAlignmentCorrXSlopeBottom; double xCorrection = xSlopeCorrection * z + - ((y > 0) ? configMFTAlignmentCorrections.fMFTAlignmentCorrXOffsetTop : configMFTAlignmentCorrections.fMFTAlignmentCorrXOffsetBottom); + ((y > 0) ? configMFTAlignmentCorrections.cfgMFTAlignmentCorrXOffsetTop : configMFTAlignmentCorrections.cfgMFTAlignmentCorrXOffsetBottom); track.setNonBendingCoor(x + xCorrection); track.setNonBendingSlope(xSlope + xSlopeCorrection); - double ySlopeCorrection = (y > 0) ? configMFTAlignmentCorrections.fMFTAlignmentCorrYSlopeTop : configMFTAlignmentCorrections.fMFTAlignmentCorrYSlopeBottom; + double ySlopeCorrection = (y > 0) ? configMFTAlignmentCorrections.cfgMFTAlignmentCorrYSlopeTop : configMFTAlignmentCorrections.cfgMFTAlignmentCorrYSlopeBottom; double yCorrection = ySlopeCorrection * z + - ((y > 0) ? configMFTAlignmentCorrections.fMFTAlignmentCorrYOffsetTop : configMFTAlignmentCorrections.fMFTAlignmentCorrYOffsetBottom); + ((y > 0) ? configMFTAlignmentCorrections.cfgMFTAlignmentCorrYOffsetTop : configMFTAlignmentCorrections.cfgMFTAlignmentCorrYOffsetBottom); track.setBendingCoor(y + yCorrection); track.setBendingSlope(ySlope + ySlopeCorrection); - /* - std::cout << std::format("[TOTO] MFT position: pos={:0.3f},{:0.3f}", x, y) << std::endl; - std::cout << std::format("[TOTO] MFT corrections: pos={:0.3f},{:0.3f} slope={:0.12f},{:0.12f} angle={:0.12f},{:0.12f}", - xCorrection, yCorrection, xSlopeCorrection, ySlopeCorrection, - std::atan2(xSlopeCorrection, 1), std::atan2(ySlopeCorrection, 1)) << std::endl; - */ - } - - void TransformMFT(o2::dataformats::GlobalFwdTrack& track) - { - auto mchTrack = FwdtoMCH(track); - - TransformMFTPar(mchTrack); - - auto transformedTrack = sExtrap.MCHtoFwd(mchTrack); - track.setParameters(transformedTrack.getParameters()); - track.setZ(transformedTrack.getZ()); - track.setCovariances(transformedTrack.getCovariances()); } void TransformMFT(o2::track::TrackParCovFwd& fwdtrack) { - o2::dataformats::GlobalFwdTrack track; - track.setParameters(fwdtrack.getParameters()); - track.setZ(fwdtrack.getZ()); - track.setCovariances(fwdtrack.getCovariances()); - - auto mchTrack = FwdtoMCH(track); + auto mchTrack = FwdtoMCH(fwdtrack); TransformMFTPar(mchTrack); - auto transformedTrack = sExtrap.MCHtoFwd(mchTrack); + auto transformedTrack = MCHtoFwd(mchTrack); fwdtrack.setParameters(transformedTrack.getParameters()); fwdtrack.setZ(transformedTrack.getZ()); fwdtrack.setCovariances(transformedTrack.getCovariances()); @@ -1136,8 +1149,8 @@ struct muonGlobalAlignment { template T UpdateTrackMomentum(const T& track, const double p, int sign) { - double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi()); - double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi()); + double px = p * std::sin(o2::constants::math::PIHalf - std::atan(track.tgl())) * std::cos(track.phi()); + double py = p * std::sin(o2::constants::math::PIHalf - std::atan(track.tgl())) * std::sin(track.phi()); double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt}; @@ -1157,8 +1170,8 @@ struct muonGlobalAlignment { template T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom) { - double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi()); - double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi()); + double px = track4mom.p() * std::sin(o2::constants::math::PIHalf - std::atan(track.tgl())) * std::cos(track.phi()); + double py = track4mom.p() * std::sin(o2::constants::math::PIHalf - std::atan(track.tgl())) * std::sin(track.phi()); double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); double sign = track4mom.getCharge(); @@ -1205,10 +1218,10 @@ struct muonGlobalAlignment { o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, z); } else if (z < absBack) { // extrapolation downstream of the absorber, correct for dipole longitudinal shift if needed - if (fDipoleZshift.value != 0) { - mchTrack.setZ(mchTrack.getZ() + fDipoleZshift.value); - o2::mch::TrackExtrap::extrapToZCov(mchTrack, z + fDipoleZshift.value); - mchTrack.setZ(mchTrack.getZ() - fDipoleZshift.value); + if (cfgDipoleZshift.value != 0) { + mchTrack.setZ(mchTrack.getZ() + cfgDipoleZshift.value); + o2::mch::TrackExtrap::extrapToZCov(mchTrack, z + cfgDipoleZshift.value); + mchTrack.setZ(mchTrack.getZ() - cfgDipoleZshift.value); } else { o2::mch::TrackExtrap::extrapToZCov(mchTrack, z); } @@ -1256,10 +1269,39 @@ struct muonGlobalAlignment { return PropagateMCH(track, z); } + o2::dataformats::GlobalFwdTrack PropagateMCHToVertex(const o2::track::TrackParCovFwd& muon, + const double vx, const double vy, const double vz, + const double covVx, const double covVy) + { + auto mchTrack = FwdtoMCH(muon); + + o2::mch::TrackExtrap::extrapToVertex(mchTrack, vx, vy, vz, covVx, covVy); + + auto proptrack = MCHtoFwd(mchTrack); + o2::dataformats::GlobalFwdTrack propmuon; + propmuon.setParameters(proptrack.getParameters()); + propmuon.setZ(proptrack.getZ()); + propmuon.setCovariances(proptrack.getCovariances()); + + return propmuon; + } + + template + o2::dataformats::GlobalFwdTrack PropagateMCHToVertex(const o2::track::TrackParCovFwd& muon, + const C& collision) + { + return PropagateMCHToVertex(muon, + collision.posX(), + collision.posY(), + collision.posZ(), + collision.covXX(), + collision.covYY()); + } + template o2::dataformats::GlobalFwdTrack PropagateMFT(const TMFT& mftTrack, float z) { - static double Bz = -10001; + // static double Bz = -10001; double chi2 = mftTrack.chi2(); SMatrix5 tpars = {mftTrack.x(), mftTrack.y(), mftTrack.phi(), mftTrack.tgl(), mftTrack.signed1Pt()}; std::vector v1{0, 0, 0, 0, 0, @@ -1277,12 +1319,12 @@ struct muonGlobalAlignment { // double centerZ[3] = {mftTrack.x() + propVec[0] / 2., // mftTrack.y() + propVec[1] / 2., // mftTrack.z() + propVec[2] / 2.}; - if (Bz < -10000) { - double centerZ[3] = {0, 0, (-45.f - 77.5f) / 2.f}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - Bz = field->getBz(centerZ); - } - fwdtrack.propagateToZ(z, Bz); + // if (Bz < -10000) { + // double centerZ[3] = {0, 0, (-45.f - 77.5f) / 2.f}; + // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + // Bz = field->getBz(centerZ); + //} + fwdtrack.propagateToZ(z, mBzAtMftCenter); propmuon.setParameters(fwdtrack.getParameters()); propmuon.setZ(fwdtrack.getZ()); @@ -1291,27 +1333,12 @@ struct muonGlobalAlignment { return propmuon; } - template - o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const C& collision, float zshift) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(o2::track::TrackParCovFwd mftTrack, + const C& collision, + float zshift) { - static double Bz = -10001; - double chi2 = mftTrack.chi2(); - double phiCorrDeg = 0; - double phiCorr = phiCorrDeg * TMath::Pi() / 180.f; - double tR = std::hypot(mftTrack.x(), mftTrack.y()); - double tphi = std::atan2(mftTrack.y(), mftTrack.x()); - double tx = std::cos(tphi + phiCorr) * tR; - double ty = std::sin(tphi + phiCorr) * tR; - SMatrix5 tpars = {tx, ty, mftTrack.phi() + phiCorr, mftTrack.tgl(), mftTrack.signed1Pt()}; - std::vector v1{0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0}; - SMatrix55 tcovs(v1.begin(), v1.end()); - o2::track::TrackParCovFwd fwdtrack{mftTrack.z(), tpars, tcovs, chi2}; - if (configMFTAlignmentCorrections.fEnableMFTAlignmentCorrections) { - TransformMFT(fwdtrack); - } - o2::dataformats::GlobalFwdTrack propmuon; + // static double Bz = -10001; // double propVec[3] = {}; // propVec[0] = collision.posX() - mftTrack.x(); @@ -1321,45 +1348,33 @@ struct muonGlobalAlignment { // double centerZ[3] = {mftTrack.x() + propVec[0] / 2., // mftTrack.y() + propVec[1] / 2., // mftTrack.z() + propVec[2] / 2.}; - if (Bz < -10000) { - double centerZ[3] = {0, 0, -45.f / 2.f}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - Bz = field->getBz(centerZ); - } - fwdtrack.propagateToZ(collision.posZ() - zshift, Bz); - - propmuon.setParameters(fwdtrack.getParameters()); - propmuon.setZ(fwdtrack.getZ()); - propmuon.setCovariances(fwdtrack.getCovariances()); - - return propmuon; + // if (Bz < -10000) { + // double centerZ[3] = {0, 0, -45.f / 2.f}; + // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + // Bz = field->getBz(centerZ); + // } + mftTrack.propagateToZ(collision.posZ() - zshift, mBzAtMftCenter); + + o2::dataformats::GlobalFwdTrack result; + result.setParameters(mftTrack.getParameters()); + result.setZ(mftTrack.getZ()); + result.setCovariances(mftTrack.getCovariances()); + + return result; } - template - o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const TMUON& mchTrack, const C& collision, float zshift) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(o2::track::TrackParCovFwd mftTrack, + const o2::track::TrackParCovFwd& mchTrack, + const C& collision, + float zshift) { - static double Bz = -10001; - double chi2 = mftTrack.chi2(); - double phiCorrDeg = 0; - double phiCorr = phiCorrDeg * TMath::Pi() / 180.f; - double tR = std::hypot(mftTrack.x(), mftTrack.y()); - double tphi = std::atan2(mftTrack.y(), mftTrack.x()); - double tx = std::cos(tphi + phiCorr) * tR; - double ty = std::sin(tphi + phiCorr) * tR; - SMatrix5 tpars = {tx, ty, mftTrack.phi() + phiCorr, mftTrack.tgl(), mftTrack.signed1Pt()}; - std::vector v1{0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0}; - SMatrix55 tcovs(v1.begin(), v1.end()); - o2::track::TrackParCovFwd fwdtrack{mftTrack.z(), tpars, tcovs, chi2}; - if (configMFTAlignmentCorrections.fEnableMFTAlignmentCorrections) { - TransformMFT(fwdtrack); - } + // static double Bz = -10001; // extrapolation with MCH tools - auto mchTrackAtMFT = FwdtoMCH(FwdToTrackPar(mchTrack)); - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z()); - UpdateTrackMomentum(fwdtrack, mchTrackAtMFT); + auto mchTrackAtMFT = FwdtoMCH(mchTrack); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.getZ()); + UpdateTrackMomentum(mftTrack, mchTrackAtMFT); // double propVec[3] = {}; // propVec[0] = collision.posX() - mftTrack.x(); @@ -1369,53 +1384,50 @@ struct muonGlobalAlignment { // double centerZ[3] = {mftTrack.x() + propVec[0] / 2., // mftTrack.y() + propVec[1] / 2., // mftTrack.z() + propVec[2] / 2.}; - if (Bz < -10000) { - double centerZ[3] = {0, 0, -45.f / 2.f}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - Bz = field->getBz(centerZ); - } - fwdtrack.propagateToZ(collision.posZ() - zshift, Bz); - - o2::dataformats::GlobalFwdTrack propmuon; - propmuon.setParameters(fwdtrack.getParameters()); - propmuon.setZ(fwdtrack.getZ()); - propmuon.setCovariances(fwdtrack.getCovariances()); - - return propmuon; + // if (Bz < -10000) { + // double centerZ[3] = {0, 0, -45.f / 2.f}; + // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + // Bz = field->getBz(centerZ); + // } + mftTrack.propagateToZ(collision.posZ() - zshift, mBzAtMftCenter); + + o2::dataformats::GlobalFwdTrack result; + result.setParameters(mftTrack.getParameters()); + result.setZ(mftTrack.getZ()); + result.setCovariances(mftTrack.getCovariances()); + + return result; } - template - o2::dataformats::GlobalFwdTrack PropagateMFTtoMCH(const TMFT& mftTrack, const o2::mch::TrackParam& mchTrackPar, const double z) + o2::dataformats::GlobalFwdTrack PropagateMFTtoMCH(o2::track::TrackParCovFwd mftTrackPar, + const o2::mch::TrackParam& mchTrackPar, + const double z) { // extrapolation with MCH tools auto mchTrackAtMFT = mchTrackPar; - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z()); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrackPar.getZ()); - auto mftTrackPar = FwdToTrackPar(mftTrack); - if (configMFTAlignmentCorrections.fEnableMFTAlignmentCorrections) { - TransformMFT(mftTrackPar); - } auto mftTrackProp = FwdtoMCH(mftTrackPar); UpdateTrackMomentum(mftTrackProp, mchTrackAtMFT); - if (z < -505.f) { - o2::mch::TrackExtrap::extrapToZ(mftTrackProp, -466.f); + if (z < AbsorberBackZ) { + o2::mch::TrackExtrap::extrapToZ(mftTrackProp, BransonPlaneZ); UpdateTrackMomentum(mftTrackProp, mchTrackPar); } - if (fDipoleZshift.value != 0) { + if (cfgDipoleZshift.value != 0) { // extrapolate to the back of the absorber, taking into account the dipole shift, // to avoid that the correction bring the track starting point back into the absorber - if (fDipoleZshift.value < 0) { - o2::mch::TrackExtrap::extrapToZ(mftTrackProp, -505.f); - } else if (fDipoleZshift.value > 0) { - o2::mch::TrackExtrap::extrapToZ(mftTrackProp, -505.f - fDipoleZshift.value); + if (cfgDipoleZshift.value < 0) { + o2::mch::TrackExtrap::extrapToZ(mftTrackProp, AbsorberBackZ); + } else if (cfgDipoleZshift.value > 0) { + o2::mch::TrackExtrap::extrapToZ(mftTrackProp, AbsorberBackZ - cfgDipoleZshift.value); } // shift the track starting point - mftTrackProp.setZ(mftTrackProp.getZ() + fDipoleZshift.value); + mftTrackProp.setZ(mftTrackProp.getZ() + cfgDipoleZshift.value); // extrapolate to the final z, corrected for the dipole shift - o2::mch::TrackExtrap::extrapToZ(mftTrackProp, z + fDipoleZshift.value); + o2::mch::TrackExtrap::extrapToZ(mftTrackProp, z + cfgDipoleZshift.value); // remove the shift from the extrapolated track - mftTrackProp.setZ(mftTrackProp.getZ() - fDipoleZshift.value); + mftTrackProp.setZ(mftTrackProp.getZ() - cfgDipoleZshift.value); } else { o2::mch::TrackExtrap::extrapToZ(mftTrackProp, z); } @@ -1423,82 +1435,385 @@ struct muonGlobalAlignment { return MCHtoFwd(mftTrackProp); } - void FillDCAPlots(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - MyMFTs const& mftTracks, - const std::map& collisionInfos) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToVertex(const o2::track::TrackParCovFwd& mftTrackPar, + const o2::track::TrackParCovFwd& mchTrackPar, + const C& collision) { - // outer loop over collisions - for (auto& [collisionIndex, collisionInfo] : collisionInfos) { - auto const& collision = collisions.rawIteratorAt(collisionIndex); - const auto& bc = bcs.rawIteratorAt(collision.bcId()); + // extrapolation with MCH tools + auto mchTrackAtMFT = FwdtoMCH(mchTrackPar); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrackPar.getZ()); - // remove TF/ROF borders and ambiguous collisions - if (!bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || - !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) - continue; + auto mftTrackProp = FwdtoMCH(mftTrackPar); - registry.get(HIST("vertex_y_vs_x"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vertex_z"))->Fill(collision.posZ()); + // update global track momentum from the MCH track + double pRatio = mftTrackProp.p() / mchTrackAtMFT.p(); + double newInvBendMom = mftTrackProp.getInverseBendingMomentum() * pRatio; + mftTrackProp.setInverseBendingMomentum(newInvBendMom); + mftTrackProp.setCharge(mchTrackAtMFT.getCharge()); - if (fEnableVertexShiftAnalysis || fEnableMftDcaAnalysis) { - registry.get(HIST("DCA/MFT/nTracksMFT"))->Fill(collisionInfo.mftTracks.size()); - } + o2::mch::TrackExtrap::extrapToVertex(mftTrackProp, + collision.posX(), + collision.posY(), + collision.posZ(), + collision.covXX(), + collision.covYY()); - if (fEnableVertexShiftAnalysis || fEnableMftDcaAnalysis) { - // loop over MFT tracks - auto mftTrackIds = collisionInfo.mftTracks; - if (fMftTracksMultiplicityMax > 0 && mftTrackIds.size() > fMftTracksMultiplicityMax) { - auto rng = std::default_random_engine{}; - std::shuffle(std::begin(mftTrackIds), std::end(mftTrackIds), rng); - mftTrackIds.resize(fMftTracksMultiplicityMax); - } + return MCHtoFwd(mftTrackProp); + } - for (auto mftIndex : mftTrackIds) { - auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + void getMuonPairs(const CollisionInfo& collisionInfo, + std::vector& muonPairs) + { + // outer loop over muon tracks + for (const auto& mchIndex1 : collisionInfo.mchTracks) { + // inner loop over muon tracks + for (const auto& mchIndex2 : collisionInfo.mchTracks) { + // avoid double-counting of muon pairs + if (mchIndex2 <= mchIndex1) { + continue; + } - if (mftTrack.isCA()) { - continue; - } + muonPairs.emplace_back(mchIndex1, mchIndex2); + } + } + } - bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5); - if (!isGoodMFT) - continue; + ROOT::Math::PxPyPzMVector getMuMu4Momentum(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + { + ROOT::Math::PxPyPzMVector muon1{ + track1.getPx(), + track1.getPy(), + track1.getPz(), + o2::constants::physics::MassMuon}; + + ROOT::Math::PxPyPzMVector muon2{ + track2.getPx(), + track2.getPy(), + track2.getPz(), + o2::constants::physics::MassMuon}; + + return muon1 + muon2; + } - auto mftTrackAtDCA = PropagateMFTToDCA(mftTrack, collision, fVertexZshift); - double dcax = mftTrackAtDCA.getX() - collision.posX(); - double dcay = mftTrackAtDCA.getY() - collision.posY(); - double phi = mftTrack.phi() * 180 / TMath::Pi(); - int mftNclusters = mftTrack.nClusters(); - double chi2NDF = static_cast(mftNclusters) * 2 - 5; + double getMuMuAngle(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + { + ROOT::Math::XYZVector muon1{ + track1.getPx(), + track1.getPy(), + track1.getPz()}; - const int nMftLayers = 10; - std::array firedLayers; - for (int layer = 0; layer < nMftLayers; layer++) { - if ((mftTrack.mftClusterSizesAndTrackFlags() >> (layer * 6)) & 0x3F) { - firedLayers[layer] = true; - } else { - firedLayers[layer] = false; - } - } + ROOT::Math::XYZVector muon2{ + track2.getPx(), + track2.getPy(), + track2.getPz()}; - if (fEnableMftDcaAnalysis) { - const int nMftLayers = 10; - if (fEnableMftDcaExtraPlots) { - for (int i = 0; i < nMftLayers; i++) { - if (firedLayers[i]) { - registry.get(HIST("DCA/MFT/trackChi2"))->Fill(mftTrack.chi2() / chi2NDF, mftTrack.x(), mftTrack.y(), mftNclusters, i); - } - } - } + return std::acos(muon1.Unit().Dot(muon2.Unit())); + } - if (mftTrack.chi2() <= fTrackChi2MftUp) { - registry.get(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay); - registry.get(HIST("DCA/MFT/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); - registry.get(HIST("DCA/MFT/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + double getMuMuInvariantMass(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + { + return getMuMu4Momentum(track1, track2).M(); + } - if (cfgProduceMFTTable) { + template + bool MchRealignTrack(const TMUON& mchTrack, const TCLUS& clusters, TrackRealigned& convertedTrack, bool applyCorrections) + { + auto mchTrackPar = FwdtoMCH(TrackToGlobalFwd(mchTrack)); + + // loop over attached clusters + int clIndex = -1; + auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id + for (auto const& cluster : clustersSliced) { + clIndex += 1; + + int deId = cluster.deId(); + int chamber = deId / 100 - 1; + if (chamber < 0 || chamber >= NMchChambers) { + continue; + } + + math_utils::Point3D local; + math_utils::Point3D master; + + master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); + + if (configRealign.cfgEnableMCHRealign) { + // Transformation from reference geometry frame to new geometry frame + transformRef[cluster.deId()].MasterToLocal(master, local); + transformNew[cluster.deId()].LocalToMaster(local, master); + } + + // shift the clusters to correct the longitudinal shift of the dipole + if (cfgDipoleZshift.value != 0) { + master.SetZ(master.z() + cfgDipoleZshift.value); + } + + if (applyCorrections) { + auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId()); + if (correctionsIt != mMchAlignmentCorrections.end()) { + const auto& corrections = correctionsIt->second; + master.SetX(master.x() + corrections.x); + master.SetY(master.y() + corrections.y); + master.SetZ(master.z() + corrections.z); + } + } + + // realigned MCH cluster + auto clusterMCH = new mch::Cluster(); + clusterMCH->x = master.x(); + clusterMCH->y = master.y(); + clusterMCH->z = master.z(); + + uint32_t ClUId = mch::Cluster::buildUniqueId(static_cast(cluster.deId() / 100) - 1, cluster.deId(), clIndex); + clusterMCH->uid = ClUId; + clusterMCH->ex = cluster.isGoodX() ? 0.2 : 10.0; + clusterMCH->ey = cluster.isGoodY() ? 0.2 : 10.0; + + // Add transformed cluster into temporary variable + convertedTrack.createParamAtCluster(*clusterMCH); + } + + bool removable{false}; + // Refit the re-aligned track + if (convertedTrack.getNClusters() != 0) { + removable = RemoveTrack(convertedTrack); + } else { + LOGF(fatal, "Muon track %d has no associated clusters.", mchTrack.globalIndex()); + } + + // subtract the longitudinal shift of the dipole from the track z + if (cfgDipoleZshift.value != 0) { + auto& trackParam = *(convertedTrack.begin()); + trackParam.setZ(trackParam.getZ() - cfgDipoleZshift.value); + } + + return !removable; + } + + template + void InitCollisions(COLL const& collisions, + BC const& bcs, + TMUON const& muonTracks, + aod::FwdTrkCls const& clusters, + std::map& collisionInfos) + { + mMchTrackPars.clear(); + mMchTrackParsNew.clear(); + + // fill collision information for global muon tracks (MFT-MCH-MID matches) + for (const auto& muonTrack : muonTracks) { + if (!muonTrack.has_collision()) { + continue; + } + + auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); + + if (cfgRequireGoodRCT && !rctChecker(collision)) { + continue; + } + + uint64_t collisionIndex = collision.globalIndex(); + + auto bc = bcs.rawIteratorAt(collision.bcId()); + + auto& collisionInfo = collisionInfos[collisionIndex]; + collisionInfo.bc = bc.globalBC(); + collisionInfo.zVertex = collision.posZ(); + + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { + // standalone MCH or MCH-MID tracks + uint64_t mchTrackIndex = muonTrack.globalIndex(); + collisionInfo.mchTracks.push_back(mchTrackIndex); + + // initialize the original MCH track parameters + mMchTrackPars.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); + + // refit MCH track if requested + if (configRealign.cfgEnableMCHRefit || configRealign.cfgEnableMCHRealign) { + TrackRealigned convertedTrack; + bool convertedTrackOk = MchRealignTrack(muonTrack, clusters, convertedTrack, !mMchAlignmentCorrections.empty()); + + // Get the re-aligned track parameters: track param at the first cluster + mch::TrackParam trackParam = mch::TrackParam(convertedTrack.first()); + + auto mchTrackParIt = mMchTrackParsNew.try_emplace(mchTrackIndex, TrackParExt(MCHtoFwd(trackParam), convertedTrack.getNClusters())); + if (mchTrackParIt.second) { + // the insertion succeeded + mchTrackParIt.first->second.setTrackChi2(trackParam.getTrackChi2() / convertedTrack.getNDF()); + if (!convertedTrackOk) { + mchTrackParIt.first->second.setRemovable(); + } + } + } else { + // initialize the new MCH track parameters with the original ones, without refitting + mMchTrackParsNew.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); + } + } else { + // global muon tracks (MFT-MCH or MFT-MCH-MID) + uint64_t muonTrackIndex = muonTrack.globalIndex(); + auto const& mchTrack = muonTrack.template matchMCHTrack_as(); + uint64_t mchTrackIndex = mchTrack.globalIndex(); + + // check if a vector of global muon candidates is already available for the current MCH index + // if not, initialize a new one and add the current global muon track + // bool globalMuonTrackFound = false; + auto matchingCandidateIterator = collisionInfo.globalMuonTracks.find(mchTrackIndex); + if (matchingCandidateIterator != collisionInfo.globalMuonTracks.end()) { + matchingCandidateIterator->second.push_back(muonTrackIndex); + // globalMuonTrackFound = true; + } else { + collisionInfo.globalMuonTracks[mchTrackIndex].push_back(muonTrackIndex); + } + } + } + + // sort the vectors of matching candidates in ascending order based on the matching chi2 value + auto compareChi2 = [&muonTracks](uint64_t trackIndex1, uint64_t trackIndex2) -> bool { + auto const& track1 = muonTracks.rawIteratorAt(trackIndex1); + auto const& track2 = muonTracks.rawIteratorAt(trackIndex2); + + return (track1.chi2MatchMCHMFT() < track2.chi2MatchMCHMFT()); + }; + + for (auto& [collisionIndex, collisionInfo] : collisionInfos) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + for (auto& [mchIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareChi2); + } + } + } + + void InitCollisions(MyEvents const& collisions, + MyBCs const& bcs, + MyMuonsWithCov const& muonTracks, + aod::FwdTrkCls const& clusters, + MyMFTs const& mftTracks, + std::map& collisionInfos) + { + InitCollisions(collisions, bcs, muonTracks, clusters, collisionInfos); + + mMftTrackPars.clear(); + mMftTrackParsNew.clear(); + + // fill collision information for MFT standalone tracks + for (const auto& mftTrack : mftTracks) { + if (!mftTrack.has_collision()) { + continue; + } + + auto collision = collisions.rawIteratorAt(mftTrack.collisionId()); + uint64_t collisionIndex = collision.globalIndex(); + + auto bc = bcs.rawIteratorAt(collision.bcId()); + + uint64_t mftTrackIndex = mftTrack.globalIndex(); + + auto& collisionInfo = collisionInfos[collisionIndex]; + collisionInfo.bc = bc.globalBC(); + collisionInfo.zVertex = collision.posZ(); + + collisionInfo.mftTracks.push_back(mftTrackIndex); + + // initialize the original MFT track parameters + auto mftTrackFwd = TrackToParCovFwd(mftTrack); + mMftTrackPars.try_emplace(mftTrackIndex, TrackParExt(mftTrackFwd, mftTrack.nClusters())); + + // initialize the corrected MFT track parameters, if requested + if (configMFTAlignmentCorrections.cfgEnableMFTAlignmentCorrections) { + TransformMFT(mftTrackFwd); + mMftTrackParsNew.try_emplace(mftTrackIndex, TrackParExt(mftTrackFwd, mftTrack.nClusters())); + } else { + // initialize the new MFT track parameters with the original ones, without corrections + mMftTrackParsNew.try_emplace(mftTrackIndex, TrackParExt(mftTrackFwd, mftTrack.nClusters())); + } + } + } + + void FillMftPlots(MyEvents const& collisions, + MyBCs const& bcs, + MyMuonsWithCov const& muonTracks, + MyMFTs const& mftTracks, + const std::map& collisionInfos) + { + // outer loop over collisions + for (const auto& [collisionIndex, collisionInfo] : collisionInfos) { + auto const& collision = collisions.rawIteratorAt(collisionIndex); + const auto& bc = bcs.rawIteratorAt(collision.bcId()); + + // remove TF/ROF borders and ambiguous collisions + if (!bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || + !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + continue; + } + + registry.get(HIST("vertex_y_vs_x"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("vertex_z"))->Fill(collision.posZ()); + + if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { + registry.get(HIST("DCA/MFT/nTracksMFT"))->Fill(collisionInfo.mftTracks.size()); + } + + if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { + // loop over MFT tracks + auto mftTrackIds = collisionInfo.mftTracks; + if (cfgMftTracksMultiplicityMax > 0 && mftTrackIds.size() > cfgMftTracksMultiplicityMax) { + auto rng = std::default_random_engine{}; + std::shuffle(std::begin(mftTrackIds), std::end(mftTrackIds), rng); + mftTrackIds.resize(cfgMftTracksMultiplicityMax); + } + + for (const auto& mftIndex : mftTrackIds) { + auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + + if (mftTrack.isCA()) { + continue; + } + + bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5); + if (!isGoodMFT) { + continue; + } + + // get the pre-stored MFT track parameters after corrections + // if MFT corrections are not enabled, the original MFT track parameters are retrieved + const auto mftTrackParIt = mMftTrackParsNew.find(mftIndex); + if (mftTrackParIt == mMftTrackParsNew.end()) { + continue; + } + const auto& mftTrackPar = mftTrackParIt->second; + + auto mftTrackAtDCA = PropagateMFTToDCA(mftTrackPar, collision, cfgVertexZshift); + double dcax = mftTrackAtDCA.getX() - collision.posX(); + double dcay = mftTrackAtDCA.getY() - collision.posY(); + double phi = mftTrack.phi() * o2::constants::math::Rad2Deg; + int mftNclusters = mftTrack.nClusters(); + double chi2NDF = static_cast(mftNclusters) * 2 - 5; + + const int nMftLayers = 10; + std::array firedLayers{false}; + for (int layer = 0; layer < nMftLayers; layer++) { + if (((mftTrack.mftClusterSizesAndTrackFlags() >> (layer * 6)) & 0x3F) != 0) { + firedLayers[layer] = true; + } else { + firedLayers[layer] = false; + } + } + + if (cfgEnableMftDcaAnalysis) { + if (cfgEnableMftDcaExtraPlots) { + for (int i = 0; i < nMftLayers; i++) { + if (firedLayers[i]) { + registry.get(HIST("DCA/MFT/trackChi2"))->Fill(mftTrack.chi2() / chi2NDF, mftTrack.x(), mftTrack.y(), mftNclusters, i); + } + } + } + + if (mftTrack.chi2() <= cfgTrackChi2MftUp) { + registry.get(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay); + registry.get(HIST("DCA/MFT/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + registry.get(HIST("DCA/MFT/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + + if (cfgProduceMFTTable) { mftTable(collision.posX(), collision.posY(), collision.posZ(), mftTrack.signed1Pt(), mftTrack.tgl(), mftTrack.phi(), dcax, dcay, mftTrackAtDCA.getSigma2X(), mftTrackAtDCA.getSigma2Y(), mftTrackAtDCA.getSigmaXY(), @@ -1506,8 +1821,9 @@ struct muonGlobalAlignment { mftTrack.x(), mftTrack.y(), mftTrack.z()); } - if (fEnableMftDcaExtraPlots) { - if (mftNclusters >= 6) { + if (cfgEnableMftDcaExtraPlots) { + static constexpr int nMftClustersMin = 6; + if (mftNclusters >= nMftClustersMin) { for (int i = 0; i < nMftLayers; i++) { auto mftTrackAtLayer = PropagateMFT(mftTrack, o2::mft::constants::mft::LayerZCoordinate()[i]); std::get>(mMftTrackEffDen[i])->Fill(mftTrackAtLayer.getX(), mftTrackAtLayer.getY()); @@ -1526,13 +1842,15 @@ struct muonGlobalAlignment { } } - if (fEnableVertexShiftAnalysis) { - if (mftTrack.chi2() <= fTrackChi2MftUp && std::fabs(collision.posZ()) < 1.f && mftNclusters >= 6) { - float zshift[21] = {// in millimeters - -5.0, -4.5, -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, - 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}; - for (int zi = 0; zi < 21; zi++) { - auto mftTrackAtDCAshifted = PropagateMFTToDCA(mftTrack, collision, zshift[zi] / 10.f); + if (cfgEnableVertexShiftAnalysis) { + static constexpr int nMftClustersMin = 6; + if (mftTrack.chi2() <= cfgTrackChi2MftUp && std::fabs(collision.posZ()) < 1.f && mftNclusters >= nMftClustersMin) { + static constexpr int nPoints = 21; + const std::array zshift{// in millimeters + -5.0, -4.5, -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, + 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}; + for (int zi = 0; zi < nPoints; zi++) { + auto mftTrackAtDCAshifted = PropagateMFTToDCA(mftTrackPar, collision, zshift[zi] / 10.f); double dcaxShifted = mftTrackAtDCAshifted.getX() - collision.posX(); double dcayShifted = mftTrackAtDCAshifted.getY() - collision.posY(); registry.get(HIST("DCA/MFT/DCA_x_vs_phi_vs_zshift"))->Fill(zshift[zi], phi, dcaxShifted); @@ -1551,14 +1869,17 @@ struct muonGlobalAlignment { } } - if (fEnableGlobalFwdDcaAnalysis) { + if (cfgEnableGlobalFwdDcaAnalysis) { // loop over global muon tracks - for (auto& [muonIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { + for (const auto& [muonIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0]); const auto& mchTrack = muonTrack.template matchMCHTrack_as(); const auto& mftTrack = muonTrack.template matchMFTTrack_as(); - if (muonTrack.chi2MatchMCHMFT() < 50) { + auto mchIndex = mchTrack.globalIndex(); + auto mftIndex = mftTrack.globalIndex(); + + if (muonTrack.chi2MatchMCHMFT() < cfgMftDcaMatchChi2Up.value) { continue; } @@ -1566,7 +1887,7 @@ struct muonGlobalAlignment { auto const& muonTrack2 = muonTracks.rawIteratorAt(globalTracksVector[1]); double dchi2 = muonTrack2.chi2MatchMCHMFT() - muonTrack.chi2MatchMCHMFT(); - if (dchi2 < 50) { + if (dchi2 < cfgMftDcaMatchChi2Up.value) { continue; } } @@ -1575,21 +1896,37 @@ struct muonGlobalAlignment { continue; } - bool isGoodMFT = IsGoodMFT(mftTrack, fTrackChi2MftUp, 5); + bool isGoodMFT = IsGoodMFT(mftTrack, cfgTrackChi2MftUp, 5); if (!isGoodMFT) { continue; } - bool isGoodMuon = IsGoodMuon(mchTrack, collision, fTrackChi2MchUp, 0.f, fPtMchLow, {fEtaMftLow, fEtaMftUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp); - if (!isGoodMuon) + bool isGoodMuon = IsGoodMuon(mchTrack, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + if (!isGoodMuon) { continue; + } + + // get the pre-stored MFT track parameters after corrections + // if MFT corrections are not enabled, the original MFT track parameters are retrieved + const auto mftTrackParIt = mMftTrackParsNew.find(mftIndex); + if (mftTrackParIt == mMftTrackParsNew.end()) { + continue; + } + const auto& mftTrackPar = mftTrackParIt->second; + + // get the pre-stored MCH track parameters + const auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; + } + const auto& mchTrackPar = mchTrackParIt->second; - auto mftTrackAtDCA = PropagateMFTToDCA(mftTrack, mchTrack, collision, fVertexZshift); + auto mftTrackAtDCA = PropagateMFTToDCA(mftTrackPar, mchTrackPar, collision, cfgVertexZshift); double dcax = mftTrackAtDCA.getX() - collision.posX(); double dcay = mftTrackAtDCA.getY() - collision.posY(); int mftNclusters = mftTrack.nClusters(); - if (mftTrack.chi2() <= fTrackChi2MftUp) { + if (mftTrack.chi2() <= cfgTrackChi2MftUp) { registry.get(HIST("DCA/GlobalFwd/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); registry.get(HIST("DCA/GlobalFwd/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); } @@ -1598,157 +1935,105 @@ struct muonGlobalAlignment { } } - template - bool MchRealignTrack(const TMUON& mchTrack, const TCLUS& clusters, TrackRealigned& convertedTrack, bool applyCorrections) - { - // loop over attached clusters - int clIndex = -1; - auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id - for (auto const& cluster : clustersSliced) { - clIndex += 1; - - int deId = cluster.deId(); - int chamber = deId / 100 - 1; - if (chamber < 0 || chamber > 9) { - continue; - } - - math_utils::Point3D local; - math_utils::Point3D master; - - master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); - - if (configRealign.fEnableMCHRealign) { - // Transformation from reference geometry frame to new geometry frame - transformRef[cluster.deId()].MasterToLocal(master, local); - transformNew[cluster.deId()].LocalToMaster(local, master); - } - - // shift the clusters to correct the longitudinal shift of the dipole - if (fDipoleZshift.value != 0) { - master.SetZ(master.z() + fDipoleZshift.value); - } - - if (applyCorrections) { - auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId()); - if (correctionsIt != mMchAlignmentCorrections.end()) { - const auto& corrections = correctionsIt->second; - master.SetX(master.x() + corrections.x); - master.SetY(master.y() + corrections.y); - master.SetZ(master.z() + corrections.z); - } - } - - // realigned MCH cluster - mch::Cluster* clusterMCH = new mch::Cluster(); - clusterMCH->x = master.x(); - clusterMCH->y = master.y(); - clusterMCH->z = master.z(); - - uint32_t ClUId = mch::Cluster::buildUniqueId(static_cast(cluster.deId() / 100) - 1, cluster.deId(), clIndex); - clusterMCH->uid = ClUId; - clusterMCH->ex = cluster.isGoodX() ? 0.2 : 10.0; - clusterMCH->ey = cluster.isGoodY() ? 0.2 : 10.0; - - // Add transformed cluster into temporary variable - convertedTrack.createParamAtCluster(*clusterMCH); - } - - bool removable{false}; - // Refit the re-aligned track - if (convertedTrack.getNClusters() != 0) { - removable = RemoveTrack(convertedTrack); - } else { - LOGF(fatal, "Muon track %d has no associated clusters.", mchTrack.globalIndex()); - } - - // subtract the longitudinal shift of the dipole from the track z - if (fDipoleZshift.value != 0) { - auto& trackParam = *(convertedTrack.begin()); - trackParam.setZ(trackParam.getZ() - fDipoleZshift.value); - } - - return !removable; - } - - void FillResidualsPlots(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - aod::FwdTrkCls const& clusters, - const std::map& collisionInfos) + void FillMchPlots(MyEvents const& collisions, + MyBCs const& bcs, + MyMuonsWithCov const& muonTracks, + aod::FwdTrkCls const& clusters, + const std::map& collisionInfos) { - if (!fEnableMftMchResidualsAnalysis && !fEnableMftMchMatchingAnalysis) { + if (!cfgEnableMftMchResidualsAnalysis && !cfgEnableMftMchMatchingAnalysis) { return; } // loop over collisions - for (auto& [collisionIndex, collisionInfo] : collisionInfos) { + for (const auto& [collisionIndex, collisionInfo] : collisionInfos) { auto const& collision = collisions.rawIteratorAt(collisionIndex); const auto& bc = bcs.rawIteratorAt(collision.bcId()); // remove TF/ROF borders and ambiguous collisions if (!bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || - !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) + !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { continue; + } // loop over global muon tracks - for (auto& [muonIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { + for (const auto& [muonIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0]); const auto& mchTrack = muonTrack.template matchMCHTrack_as(); const auto& mftTrack = muonTrack.template matchMFTTrack_as(); // int quadrant = GetQuadrant(mchTrack); int quadrant = GetQuadrant(mftTrack); + //int quadrantMCH = GetQuadrant(mchTrack); int posNeg = (mchTrack.sign() >= 0) ? 0 : 1; - bool isGoodMuon = IsGoodMuon(mchTrack, collision, fTrackChi2MchUp, fMftMchResidualsPLow, fMftMchResidualsPtLow, {fEtaMftLow, fEtaMftUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp); - if (!isGoodMuon) + auto mchIndex = mchTrack.globalIndex(); + auto mftIndex = mftTrack.globalIndex(); + + bool isGoodMuon = IsGoodMuon(mchTrack, collision, cfgTrackChi2MchUp, cfgMftMchResidualsPLow, cfgMftMchResidualsPtLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + if (!isGoodMuon) { + continue; + } + + bool isGoodMFT = IsGoodMFT(mftTrack, cfgTrackChi2MftUp, cfgTrackNClustMftLow); + if (!isGoodMFT) { continue; + } - bool isGoodMFT = IsGoodMFT(mftTrack, fTrackChi2MftUp, fTrackNClustMftLow); - if (!isGoodMFT) + // get the pre-stored MFT track parameters + const auto mftTrackParIt = mMftTrackPars.find(mftIndex); + if (mftTrackParIt == mMftTrackPars.end()) { continue; + } + const auto& mftTrackPar = mftTrackParIt->second; - double matchChi2 = muonTrack.chi2MatchMCHMFT() / 5.f; - if (matchChi2 > 10.f) + // get the pre-stored MFT track parameters after corrections + // if MFT corrections are not enabled, the original MFT track parameters are retrieved + const auto mftTrackParNewIt = mMftTrackParsNew.find(mftIndex); + if (mftTrackParNewIt == mMftTrackParsNew.end()) { continue; + } + const auto& mftTrackParNew = mftTrackParNewIt->second; - // refit MCH track if enabled - TrackRealigned convertedTrack; - bool convertedTrackOk = false; - if (configRealign.fEnableMCHRealign) { - convertedTrackOk = MchRealignTrack(mchTrack, clusters, convertedTrack, false); + // get the pre-stored MCH track parameters + const auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; } + const auto& mchTrackPar = mchTrackParIt->second; - // apply alignment corrections if available - TrackRealigned convertedTrackWithCorr; - bool convertedTrackWithCorrOk = false; - if (!mMchAlignmentCorrections.empty()) { - convertedTrackWithCorrOk = MchRealignTrack(mchTrack, clusters, convertedTrackWithCorr, true); + // get the pre-stored MCH track parameters after refit + const auto mchTrackParNewIt = mMchTrackParsNew.find(mchIndex); + if (mchTrackParNewIt == mMchTrackParsNew.end()) { + continue; } + const auto& mchTrackParNew = mchTrackParNewIt->second; - if (fEnableMftMchResidualsAnalysis) { + double matchChi2 = muonTrack.chi2MatchMCHMFT() / 5.f; + + // Residuals analysis between MFT tracks and MCH clusters + if (cfgEnableMftMchResidualsAnalysis && (matchChi2 <= 10.f)) { // loop over attached clusters auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id for (auto const& cluster : clustersSliced) { int deId = cluster.deId(); int chamber = deId / 100 - 1; - if (chamber < 0 || chamber > 9) + if (chamber < 0 || chamber >= NMchChambers) { continue; + } int deIndex = getDEindex(deId); math_utils::Point3D local; - math_utils::Point3D master; - math_utils::Point3D masterWithCorr; + math_utils::Point3D master; // original cluster position + math_utils::Point3D masterRealign; // cluster position after realignment master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); - masterWithCorr.SetXYZ(cluster.x(), cluster.y(), cluster.z()); + masterRealign.SetXYZ(cluster.x(), cluster.y(), cluster.z()); // apply realignment to MCH cluster - if (configRealign.fEnableMCHRealign) { + if (configRealign.cfgEnableMCHRealign) { // Transformation from reference geometry frame to new geometry frame transformRef[cluster.deId()].MasterToLocal(master, local); - transformNew[cluster.deId()].LocalToMaster(local, master); - transformNew[cluster.deId()].LocalToMaster(local, masterWithCorr); + transformNew[cluster.deId()].LocalToMaster(local, masterRealign); } // apply alignment corrections to MCH cluster (if available) @@ -1756,39 +2041,35 @@ struct muonGlobalAlignment { auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId()); if (correctionsIt != mMchAlignmentCorrections.end()) { const auto& corrections = correctionsIt->second; - masterWithCorr.SetX(masterWithCorr.x() + corrections.x); - masterWithCorr.SetY(masterWithCorr.y() + corrections.y); - masterWithCorr.SetZ(masterWithCorr.z() + corrections.z); + masterRealign.SetX(masterRealign.x() + corrections.x); + masterRealign.SetY(masterRealign.y() + corrections.y); + masterRealign.SetZ(masterRealign.z() + corrections.z); } } - // MFT-MCH residuals (MCH cluster is realigned if enabled) - // if the realignment is enabled and successful, the MFT track is extrpolated - // by taking the momentum from the MCH track refitted with the new alignment - if (!configRealign.fEnableMCHRealign || convertedTrackOk) { - auto mftTrackAtCluster = configRealign.fEnableMCHRealign ? PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrack.first()), master.z()) : PropagateMFTtoMCH(mftTrack, FwdtoMCH(FwdToTrackPar(mchTrack)), master.z()); - auto mftTrackParamAtCluster = FwdtoMCH(mftTrackAtCluster); + // MFT-MCH residuals from original alignment + const auto mftTrackAtCluster = PropagateMFTtoMCH(mftTrackPar, FwdtoMCH(mchTrackPar), master.z()); + const auto mftTrackParamAtCluster = FwdtoMCH(mftTrackAtCluster); - std::array xPos{master.x(), mftTrackAtCluster.getX()}; - std::array yPos{master.y(), mftTrackAtCluster.getY()}; + const std::array xPos{master.x(), mftTrackAtCluster.getX()}; + const std::array yPos{master.y(), mftTrackAtCluster.getY()}; - registry.get(HIST("residuals/dx_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]); - registry.get(HIST("residuals/dy_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]); + registry.get(HIST("residuals/dx_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]); + registry.get(HIST("residuals/dy_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]); - registry.get(HIST("residuals/dx_vs_de"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getNonBendingSlope()); - registry.get(HIST("residuals/dy_vs_de"))->Fill(yPos[0] - yPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getBendingSlope()); - } + registry.get(HIST("residuals/dx_vs_de"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getNonBendingSlope()); + registry.get(HIST("residuals/dy_vs_de"))->Fill(yPos[0] - yPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getBendingSlope()); // MFT-MCH residuals with realigned and/or corrected MCH clusters // if the alignment corrections are available and the refitting is successful, the MFT track is extrpolated // by taking the momentum from the MCH track refitted with the alignment corrections and the new // alignment (if realignment is enabled) - if (convertedTrackWithCorrOk) { - auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrackWithCorr.first()), masterWithCorr.z()); - auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr); + if (!mchTrackParNew.isRemovable()) { + const auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrackParNew, FwdtoMCH(mchTrackParNew), masterRealign.z()); + const auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr); - std::array xPos{masterWithCorr.x(), mftTrackAtClusterWithCorr.getX()}; - std::array yPos{masterWithCorr.y(), mftTrackAtClusterWithCorr.getY()}; + const std::array xPos{masterRealign.x(), mftTrackAtClusterWithCorr.getX()}; + const std::array yPos{masterRealign.y(), mftTrackAtClusterWithCorr.getY()}; registry.get(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]); registry.get(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]); @@ -1798,75 +2079,259 @@ struct muonGlobalAlignment { } } - if (!configRealign.fEnableMCHRealign || convertedTrackOk) { - auto mchTrackAtDCA = configRealign.fEnableMCHRealign ? PropagateMCHRealigned(convertedTrack, collision.posZ()) : PropagateMCH(mchTrack, collision.posZ()); - auto dcax = mchTrackAtDCA.getX() - collision.posX(); - auto dcay = mchTrackAtDCA.getY() - collision.posY(); - - registry.get(HIST("DCA/MCH/DCA_y_vs_x"))->Fill(dcax, dcay); - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcay); - - if (fEnableMftMchResidualsExtraPlots) { - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcay); - auto mchTrackAtMFT = configRealign.fEnableMCHRealign ? PropagateMCHRealigned(convertedTrack, mftTrack.z()) : PropagateMCH(mchTrack, mftTrack.z()); - double deltaPhi = mchTrackAtMFT.getPhi() - mftTrack.phi(); - registry.get(HIST("residuals/dphi_at_mft"))->Fill(deltaPhi, mftTrack.x(), mftTrack.y(), posNeg, mchTrackAtMFT.getP()); - } + const auto mchTrackAtDCA = PropagateMCHParam(FwdtoMCH(mchTrackPar), collision.posZ()); + const auto dcax = mchTrackAtDCA.getX() - collision.posX(); + const auto dcay = mchTrackAtDCA.getY() - collision.posY(); + + registry.get(HIST("DCA/MCH/DCA_y_vs_x"))->Fill(dcax, dcay); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrackPar.getP(), quadrant, posNeg, dcax); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrackPar.getP(), quadrant, posNeg, dcay); + + if (cfgEnableMftMchResidualsExtraPlots) { + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcax); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcay); + const auto mchTrackAtMFT = PropagateMCHParam(FwdtoMCH(mchTrackPar), mftTrack.z()); + const double deltaPhi = mchTrackAtMFT.getPhi() - mftTrack.phi(); + registry.get(HIST("residuals/dphi_at_mft"))->Fill(deltaPhi, mftTrack.x(), mftTrack.y(), posNeg, mchTrackAtMFT.getP()); } - if (convertedTrackWithCorrOk) { - auto mchTrackAtDCA = PropagateMCHRealigned(convertedTrackWithCorr, collision.posZ()); - auto dcax = mchTrackAtDCA.getX() - collision.posX(); - auto dcay = mchTrackAtDCA.getY() - collision.posY(); + if (!mchTrackParNew.isRemovable()) { + const auto mchTrackAtDCAWithCorr = PropagateMCHParam(FwdtoMCH(mchTrackParNew), collision.posZ()); + const auto dcax = mchTrackAtDCAWithCorr.getX() - collision.posX(); + const auto dcay = mchTrackAtDCAWithCorr.getY() - collision.posY(); - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrant, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrant, posNeg, dcay); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcax); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcay); } } // MFT-MCH track residuals analysis - if (fEnableMftMchMatchingAnalysis && convertedTrackWithCorrOk) { - double refPlaneZ[2] = {fRefPlaneZMFT, fRefPlaneZMCH}; - - std::shared_ptr dxPlots[2]{registry.get(HIST("matching/dxAtMFT")), registry.get(HIST("matching/dxAtMCH"))}; - std::shared_ptr dyPlots[2]{registry.get(HIST("matching/dyAtMFT")), registry.get(HIST("matching/dyAtMCH"))}; - std::shared_ptr dsxPlots[2]{registry.get(HIST("matching/dsxAtMFT")), registry.get(HIST("matching/dsxAtMCH"))}; - std::shared_ptr dsyPlots[2]{registry.get(HIST("matching/dsyAtMFT")), registry.get(HIST("matching/dsyAtMCH"))}; - std::shared_ptr dphiPlots[2]{registry.get(HIST("matching/dphiAtMFT")), registry.get(HIST("matching/dphiAtMCH"))}; - - for (int iRefPlane = 0; iRefPlane < 2; iRefPlane++) { - const auto mftTrackAtRefPlane = configRealign.fEnableMCHRealign ? PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrackWithCorr.first()), refPlaneZ[iRefPlane]) : PropagateMFTtoMCH(mftTrack, FwdtoMCH(FwdToTrackPar(mchTrack)), refPlaneZ[iRefPlane]); - const auto mchTrackAtRefPlane = configRealign.fEnableMCHRealign ? PropagateMCHRealigned(convertedTrackWithCorr, refPlaneZ[iRefPlane]) : PropagateMCH(mchTrack, refPlaneZ[iRefPlane]); + if (cfgEnableMftMchMatchingAnalysis && !mchTrackParNew.isRemovable() && (matchChi2 <= 10.f)) { + static constexpr int nRefPlanes = 2; + const std::array refPlaneZ{cfgRefPlaneZMFT, cfgRefPlaneZMCH}; + + std::array, 2> dxPlots{registry.get(HIST("matching/dxAtMFT")), registry.get(HIST("matching/dxAtMCH"))}; + std::array, 2> dyPlots{registry.get(HIST("matching/dyAtMFT")), registry.get(HIST("matching/dyAtMCH"))}; + std::array, 2> dsxPlots{registry.get(HIST("matching/dsxAtMFT")), registry.get(HIST("matching/dsxAtMCH"))}; + std::array, 2> dsyPlots{registry.get(HIST("matching/dsyAtMFT")), registry.get(HIST("matching/dsyAtMCH"))}; + std::array, 2> dphiPlots{registry.get(HIST("matching/dphiAtMFT")), registry.get(HIST("matching/dphiAtMCH"))}; + + for (int iRefPlane = 0; iRefPlane < nRefPlanes; iRefPlane++) { + const auto mftTrackAtRefPlane = PropagateMFTtoMCH(mftTrackParNew, FwdtoMCH(mchTrackParNew), refPlaneZ[iRefPlane]); + const auto mchTrackAtRefPlane = PropagateMCHParam(FwdtoMCH(mchTrackParNew), refPlaneZ[iRefPlane]); const auto& refTrackAtRefPlane = (iRefPlane == 0) ? mftTrackAtRefPlane : mchTrackAtRefPlane; auto dx = mchTrackAtRefPlane.getX() - mftTrackAtRefPlane.getX(); - dxPlots[iRefPlane]->Fill(dx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dxPlots[iRefPlane]->Fill(dx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); auto dy = mchTrackAtRefPlane.getY() - mftTrackAtRefPlane.getY(); - dyPlots[iRefPlane]->Fill(dy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dyPlots[iRefPlane]->Fill(dy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); - auto mftParamAtRefPlane = FwdtoMCH(mftTrackAtRefPlane); - auto mchParamAtRefPlane = FwdtoMCH(mchTrackAtRefPlane); + const auto mftParamAtRefPlane = FwdtoMCH(mftTrackAtRefPlane); + const auto mchParamAtRefPlane = FwdtoMCH(mchTrackAtRefPlane); auto dsx = mchParamAtRefPlane.getNonBendingSlope() - mftParamAtRefPlane.getNonBendingSlope(); - dsxPlots[iRefPlane]->Fill(dsx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dsxPlots[iRefPlane]->Fill(dsx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); auto dsy = mchParamAtRefPlane.getBendingSlope() - mftParamAtRefPlane.getBendingSlope(); - dsyPlots[iRefPlane]->Fill(dsy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dsyPlots[iRefPlane]->Fill(dsy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); - auto dphi = mchTrackAtRefPlane.getPhi() - mftTrackAtRefPlane.getPhi(); - if (dphi < -TMath::Pi()) { - dphi += TMath::Pi() * 2.0; - } else if (dphi > TMath::Pi()) { - dphi -= TMath::Pi() * 2.0; - } - dphiPlots[iRefPlane]->Fill(dphi, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + auto dphi = RecoDecay::constrainAngle(mchTrackAtRefPlane.getPhi() - mftTrackAtRefPlane.getPhi(), -o2::constants::math::PI); + dphiPlots[iRefPlane]->Fill(dphi, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); } } } } } +#define FILL_DIMUON_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, histName) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double mass = mumu4mom.M(); \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histName))->Fill(mass, p, pT, quadrant1, quadrant2); \ + } + +#define FILL_DIMUON_DCA_PLOTS(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, trackPar1AtDca, trackPar2AtDca, histNameX, histNameY) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double dcax = trackPar1AtDca.getX() - trackPar2AtDca.getX(); \ + double dcay = trackPar1AtDca.getY() - trackPar2AtDca.getY(); \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histNameX))->Fill(dcax, p, pT, quadrant1, quadrant2); \ + registry.get(HIST(histNameY))->Fill(dcay, p, pT, quadrant1, quadrant2); \ + } + +#define FILL_DIMUON_ANGLE_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, mchAngle, fwdAngle, histName) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double dAngle = mchAngle - fwdAngle; \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histName))->Fill(dAngle, fwdAngle, p, pT, quadrant1, quadrant2); \ + } + + void FillDimuonPlots(MyEvents const& collisions, + MyMuonsWithCov const& muonTracks, + const std::map& collisionInfos) + { + if (!cfgEnableDimuonAnalysis) { + return; + } + + for (const auto& [collisionIndex, collisionInfo] : collisionInfos) { + auto const& collision = collisions.rawIteratorAt(collisionIndex); + + std::vector muonPairs; + getMuonPairs(collisionInfo, muonPairs); + + for (const auto& [mchIndex1, mchIndex2] : muonPairs) { + + auto const& muonTrack1 = muonTracks.rawIteratorAt(mchIndex1); + auto const& muonTrack2 = muonTracks.rawIteratorAt(mchIndex2); + int sign1 = muonTrack1.sign(); + int sign2 = muonTrack2.sign(); + + // only consider opposite-sign pairs + if ((sign1 * sign2) >= 0) + continue; + + bool isGoodMuon1 = IsGoodMuon(muonTrack1, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool isGoodMuon2 = IsGoodMuon(muonTrack2, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool goodMuonTracks = (isGoodMuon1 && isGoodMuon2); + + if (!goodMuonTracks) { + continue; + } + + + // get the pre-stored MCH track parameters + const auto mchTrackParIt1 = mMchTrackPars.find(mchIndex1); + if (mchTrackParIt1 == mMchTrackPars.end()) { + continue; + } + const auto mchTrackParIt2 = mMchTrackPars.find(mchIndex2); + if (mchTrackParIt2 == mMchTrackPars.end()) { + continue; + } + const auto& mchTrackPar1 = mchTrackParIt1->second; + auto mchTrackPar1AtDca = PropagateMCHParam(FwdtoMCH(mchTrackPar1), collision.posZ()); + auto mchTrackPar1AtVertex = PropagateMCHToVertex(mchTrackPar1, collision); + const auto& mchTrackPar2 = mchTrackParIt2->second; + auto mchTrackPar2AtDca = PropagateMCHParam(FwdtoMCH(mchTrackPar2), collision.posZ()); + auto mchTrackPar2AtVertex = PropagateMCHToVertex(mchTrackPar2, collision); + + // get the pre-stored MCH track parameters after refit + const auto mchTrackParNewIt1 = mMchTrackParsNew.find(mchIndex1); + if (mchTrackParNewIt1 == mMchTrackParsNew.end()) { + continue; + } + const auto mchTrackParNewIt2 = mMchTrackParsNew.find(mchIndex2); + if (mchTrackParNewIt2 == mMchTrackParsNew.end()) { + continue; + } + const auto& mchTrackParNew1 = mchTrackParNewIt1->second; + auto mchTrackParNew1AtDca = PropagateMCHParam(FwdtoMCH(mchTrackParNew1), collision.posZ()); + auto mchTrackParNew1AtVertex = PropagateMCHToVertex(mchTrackParNew1, collision); + const auto& mchTrackParNew2 = mchTrackParNewIt2->second; + auto mchTrackParNew2AtDca = PropagateMCHParam(FwdtoMCH(mchTrackParNew2), collision.posZ()); + auto mchTrackParNew2AtVertex = PropagateMCHToVertex(mchTrackParNew2, collision); + + FILL_DIMUON_PLOT(mchTrackPar1, mchTrackPar2, mchTrackPar1AtVertex, mchTrackPar2AtVertex, "dimuon/invariantMass_MuonKine_MuonCuts"); + FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, "dimuon/realign/invariantMass_MuonKine_MuonCuts"); + + FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, + mchTrackPar1AtVertex, mchTrackPar2AtVertex, + mchTrackPar1AtDca, mchTrackPar2AtDca, + "dimuon/dcax_MuonKine_MuonCuts", "dimuon/dcay_MuonKine_MuonCuts"); + FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, + mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, + mchTrackParNew1AtDca, mchTrackParNew2AtDca, + "dimuon/realign/dcax_MuonKine_MuonCuts", "dimuon/realign/dcay_MuonKine_MuonCuts"); + + double mchAngle = getMuMuAngle(mchTrackPar1AtVertex, mchTrackPar2AtVertex); + double mchAngleNew = getMuMuAngle(mchTrackParNew1AtVertex, mchTrackParNew2AtVertex); + + try { + const auto& candidates1 = collisionInfo.globalMuonTracks.at(mchIndex1); + const auto& candidates2 = collisionInfo.globalMuonTracks.at(mchIndex2); + + auto fwdIndex1 = candidates1[0]; + auto fwdIndex2 = candidates2[0]; + + auto const& fwdTrack1 = muonTracks.rawIteratorAt(fwdIndex1); + auto const& fwdTrack2 = muonTracks.rawIteratorAt(fwdIndex2); + + bool isGoodGlobalMuon1 = IsGoodMuon(muonTrack1, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool isGoodGlobalMuon2 = IsGoodMuon(muonTrack2, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool goodGlobalMuonTracks = (isGoodGlobalMuon1 && isGoodGlobalMuon2); + + bool isGoodMatch1 = isGoodGlobalMatching(fwdTrack1, 50.f); + bool isGoodMatch2 = isGoodGlobalMatching(fwdTrack2, 50.f); + bool goodGlobalMuonMatches = (isGoodMatch1 && isGoodMatch2); + + if (!goodGlobalMuonTracks || !goodGlobalMuonMatches) { + continue; + } + + FILL_DIMUON_PLOT(mchTrackPar1, mchTrackPar2, mchTrackPar1AtVertex, mchTrackPar2AtVertex, "dimuon/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, "dimuon/realign/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches"); + + FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, + mchTrackPar1AtVertex, mchTrackPar2AtVertex, + mchTrackPar1AtDca, mchTrackPar2AtDca, + "dimuon/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, + mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, + mchTrackParNew1AtDca, mchTrackParNew2AtDca, + "dimuon/realign/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); + + auto mftIndex1 = fwdTrack1.matchMFTTrackId(); + auto mftIndex2 = fwdTrack2.matchMFTTrackId(); + + const auto mftTrackPar1 = mMftTrackPars.at(mftIndex1); + auto fwdTrackPar1AtDca = PropagateMFTToDCA(mftTrackPar1, mchTrackPar1, collision, cfgVertexZshift); + auto fwdTrackPar1AtVertex = PropagateMFTToVertex(mftTrackPar1, mchTrackPar1, collision); + const auto mftTrackPar2 = mMftTrackPars.at(mftIndex2); + auto fwdTrackPar2AtDca = PropagateMFTToDCA(mftTrackPar2, mchTrackPar2, collision, cfgVertexZshift); + auto fwdTrackPar2AtVertex = PropagateMFTToVertex(mftTrackPar2, mchTrackPar2, collision); + const auto mftTrackParNew1 = mMftTrackParsNew.at(mftIndex1); + auto fwdTrackParNew1AtDca = PropagateMFTToDCA(mftTrackParNew1, mchTrackParNew1, collision, cfgVertexZshift); + auto fwdTrackParNew1AtVertex = PropagateMFTToVertex(mftTrackParNew1, mchTrackParNew1, collision); + const auto mftTrackParNew2 = mMftTrackParsNew.at(mftIndex2); + auto fwdTrackParNew2AtDca = PropagateMFTToDCA(mftTrackParNew2, mchTrackParNew2, collision, cfgVertexZshift); + auto fwdTrackParNew2AtVertex = PropagateMFTToVertex(mftTrackParNew2, mchTrackParNew2, collision); + + FILL_DIMUON_PLOT(mchTrackPar1, mchTrackPar2, fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, "dimuon/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, "dimuon/realign/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + + FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, + fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, + fwdTrackPar1AtDca, fwdTrackPar2AtDca, + "dimuon/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, + fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, + fwdTrackParNew1AtDca, fwdTrackParNew2AtDca, + "dimuon/realign/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + + double fwdAngle = getMuMuAngle(fwdTrackPar1AtVertex, fwdTrackPar2AtVertex); + double fwdAngleNew = getMuMuAngle(fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex); + + FILL_DIMUON_ANGLE_PLOT(mchTrackPar1, mchTrackPar2, fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, mchAngle, fwdAngle, "dimuon/angle_GlobalMuonCuts_GoodMatches") + FILL_DIMUON_ANGLE_PLOT(mchTrackParNew1, mchTrackParNew2, fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, mchAngleNew, fwdAngleNew, "dimuon/realign/angle_GlobalMuonCuts_GoodMatches") + } catch (const std::exception&) { + continue; + } + } + } + } + void processQA(MyEvents const& collisions, MyBCs const& bcs, MyMuonsWithCov const& muonTracks, @@ -1883,14 +2348,16 @@ struct muonGlobalAlignment { } std::map collisionInfos; - InitCollisions(collisions, bcs, muonTracks, mftTracks, collisionInfos); + InitCollisions(collisions, bcs, muonTracks, clusters, mftTracks, collisionInfos); + + FillMftPlots(collisions, bcs, muonTracks, mftTracks, collisionInfos); - FillDCAPlots(collisions, bcs, muonTracks, mftTracks, collisionInfos); + FillMchPlots(collisions, bcs, muonTracks, clusters, collisionInfos); - FillResidualsPlots(collisions, bcs, muonTracks, clusters, collisionInfos); + FillDimuonPlots(collisions, muonTracks, collisionInfos); } - PROCESS_SWITCH(muonGlobalAlignment, processQA, "process qa", true); + PROCESS_SWITCH(muonGlobalAlignment, processQA, "processQA", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 577fbd868b4776127ea1bcd4794bf3f486c63c53 Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Tue, 11 Aug 2026 17:53:21 +0200 Subject: [PATCH 2/2] [PWGDQ] fixed linter and clang errors --- PWGDQ/Tasks/muonGlobalAlignment.cxx | 127 ++++++++++++++-------------- 1 file changed, 64 insertions(+), 63 deletions(-) diff --git a/PWGDQ/Tasks/muonGlobalAlignment.cxx b/PWGDQ/Tasks/muonGlobalAlignment.cxx index dd8cc3bee63..13b32a0d244 100644 --- a/PWGDQ/Tasks/muonGlobalAlignment.cxx +++ b/PWGDQ/Tasks/muonGlobalAlignment.cxx @@ -18,9 +18,9 @@ #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/RCTSelectionFlags.h" #include "Common/Core/RecoDecay.h" +#include "Common/Core/fwdtrackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/fwdtrackUtilities.h" #include #include @@ -81,6 +81,7 @@ #include #include #include +#include #include using namespace o2; @@ -111,7 +112,6 @@ using MyMFTCovariance = MyMFTCovariances::iterator; using SMatrix55 = ROOT::Math::SMatrix>; using SMatrix5 = ROOT::Math::SVector; - using o2::dataformats::GlobalFwdTrack; using o2::track::TrackParCovFwd; @@ -153,7 +153,9 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc Configurable cfgTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; Configurable cfgTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; - Configurable cfgMftDcaMatchChi2Up{"cfgMftDcaMatchChi2Up", 10.f, ""}; + Configurable cfgMftDcaMatchChi2Up{"cfgMftDcaMatchChi2Up", 50.f, ""}; + Configurable cfgMftMchResidualsMatchChi2Up{"cfgMftMchResidualsMatchChi2Up", 50.f, ""}; + Configurable cfgDimuonMatchChi2Up{"cfgDimuonMatchChi2Up", 50.f, ""}; Configurable cfgMftMchResidualsPLow{"cfgMftMchResidualsPLow", 30.f, ""}; Configurable cfgMftMchResidualsPtLow{"cfgMftMchResidualsPtLow", 4.f, ""}; @@ -1963,7 +1965,7 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc const auto& mftTrack = muonTrack.template matchMFTTrack_as(); // int quadrant = GetQuadrant(mchTrack); int quadrant = GetQuadrant(mftTrack); - //int quadrantMCH = GetQuadrant(mchTrack); + // int quadrantMCH = GetQuadrant(mchTrack); int posNeg = (mchTrack.sign() >= 0) ? 0 : 1; auto mchIndex = mchTrack.globalIndex(); @@ -2008,10 +2010,10 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc } const auto& mchTrackParNew = mchTrackParNewIt->second; - double matchChi2 = muonTrack.chi2MatchMCHMFT() / 5.f; + double matchChi2 = muonTrack.chi2MatchMCHMFT(); // Residuals analysis between MFT tracks and MCH clusters - if (cfgEnableMftMchResidualsAnalysis && (matchChi2 <= 10.f)) { + if (cfgEnableMftMchResidualsAnalysis && (matchChi2 <= cfgMftMchResidualsMatchChi2Up.value)) { // loop over attached clusters auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id for (auto const& cluster : clustersSliced) { @@ -2068,14 +2070,14 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc const auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrackParNew, FwdtoMCH(mchTrackParNew), masterRealign.z()); const auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr); - const std::array xPos{masterRealign.x(), mftTrackAtClusterWithCorr.getX()}; - const std::array yPos{masterRealign.y(), mftTrackAtClusterWithCorr.getY()}; + const std::array xPosWithCorr{masterRealign.x(), mftTrackAtClusterWithCorr.getX()}; + const std::array yPosWithCorr{masterRealign.y(), mftTrackAtClusterWithCorr.getY()}; - registry.get(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]); - registry.get(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]); + registry.get(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, xPosWithCorr[0] - xPosWithCorr[1]); + registry.get(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, yPosWithCorr[0] - yPosWithCorr[1]); - registry.get(HIST("residuals/dx_vs_de_corr"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getNonBendingSlope()); - registry.get(HIST("residuals/dy_vs_de_corr"))->Fill(yPos[0] - yPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getBendingSlope()); + registry.get(HIST("residuals/dx_vs_de_corr"))->Fill(xPosWithCorr[0] - xPosWithCorr[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getNonBendingSlope()); + registry.get(HIST("residuals/dy_vs_de_corr"))->Fill(yPosWithCorr[0] - yPosWithCorr[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getBendingSlope()); } } @@ -2097,16 +2099,16 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc if (!mchTrackParNew.isRemovable()) { const auto mchTrackAtDCAWithCorr = PropagateMCHParam(FwdtoMCH(mchTrackParNew), collision.posZ()); - const auto dcax = mchTrackAtDCAWithCorr.getX() - collision.posX(); - const auto dcay = mchTrackAtDCAWithCorr.getY() - collision.posY(); + const auto dcaxWithCorr = mchTrackAtDCAWithCorr.getX() - collision.posX(); + const auto dcayWithCorr = mchTrackAtDCAWithCorr.getY() - collision.posY(); - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcay); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcaxWithCorr); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcayWithCorr); } } // MFT-MCH track residuals analysis - if (cfgEnableMftMchMatchingAnalysis && !mchTrackParNew.isRemovable() && (matchChi2 <= 10.f)) { + if (cfgEnableMftMchMatchingAnalysis && !mchTrackParNew.isRemovable() && (matchChi2 <= cfgMftMchResidualsMatchChi2Up.value)) { static constexpr int nRefPlanes = 2; const std::array refPlaneZ{cfgRefPlaneZMFT, cfgRefPlaneZMCH}; @@ -2142,39 +2144,39 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc } } -#define FILL_DIMUON_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, histName) \ - { \ - auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ - double p = mumu4mom.P(); \ - double pT = mumu4mom.Pt(); \ - double mass = mumu4mom.M(); \ +#define FILL_DIMUON_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, histName) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double mass = mumu4mom.M(); \ int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ - registry.get(HIST(histName))->Fill(mass, p, pT, quadrant1, quadrant2); \ + registry.get(HIST(histName))->Fill(mass, p, pT, quadrant1, quadrant2); \ } #define FILL_DIMUON_DCA_PLOTS(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, trackPar1AtDca, trackPar2AtDca, histNameX, histNameY) \ - { \ - auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ - double p = mumu4mom.P(); \ - double pT = mumu4mom.Pt(); \ - double dcax = trackPar1AtDca.getX() - trackPar2AtDca.getX(); \ - double dcay = trackPar1AtDca.getY() - trackPar2AtDca.getY(); \ - int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ - int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ - registry.get(HIST(histNameX))->Fill(dcax, p, pT, quadrant1, quadrant2); \ - registry.get(HIST(histNameY))->Fill(dcay, p, pT, quadrant1, quadrant2); \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double dcax = trackPar1AtDca.getX() - trackPar2AtDca.getX(); \ + double dcay = trackPar1AtDca.getY() - trackPar2AtDca.getY(); \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histNameX))->Fill(dcax, p, pT, quadrant1, quadrant2); \ + registry.get(HIST(histNameY))->Fill(dcay, p, pT, quadrant1, quadrant2); \ } #define FILL_DIMUON_ANGLE_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, mchAngle, fwdAngle, histName) \ - { \ - auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ - double p = mumu4mom.P(); \ - double pT = mumu4mom.Pt(); \ - double dAngle = mchAngle - fwdAngle; \ - int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ - int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ - registry.get(HIST(histName))->Fill(dAngle, fwdAngle, p, pT, quadrant1, quadrant2); \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double dAngle = mchAngle - fwdAngle; \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histName))->Fill(dAngle, fwdAngle, p, pT, quadrant1, quadrant2); \ } void FillDimuonPlots(MyEvents const& collisions, @@ -2210,7 +2212,6 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc continue; } - // get the pre-stored MCH track parameters const auto mchTrackParIt1 = mMchTrackPars.find(mchIndex1); if (mchTrackParIt1 == mMchTrackPars.end()) { @@ -2247,13 +2248,13 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, "dimuon/realign/invariantMass_MuonKine_MuonCuts"); FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, - mchTrackPar1AtVertex, mchTrackPar2AtVertex, - mchTrackPar1AtDca, mchTrackPar2AtDca, - "dimuon/dcax_MuonKine_MuonCuts", "dimuon/dcay_MuonKine_MuonCuts"); + mchTrackPar1AtVertex, mchTrackPar2AtVertex, + mchTrackPar1AtDca, mchTrackPar2AtDca, + "dimuon/dcax_MuonKine_MuonCuts", "dimuon/dcay_MuonKine_MuonCuts"); FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, - mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, - mchTrackParNew1AtDca, mchTrackParNew2AtDca, - "dimuon/realign/dcax_MuonKine_MuonCuts", "dimuon/realign/dcay_MuonKine_MuonCuts"); + mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, + mchTrackParNew1AtDca, mchTrackParNew2AtDca, + "dimuon/realign/dcax_MuonKine_MuonCuts", "dimuon/realign/dcay_MuonKine_MuonCuts"); double mchAngle = getMuMuAngle(mchTrackPar1AtVertex, mchTrackPar2AtVertex); double mchAngleNew = getMuMuAngle(mchTrackParNew1AtVertex, mchTrackParNew2AtVertex); @@ -2272,8 +2273,8 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc bool isGoodGlobalMuon2 = IsGoodMuon(muonTrack2, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); bool goodGlobalMuonTracks = (isGoodGlobalMuon1 && isGoodGlobalMuon2); - bool isGoodMatch1 = isGoodGlobalMatching(fwdTrack1, 50.f); - bool isGoodMatch2 = isGoodGlobalMatching(fwdTrack2, 50.f); + bool isGoodMatch1 = isGoodGlobalMatching(fwdTrack1, cfgDimuonMatchChi2Up.value); + bool isGoodMatch2 = isGoodGlobalMatching(fwdTrack2, cfgDimuonMatchChi2Up.value); bool goodGlobalMuonMatches = (isGoodMatch1 && isGoodMatch2); if (!goodGlobalMuonTracks || !goodGlobalMuonMatches) { @@ -2284,13 +2285,13 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, "dimuon/realign/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches"); FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, - mchTrackPar1AtVertex, mchTrackPar2AtVertex, - mchTrackPar1AtDca, mchTrackPar2AtDca, - "dimuon/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); + mchTrackPar1AtVertex, mchTrackPar2AtVertex, + mchTrackPar1AtDca, mchTrackPar2AtDca, + "dimuon/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, - mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, - mchTrackParNew1AtDca, mchTrackParNew2AtDca, - "dimuon/realign/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); + mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, + mchTrackParNew1AtDca, mchTrackParNew2AtDca, + "dimuon/realign/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); auto mftIndex1 = fwdTrack1.matchMFTTrackId(); auto mftIndex2 = fwdTrack2.matchMFTTrackId(); @@ -2312,13 +2313,13 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, "dimuon/realign/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches"); FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, - fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, - fwdTrackPar1AtDca, fwdTrackPar2AtDca, - "dimuon/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, + fwdTrackPar1AtDca, fwdTrackPar2AtDca, + "dimuon/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, - fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, - fwdTrackParNew1AtDca, fwdTrackParNew2AtDca, - "dimuon/realign/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, + fwdTrackParNew1AtDca, fwdTrackParNew2AtDca, + "dimuon/realign/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); double fwdAngle = getMuMuAngle(fwdTrackPar1AtVertex, fwdTrackPar2AtVertex); double fwdAngleNew = getMuMuAngle(fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex);