diff --git a/Common/DataModel/MatchMFTMuonData.h b/Common/DataModel/MatchMFTMuonData.h new file mode 100644 index 00000000000..645f26559b0 --- /dev/null +++ b/Common/DataModel/MatchMFTMuonData.h @@ -0,0 +1,138 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +// \file MatchMFTFT0.h +// \author Sarah Herrmann +// +// \brief Declaration of tables useful for the matching of MFT tracks and MUON tracks +// \date 11/12/24 + +#include "Framework/AnalysisDataModel.h" +namespace o2::aod +{ +namespace matching_params +{ +// matching parameters +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); +} // namespace matching_params + +DECLARE_SOA_TABLE(MatchParams, "AOD", "MATCHING", + matching_params::GMuonPt, + matching_params::GMuonEta, + matching_params::PairQ, + matching_params::DeltaPt, + matching_params::DeltaX, + matching_params::DeltaY, + matching_params::DeltaEta, + matching_params::DeltaPhi, + matching_params::IsCorrectMatch); + +namespace tag_matching_params +{ +// matching parameters +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); +} // namespace tag_matching_params + +DECLARE_SOA_TABLE(TagMatchParams, "AOD", "TAGMATCHING", + tag_matching_params::GMuonPt, + tag_matching_params::GMuonEta, + tag_matching_params::PairQ, + tag_matching_params::DeltaPt, + tag_matching_params::DeltaX, + tag_matching_params::DeltaY, + tag_matching_params::DeltaEta, + tag_matching_params::DeltaPhi, + tag_matching_params::IsCorrectMatch); + +namespace probe_matching_params +{ +// matching parameters +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(TagGMuonPt, mTagGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); +} // namespace probe_matching_params + +DECLARE_SOA_TABLE(ProbeMatchParams, "AOD", "PROBEMATCHING", + probe_matching_params::TagGMuonPt, + probe_matching_params::GMuonPt, + probe_matching_params::GMuonEta, + probe_matching_params::PairQ, + probe_matching_params::DeltaPt, + probe_matching_params::DeltaX, + probe_matching_params::DeltaY, + probe_matching_params::DeltaEta, + probe_matching_params::DeltaPhi, + probe_matching_params::IsCorrectMatch); + +namespace mix_matching_params +{ +// matching parameters +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); +} // namespace mix_matching_params + +DECLARE_SOA_TABLE(MixMatchParams, "AOD", "MIXMATCHING", + mix_matching_params::GMuonPt, + mix_matching_params::GMuonEta, + mix_matching_params::PairQ, + mix_matching_params::DeltaPt, + mix_matching_params::DeltaX, + mix_matching_params::DeltaY, + mix_matching_params::DeltaEta, + mix_matching_params::DeltaPhi, + mix_matching_params::IsCorrectMatch); + +namespace muon_pair +{ +// matching parameters +DECLARE_SOA_COLUMN(Mass, mMass, float); +DECLARE_SOA_COLUMN(Pt, mPt, float); +DECLARE_SOA_COLUMN(Rap, mRap, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +} // namespace muon_pair + +DECLARE_SOA_TABLE(MuonPair, "AOD", "MUONPAIR", + muon_pair::PairQ, + muon_pair::Mass, + muon_pair::Pt, + muon_pair::Rap); + +} // namespace o2::aod diff --git a/Common/TableProducer/CMakeLists.txt b/Common/TableProducer/CMakeLists.txt index 91513029f73..d0a9e6f5440 100644 --- a/Common/TableProducer/CMakeLists.txt +++ b/Common/TableProducer/CMakeLists.txt @@ -134,3 +134,8 @@ o2physics_add_dpl_workflow(mftmch-matching-data SOURCES match-mft-mch-data.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(mftmch-matching-data-mc + SOURCES match-mft-mch-data-mc.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) diff --git a/Common/TableProducer/match-mft-mch-data-mc.cxx b/Common/TableProducer/match-mft-mch-data-mc.cxx new file mode 100644 index 00000000000..3cfd0c3409f --- /dev/null +++ b/Common/TableProducer/match-mft-mch-data-mc.cxx @@ -0,0 +1,841 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +#include +#include +#include +#include +#include + +#include "CCDB/BasicCCDBManager.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Centrality.h" +#include "Common/CCDB/TriggerAliases.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "Common/DataModel/MftmchMatchingML.h" +#include "Common/DataModel/MatchMFTMuonData.h" +#include "Common/Core/trackUtilities.h" +#include "PWGDQ/DataModel/ReducedInfoTables.h" +#include "PWGDQ/Core/VarManager.h" +#include "PWGDQ/Core/HistogramManager.h" +#include "PWGDQ/Core/AnalysisCut.h" +#include "PWGDQ/Core/AnalysisCompositeCut.h" +#include "PWGDQ/Core/HistogramsLibrary.h" +#include "PWGDQ/Core/CutsLibrary.h" +#include "DataFormatsGlobalTracking/RecoContainerCreateTracksVariadic.h" +#include "DetectorsVertexing/VertexTrackMatcher.h" +#include "ReconstructionDataFormats/PrimaryVertex.h" +#include "ReconstructionDataFormats/VtxTrackIndex.h" +#include "ReconstructionDataFormats/VtxTrackRef.h" +#include "DataFormatsITSMFT/ROFRecord.h" +#include "CommonDataFormat/InteractionRecord.h" +#include "DetectorsVertexing/PVertexerParams.h" +#include "MathUtils/Primitive2D.h" +#include "DataFormatsGlobalTracking/RecoContainer.h" +#include "Common/DataModel/CollisionAssociationTables.h" +#include "Common/DataModel/MatchMFTFT0.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DataFormatsParameters/GRPObject.h" +#include "Field/MagneticField.h" +#include "TGeoGlobalMagField.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsBase/GeometryManager.h" +#include "EventFiltering/Zorro.h" +#include "ReconstructionDataFormats/TrackFwd.h" +#include "Math/MatrixFunctions.h" +#include "Math/SMatrix.h" +#include "MFTTracking/Tracker.h" +#include "MCHTracking/TrackParam.h" +#include "MCHTracking/TrackExtrap.h" +#include "GlobalTracking/MatchGlobalFwd.h" +#include +#include +#include "TDatabasePDG.h" + +using namespace std; + +using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; + +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/DataTypes.h" +#include "Framework/runDataProcessing.h" + +using MyCollisions = aod::Collisions; +using MyBCs = soa::Join; +using MyMUONs = soa::Join; +using MyMFTs = soa::Join; + +using MyCollision = MyCollisions::iterator; +using MyBC = MyBCs::iterator; +using MyMUON = MyMUONs::iterator; +using MyMFT = MyMFTs::iterator; + +using SMatrix55 = ROOT::Math::SMatrix>; +using SMatrix5 = ROOT::Math::SVector; + +float mMu = TDatabasePDG::Instance()->GetParticle(13)->Mass(); +int mRunNumber; +/* + TLorentzVector muon1LV; + TLorentzVector muon2LV; + TLorentzVector dimuonLV; +*/ +unordered_map> map_mfttracks; +unordered_map> map_muontracks; +unordered_map map_collisions; +unordered_map map_has_mfttracks_collisions; +unordered_map map_has_muontracks_collisions; +unordered_map map_vtxz; +unordered_map map_nmfttrack; + +struct match_mft_mch_data_mc { + + //// Variables for matching method + Configurable fMatchingMethod{"cfgMatchingMethod", 0, ""}; + + //// Variables for selecting muon tracks + Configurable fEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; + Configurable fEtaMchUp{"cfgEtaMchUp", -2.5f, ""}; + Configurable fRabsLow1{"cfgRabsLow1", 17.6f, ""}; + Configurable fRabsUp1{"cfgRabsUp1", 26.5f, ""}; + Configurable fRabsLow2{"cfgRabsLow2", 26.5f, ""}; + Configurable fRabsUp2{"cfgRabsUp2", 89.5f, ""}; + Configurable fPdcaUp1{"cfgPdcaUp1", 594.f, ""}; + Configurable fPdcaUp2{"cfgPdcaUp2", 324.f, ""}; + Configurable fTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; + Configurable fMatchingChi2MchMidUp{"cfgMatchingChi2MchMidUp", 999.f, ""}; + + //// Variables for selecting mft tracks + Configurable fEtaMftLow{"cfgEtaMftlow", -3.6f, ""}; + Configurable fEtaMftUp{"cfgEtaMftup", -2.5f, ""}; + Configurable fTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; + Configurable fTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; + + /// Variables to add preselection for the matching table + Configurable fPreselectMatchingX{"cfgPreselectMatchingX", 15.f, ""}; + Configurable fPreselectMatchingY{"cfgPreselectMatchingY", 15.f, ""}; + + /// Variables to event mixing criteria + Configurable fSaveMixedMatchingParamsRate{"cfgSaveMixedMatchingParamsRate", 0.002f, ""}; + Configurable fEventMaxDeltaNMFT{"cfgEventMaxDeltaNMFT", 1, ""}; + Configurable fEventMaxDeltaVtxZ{"cfgEventMaxDeltaVtxZ", 1.f, ""}; + + //// Variables for selecting tag muon + Configurable fTagMassWindowMin{"cfgTagMassWindowMin", 2.8f, ""}; + Configurable fTagMassWindowMax{"cfgTagMassWindowMax", 3.3f, ""}; + + //// Variables for ccdb + Configurable ccdburl{"ccdb-url", "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"}; + + //// Variables for Tag matching criteria + Configurable fSigmaXTagMuonCut{"cfgSigmaXTagMuonCut", 1.f, ""}; + Configurable fMeanXTagMuonCut{"cfgMeanXTagMuonCut", 0.f, ""}; + Configurable fSigmaYTagMuonCut{"cfgSigmaYTagMuonCut", 1.f, ""}; + Configurable fMeanYTagMuonCut{"cfgMeanYTagMuonCut", 1.f, ""}; + + Configurable fSigmaEtaTagMuonCut{"cfgSigmaEtaTagMuonCut", 0.2f, ""}; + Configurable fMeanEtaTagMuonCut{"cfgMeanEtaTagMuonCut", 0.f, ""}; + Configurable fSigmaPhiTagMuonCut{"cfgSigmaPhiTagMuonCut", 0.2f, ""}; + Configurable fMeanPhiTagMuonCut{"cfgMeanPhiTagMuonCut", 0.f, ""}; + + template + class FindTagAndProbe + { + private: + o2::dataformats::GlobalFwdTrack muontrack_at_pv[2]; + + TLorentzVector mDimuon; + MUON muontrack1; + MUON muontrack2; + Collision collision; + int tagIdx, probeIdx; + + int16_t mQ; + + inline void fillCovarianceArray(MUON const& muontrack, float cov[15]) const + { + cov[0] = muontrack.cXX(); + cov[1] = muontrack.cXY(); + cov[2] = muontrack.cYY(); + cov[3] = muontrack.cPhiX(); + cov[4] = muontrack.cPhiY(); + cov[5] = muontrack.cPhiPhi(); + cov[6] = muontrack.cTglX(); + cov[7] = muontrack.cTglY(); + cov[8] = muontrack.cTglPhi(); + cov[9] = muontrack.cTglTgl(); + cov[10] = muontrack.c1PtX(); + cov[11] = muontrack.c1PtY(); + cov[12] = muontrack.c1PtPhi(); + cov[13] = muontrack.c1PtTgl(); + cov[14] = muontrack.c1Pt21Pt2(); + } + + inline o2::dataformats::GlobalFwdTrack propagateMUONtoPV(MUON const& muontrack) const + { + const double mz = muontrack.z(); + const double mchi2 = muontrack.chi2(); + const float mx = muontrack.x(); + const float my = muontrack.y(); + const float mphi = muontrack.phi(); + const float mtgl = muontrack.tgl(); + const float m1pt = muontrack.signed1Pt(); + + float cov[15]; + fillCovarianceArray(muontrack, cov); + SMatrix5 tpars(mx, my, mphi, mtgl, m1pt); + SMatrix55 tcovs(cov, cov + 15); + + o2::track::TrackParCovFwd parcovmuontrack{mz, tpars, tcovs, mchi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + + o2::globaltracking::MatchGlobalFwd mMatching; + auto mchtrack = mMatching.FwdtoMCH(gtrack); + + o2::mch::TrackExtrap::extrapToVertex(mchtrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + o2::dataformats::GlobalFwdTrack extrap_muontrack; + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + + return extrap_muontrack; + } + + inline void setTagAndProbe() + { + if (muontrack1.pt() > muontrack2.pt()) { + tagIdx = 0; + probeIdx = 1; + } else { + tagIdx = 1; + probeIdx = 0; + } + } + + public: + inline FindTagAndProbe(const MUON& muon1, const MUON& muon2, const Collision& coll) + : muontrack_at_pv(), mDimuon(), muontrack1(muon1), muontrack2(muon2), collision(coll), tagIdx(-1), probeIdx(-1), mQ(0) + { + mQ = muontrack1.sign() + muontrack2.sign(); + setTagAndProbe(); + } + + void calcMuonPairAtPV() + { + muontrack_at_pv[0] = propagateMUONtoPV(muontrack1); + muontrack_at_pv[1] = propagateMUONtoPV(muontrack2); + TLorentzVector vMuon1, vMuon2; + vMuon1.SetPtEtaPhiM(muontrack_at_pv[0].getPt(), muontrack_at_pv[0].getEta(), muontrack_at_pv[0].getPhi(), mMu); + vMuon2.SetPtEtaPhiM(muontrack_at_pv[1].getPt(), muontrack_at_pv[1].getEta(), muontrack_at_pv[1].getPhi(), mMu); + mDimuon = vMuon1 + vMuon2; + } + inline int getTagMuonIndex() const { return tagIdx; } + inline int getProbeMuonIndex() const { return probeIdx; } + inline float getMass() const { return mDimuon.M(); } + inline float getPt() const { return mDimuon.Pt(); } + inline float getRap() const { return mDimuon.Rapidity(); } + inline int16_t getCharge() const { return mQ; } + inline const o2::dataformats::GlobalFwdTrack& getMuonAtPV(int idx) const { return muontrack_at_pv[idx]; } + }; // end of class FindTagAndProbe + + template + class MatchingParamsML + { + private: + MUON muontrack; + MFT mfttrack; + Collision collision; + + float mDX, mDY, mDPt, mDPhi, mDEta; + float mGlobalMuonPtAtDCA, mGlobalMuonEtaAtDCA, mGlobalMuonPhiAtDCA, mGlobalMuonDCAx, mGlobalMuonDCAy, mGlobalMuonQ; + int mMatchingType; + + o2::field::MagneticField* fieldB; + o2::globaltracking::MatchGlobalFwd mMatching; + + inline o2::track::TrackParCovFwd propagateMFTtoMatchingPlane() + { + double covArr[15]{0.0}; + SMatrix55 tmftcovs(covArr, covArr + 15); + + SMatrix5 tmftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); + o2::track::TrackParCovFwd extrap_mfttrack{mfttrack.z(), tmftpars, tmftcovs, mfttrack.chi2()}; + + double propVec[3] = {0.}; + float zPlane = 0.f; + if (mMatchingType == MCH_FIRST_CLUSTER) { + propVec[0] = muontrack.x() - mfttrack.x(); + propVec[1] = muontrack.y() - mfttrack.y(); + propVec[2] = muontrack.z() - mfttrack.z(); + zPlane = muontrack.z(); + } else if (mMatchingType == END_OF_ABSORBER || mMatchingType == BEGINING_OF_ABSORBER) { + auto extrap_muontrack = propagateMUONtoMatchingPlane(); + propVec[0] = extrap_muontrack.getX() - mfttrack.x(); + propVec[1] = extrap_muontrack.getY() - mfttrack.y(); + propVec[2] = extrap_muontrack.getZ() - mfttrack.z(); + zPlane = (mMatchingType == END_OF_ABSORBER) ? -505.f : -90.f; + } else { + zPlane = mfttrack.z(); + } + + double centerZ[3] = {mfttrack.x() + propVec[0] / 2., mfttrack.y() + propVec[1] / 2., mfttrack.z() + propVec[2] / 2.}; + float Bz = fieldB->getBz(centerZ); + extrap_mfttrack.propagateToZ(zPlane, Bz); // z in cm + return extrap_mfttrack; + } + + inline o2::dataformats::GlobalFwdTrack propagateMUONtoMatchingPlane() + { + float cov[15] = { + muontrack.cXX(), muontrack.cXY(), muontrack.cYY(), + muontrack.cPhiX(), muontrack.cPhiY(), muontrack.cPhiPhi(), + muontrack.cTglX(), muontrack.cTglY(), muontrack.cTglPhi(), + muontrack.cTglTgl(), muontrack.c1PtX(), muontrack.c1PtY(), + muontrack.c1PtPhi(), muontrack.c1PtTgl(), muontrack.c1Pt21Pt2()}; + + SMatrix5 tpars(muontrack.x(), muontrack.y(), muontrack.phi(), muontrack.tgl(), muontrack.signed1Pt()); + SMatrix55 tcovs(cov, cov + 15); + double chi2 = muontrack.chi2(); + + o2::track::TrackParCovFwd parcovmuontrack{muontrack.z(), tpars, tcovs, chi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + + auto mchtrack = mMatching.FwdtoMCH(gtrack); + + if (mMatchingType == MFT_LAST_CLUSTR) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchtrack, mfttrack.z()); + } else if (mMatchingType == END_OF_ABSORBER) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchtrack, -505.); + } else if (mMatchingType == BEGINING_OF_ABSORBER) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchtrack, -90.); + } + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + + o2::dataformats::GlobalFwdTrack extrap_muontrack; + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + return extrap_muontrack; + } + + inline o2::track::TrackParCovFwd propagateMFTtoDCA() + { + double covArr[15]{0.0}; + SMatrix55 tmftcovs(covArr, covArr + 15); + + SMatrix5 tmftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); + o2::track::TrackParCovFwd extrap_mfttrack{mfttrack.z(), tmftpars, tmftcovs, mfttrack.chi2()}; + + double propVec[3] = {}; + propVec[0] = collision.posX() - mfttrack.x(); + propVec[1] = collision.posY() - mfttrack.y(); + propVec[2] = collision.posZ() - mfttrack.z(); + + double centerZ[3] = {mfttrack.x() + propVec[0] / 2., mfttrack.y() + propVec[1] / 2., mfttrack.z() + propVec[2] / 2.}; + float Bz = fieldB->getBz(centerZ); + extrap_mfttrack.propagateToZ(collision.posZ(), Bz); // z in cm + return extrap_mfttrack; + } + + inline o2::dataformats::GlobalFwdTrack propagateMUONtoPV() + { + float cov[15] = { + muontrack.cXX(), muontrack.cXY(), muontrack.cYY(), + muontrack.cPhiX(), muontrack.cPhiY(), muontrack.cPhiPhi(), + muontrack.cTglX(), muontrack.cTglY(), muontrack.cTglPhi(), + muontrack.cTglTgl(), muontrack.c1PtX(), muontrack.c1PtY(), + muontrack.c1PtPhi(), muontrack.c1PtTgl(), muontrack.c1Pt21Pt2()}; + + SMatrix5 tpars(muontrack.x(), muontrack.y(), muontrack.phi(), muontrack.tgl(), muontrack.signed1Pt()); + SMatrix55 tcovs(cov, cov + 15); + double chi2 = muontrack.chi2(); + + o2::track::TrackParCovFwd parcovmuontrack{muontrack.z(), tpars, tcovs, chi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + + auto mchtrack = mMatching.FwdtoMCH(gtrack); + o2::mch::TrackExtrap::extrapToVertex(mchtrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + o2::dataformats::GlobalFwdTrack extrap_muontrack; + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + + return extrap_muontrack; + } + + public: + enum MATCHING_TYPE { MCH_FIRST_CLUSTER, + MFT_LAST_CLUSTR, + END_OF_ABSORBER, + BEGINING_OF_ABSORBER }; + + MatchingParamsML(MUON const& muon, MFT const& mft, Collision const& coll, int MType, o2::field::MagneticField* field) : muontrack(muon), mfttrack(mft), collision(coll), mDX(0.f), mDY(0.f), mDPt(0.f), mDPhi(0.f), mDEta(0.f), mGlobalMuonPtAtDCA(0.f), mGlobalMuonEtaAtDCA(0.f), mGlobalMuonPhiAtDCA(0.f), mGlobalMuonDCAx(0.f), mGlobalMuonDCAy(0.f), mGlobalMuonQ(0.f), mMatchingType(MType), fieldB(field) {} + void calcMatchingParams() + { + auto mfttrack_on_matchingP = propagateMFTtoMatchingPlane(); + auto muontrack_on_matchingP = propagateMUONtoMatchingPlane(); + + float dphiRaw = mfttrack_on_matchingP.getPhi() - muontrack_on_matchingP.getPhi(); + float dphi = TVector2::Phi_mpi_pi(dphiRaw); + float deta = mfttrack_on_matchingP.getEta() - muontrack_on_matchingP.getEta(); + + mDX = mfttrack_on_matchingP.getX() - muontrack_on_matchingP.getX(); + mDY = mfttrack_on_matchingP.getY() - muontrack_on_matchingP.getY(); + mDPt = mfttrack_on_matchingP.getPt() - muontrack_on_matchingP.getPt(); + mDPhi = dphi; + mDEta = deta; + } + + void calcGlobalMuonParams() + { + auto mfttrack_at_dca = propagateMFTtoDCA(); + auto muontrack_at_pv = propagateMUONtoPV(); + + float momentum = muontrack_at_pv.getP(); + float theta = mfttrack_at_dca.getTheta(); + float phiTrack = mfttrack_at_dca.getPhi(); + float px = momentum * std::sin(theta) * std::cos(phiTrack); + float py = momentum * std::sin(theta) * std::sin(phiTrack); + + mGlobalMuonQ = muontrack.sign() + mfttrack.sign(); + mGlobalMuonPtAtDCA = std::sqrt(px * px + py * py); + mGlobalMuonEtaAtDCA = mfttrack_at_dca.getEta(); + mGlobalMuonPhiAtDCA = mfttrack_at_dca.getPhi(); + mGlobalMuonDCAx = mfttrack_at_dca.getX() - collision.posX(); + mGlobalMuonDCAy = mfttrack_at_dca.getY() - collision.posY(); + } + + inline float getDx() const { return mDX; } + inline float getDy() const { return mDY; } + inline float getDphi() const { return mDPhi; } + inline float getDeta() const { return mDEta; } + inline float getDpt() const { return mDPt; } + inline float getGMPtAtDCA() const { return mGlobalMuonPtAtDCA; } + inline float getGMEtaAtDCA() const { return mGlobalMuonEtaAtDCA; } + inline float getGMPhiAtDCA() const { return mGlobalMuonPhiAtDCA; } + inline float getGMDcaX() const { return mGlobalMuonDCAx; } + inline float getGMDcaY() const { return mGlobalMuonDCAy; } + inline float getGMDcaXY() const { return std::sqrt(mGlobalMuonDCAx * mGlobalMuonDCAx + mGlobalMuonDCAy * mGlobalMuonDCAy); } + inline int16_t getGMQ() const { return static_cast(mGlobalMuonQ); } + + }; // end of class MatchingParamsML + + template + o2::dataformats::GlobalFwdTrack propagateMUONtoPV(MUON const& muontrack, Collisions const& collisions) + { + auto collision = collisions.rawIteratorAt(muontrack.collisionId()); + o2::globaltracking::MatchGlobalFwd mMatching; + o2::dataformats::GlobalFwdTrack extrap_muontrack; + + SMatrix5 tpars(muontrack.x(), muontrack.y(), muontrack.phi(), muontrack.tgl(), muontrack.signed1Pt()); + std::vector v1{muontrack.cXX(), muontrack.cXY(), muontrack.cYY(), + muontrack.cPhiX(), muontrack.cPhiY(), muontrack.cPhiPhi(), + muontrack.cTglX(), muontrack.cTglY(), muontrack.cTglPhi(), + muontrack.cTglTgl(), muontrack.c1PtX(), muontrack.c1PtY(), + muontrack.c1PtPhi(), muontrack.c1PtTgl(), muontrack.c1Pt21Pt2()}; + SMatrix55 tcovs(v1.begin(), v1.end()); + double chi2 = muontrack.chi2(); + o2::track::TrackParCovFwd parcovmuontrack{muontrack.z(), tpars, tcovs, chi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + auto mchtrack = mMatching.FwdtoMCH(gtrack); + + o2::mch::TrackExtrap::extrapToVertex(mchtrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + + return extrap_muontrack; + } + + inline bool isGoodTagDimuon(float M) + { + return !(M < fTagMassWindowMin || M > fTagMassWindowMax); + } + + inline bool isGoodTagMatching(float mDX, float mDY, float mDEta, float mDPhi) + { + float dxNorm = (mDX - fMeanXTagMuonCut) / (fSigmaXTagMuonCut * 3); + float dyNorm = (mDY - fMeanYTagMuonCut) / (fSigmaYTagMuonCut * 3); + float detaNorm = (mDEta - fMeanEtaTagMuonCut) / (fSigmaEtaTagMuonCut * 3); + float dphiNorm = (mDPhi - fMeanPhiTagMuonCut) / (fSigmaPhiTagMuonCut * 3); + + float rTagXY = dxNorm * dxNorm + dyNorm * dyNorm; + float rTagEtaPhi = detaNorm * detaNorm + dphiNorm * dphiNorm; + + return (rTagXY < 1.f && rTagEtaPhi > 0.f); + } + + template + bool isCorrectMatching(MUON const& muontrack, MFT const& mfttrack) + { + + int idmuon = muontrack.mcParticleId(); + int idmft = mfttrack.mcParticleId(); + + if (idmuon == -1 || idmft == -1) + return false; + if (idmuon != idmft) + return false; + else + return true; + }; + + template + bool isGoodMuonQuality(MUON muontrack) + { + if (!muontrack.has_collision()) + return false; + if (muontrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) + return false; + if (muontrack.chi2() > fTrackChi2MchUp) + return false; + if (fRabsLow1 > muontrack.rAtAbsorberEnd() || muontrack.rAtAbsorberEnd() > fRabsUp2) + return false; + if (muontrack.rAtAbsorberEnd() < fRabsUp1 && fPdcaUp1 < muontrack.pDca()) + return false; + if (muontrack.rAtAbsorberEnd() > fRabsLow2 && fPdcaUp2 < muontrack.pDca()) + return false; + return true; + } + + template + bool isGoodMuonKine(MUON muontrack) + { + if (fEtaMchLow > muontrack.getEta() || muontrack.getEta() > fEtaMchUp) + return false; + return true; + } + + template + bool isGoodMFTQuality(MFT mfttrack) + { + if (!mfttrack.has_collision()) + return false; + if (mfttrack.chi2() > fTrackChi2MftUp) + return false; + if (mfttrack.nClusters() < fTrackNClustMftLow) + return false; + return true; + } + + template + bool isGoodMFTKine(MFT mfttrack) + { + if (fEtaMftLow > mfttrack.getEta() || mfttrack.getEta() > fEtaMftUp) + return false; + return true; + } + + inline bool isPassMatchingPreselection(float Dx, float Dy) + { + return !(std::abs(Dx) > fPreselectMatchingX || std::abs(Dy) > fPreselectMatchingY); + } + + template + void setMUONs(MUONs const& muontracks, Collisions const& collisions) + { + for (auto muontrack : muontracks) { + if (!isGoodMuonQuality(muontrack)) + continue; + o2::dataformats::GlobalFwdTrack muontrack_at_pv = propagateMUONtoPV(muontrack, collisions); + if (!isGoodMuonKine(muontrack_at_pv)) + continue; + + auto collision = collisions.rawIteratorAt(muontrack.collisionId()); + + bool& has = map_has_muontracks_collisions[muontrack.collisionId()]; + has = true; + + vector& arr_muontracks = map_muontracks[collision.globalIndex()]; + arr_muontracks.push_back(muontrack.globalIndex()); + } + } + + template + o2::track::TrackParCovFwd PropagateMFTtoDCA(MFT const& mfttrack, Collisions const& collisions, o2::field::MagneticField* field) + { + auto collision = collisions.rawIteratorAt(mfttrack.collisionId()); + std::vector mftv1; + SMatrix55 mftcovs{mftv1.begin(), mftv1.end()}; + SMatrix5 mftpars = {mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()}; + o2::track::TrackParCovFwd mftpartrack = {mfttrack.z(), mftpars, mftcovs, mfttrack.chi2()}; + double propVec[3] = {fabs(mfttrack.x() - collision.posX()), + fabs(mfttrack.y() - collision.posY()), + fabs(mfttrack.z() - collision.posZ())}; + double centerZ[3] = {mfttrack.x() - propVec[0] / 2., + mfttrack.y() - propVec[1] / 2., + mfttrack.z() - propVec[2] / 2.}; + float Bz = field->getBz(centerZ); + mftpartrack.propagateToZ(collision.posZ(), Bz); + return mftpartrack; + } + + template + void setMFTs(MFTs const& mfttracks, Collisions const& collisions, o2::field::MagneticField* field) + { + for (auto mfttrack : mfttracks) { + if (!isGoodMFTQuality(mfttrack)) + continue; + + o2::track::TrackParCovFwd mfttrack_at_dca = PropagateMFTtoDCA(mfttrack, collisions, field); + if (!isGoodMFTKine(mfttrack_at_dca)) + continue; + + auto collision = collisions.rawIteratorAt(mfttrack.collisionId()); + + map_vtxz[mfttrack.collisionId()] = collision.posZ(); + map_nmfttrack[mfttrack.collisionId()] += 1; + + bool& has = map_has_mfttracks_collisions[mfttrack.collisionId()]; + has = true; + + vector& arr_mfttracks = map_mfttracks[collision.globalIndex()]; + arr_mfttracks.push_back(mfttrack.globalIndex()); + } + } + + Produces tableMatchingParams; + Produces tableTagMatchingParams; + Produces tableProbeMatchingParams; + Produces tableMixMatchingParams; + Produces tableMuonPair; + + Service ccdbManager; + + o2::field::MagneticField* fieldB; + o2::ccdb::CcdbApi ccdbApi; + + template + void initCCDB(BC const& bc) + { + if (mRunNumber == bc.runNumber()) + return; + + mRunNumber = bc.runNumber(); + std::map metadata; + auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, mRunNumber); + auto ts = soreor.first; + auto grpmag = ccdbApi.retrieveFromTFileAny(grpmagPath, metadata, ts); + o2::base::Propagator::initFieldFromGRP(grpmag); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdbManager->get(geoPath); + } + o2::mch::TrackExtrap::setField(); + fieldB = static_cast(TGeoGlobalMagField::Instance()->GetField()); + } + + void init(o2::framework::InitContext&) + { + ccdbManager->setURL(ccdburl); + ccdbManager->setCaching(true); + ccdbManager->setLocalObjectValidityChecking(); + ccdbManager->setFatalWhenNull(false); + ccdbApi.init(ccdburl); + mRunNumber = 0; + } + + void process(MyCollisions const& collisions, + MyBCs const& bcs, + MyMUONs const& muontracks, + MyMFTs const& mfttracks) + { + LOG(info) << "Process() "; + map_muontracks.clear(); + map_mfttracks.clear(); + map_collisions.clear(); + map_has_muontracks_collisions.clear(); + map_has_mfttracks_collisions.clear(); + + initCCDB(bcs.begin()); + setMUONs(muontracks, collisions); + setMFTs(mfttracks, collisions, fieldB); + + for (auto map_has_muontracks_collision : map_has_muontracks_collisions) { + auto idmuontrack_collisions = map_has_muontracks_collision.first; + for (auto map_has_mfttracks_collision : map_has_mfttracks_collisions) { + auto idmfttrack_collisions = map_has_mfttracks_collision.first; + if (idmuontrack_collisions != idmfttrack_collisions) + continue; + map_collisions[idmfttrack_collisions] = true; + } + } + + for (auto const& map_collision : map_collisions) { + auto const& collision = collisions.rawIteratorAt(map_collision.first); + + for (auto const& imuontrack1 : map_muontracks[map_collision.first]) { + auto const& muontrack1 = muontracks.rawIteratorAt(imuontrack1); + + for (auto const& imfttrack1 : map_mfttracks[map_collision.first]) { + auto const& mfttrack1 = mfttracks.rawIteratorAt(imfttrack1); + + MatchingParamsML matching(muontrack1, mfttrack1, collision, fMatchingMethod, fieldB); + matching.calcMatchingParams(); + + if (!isPassMatchingPreselection(matching.getDx(), matching.getDy())) + continue; + + matching.calcGlobalMuonParams(); + + bool isTrue = isCorrectMatching(muontrack1, mfttrack1); + + tableMatchingParams(matching.getGMPtAtDCA(), + matching.getGMEtaAtDCA(), + static_cast(matching.getGMQ()), + matching.getDpt(), + matching.getDx(), + matching.getDy(), + matching.getDeta(), + matching.getDphi(), + isTrue); + } + + for (auto const& map_mfttrack : map_mfttracks) { + if (map_mfttrack.first == map_collision.first) + continue; + if (fabs(map_vtxz[map_mfttrack.first] - map_vtxz[map_collision.first]) > fEventMaxDeltaVtxZ) + continue; + if (fabs(map_nmfttrack[map_mfttrack.first] - map_nmfttrack[map_collision.first]) > fEventMaxDeltaNMFT) + continue; + + for (auto const& imfttrack1 : map_mfttrack.second) { + auto const& mfttrack1 = mfttracks.rawIteratorAt(imfttrack1); + MatchingParamsML matching(muontrack1, mfttrack1, collision, fMatchingMethod, fieldB); + matching.calcMatchingParams(); + if (!isPassMatchingPreselection(matching.getDx(), matching.getDy())) + continue; + matching.calcGlobalMuonParams(); + + bool isTrue = isCorrectMatching(muontrack1, mfttrack1); + + tableMixMatchingParams(matching.getGMPtAtDCA(), + matching.getGMEtaAtDCA(), + static_cast(matching.getGMQ()), + matching.getDpt(), + matching.getDx(), + matching.getDy(), + matching.getDeta(), + matching.getDphi(), + isTrue); + } + } + + for (auto const& imuontrack2 : map_muontracks[map_collision.first]) { + + if (imuontrack1 == imuontrack2) + continue; + + auto const& muontrack2 = muontracks.rawIteratorAt(imuontrack2); + + FindTagAndProbe tagdimuon(muontrack1, muontrack2, collision); + tagdimuon.calcMuonPairAtPV(); + tableMuonPair(tagdimuon.getCharge(), tagdimuon.getMass(), tagdimuon.getPt(), tagdimuon.getRap()); + + if (!isGoodTagDimuon(tagdimuon.getMass())) + continue; + + auto tagmuontrack = muontrack1; + auto probemuontrack = muontrack2; + + if (tagdimuon.getTagMuonIndex() == 1) { + tagmuontrack = muontrack2; + probemuontrack = muontrack1; + } + + int nTagMFTCand = 0; + int nProbeMFTCand = 0; + + unordered_map> map_tagMatchingParams; + unordered_map> map_probeMatchingParams; + + for (auto const& imfttrack1 : map_mfttracks[map_collision.first]) { + auto const& mfttrack1 = mfttracks.rawIteratorAt(imfttrack1); + MatchingParamsML matchingTag(tagmuontrack, mfttrack1, collision, fMatchingMethod, fieldB); + matchingTag.calcMatchingParams(); + matchingTag.calcGlobalMuonParams(); + if (isGoodTagMatching(matchingTag.getDx(), matchingTag.getDy(), matchingTag.getDeta(), matchingTag.getDphi()) && + isPassMatchingPreselection(matchingTag.getDx(), matchingTag.getDy())) { + bool isTrue = isCorrectMatching(tagmuontrack, mfttrack1); + tableTagMatchingParams(matchingTag.getGMPtAtDCA(), + matchingTag.getGMEtaAtDCA(), + matchingTag.getGMQ(), + matchingTag.getDpt(), + matchingTag.getDx(), + matchingTag.getDy(), + matchingTag.getDeta(), + matchingTag.getDphi(), + isTrue); + ++nTagMFTCand; + } + + if (nTagMFTCand == 1) { + MatchingParamsML matchingProbe(probemuontrack, mfttrack1, collision, fMatchingMethod, fieldB); + matchingProbe.calcMatchingParams(); + matchingProbe.calcGlobalMuonParams(); + if (isPassMatchingPreselection(matchingProbe.getDx(), matchingProbe.getDy())) { + bool isTrue = isCorrectMatching(probemuontrack, mfttrack1); + tableProbeMatchingParams(matchingTag.getGMPtAtDCA(), + matchingProbe.getGMPtAtDCA(), + matchingProbe.getGMEtaAtDCA(), + matchingProbe.getGMQ(), + matchingProbe.getDpt(), + matchingProbe.getDx(), + matchingProbe.getDy(), + matchingProbe.getDeta(), + matchingProbe.getDphi(), + isTrue); + } + } // nTagMFTCand + } // end of loop imfttrack1 + } // end of loop imuontrack2 + } // end of loop imuontrack1 + } // end of loop map_collision + + } // end of processMC +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/Common/TableProducer/match-mft-mch-data.cxx b/Common/TableProducer/match-mft-mch-data.cxx index d87504ccfef..e7427f81b83 100644 --- a/Common/TableProducer/match-mft-mch-data.cxx +++ b/Common/TableProducer/match-mft-mch-data.cxx @@ -14,11 +14,6 @@ #include #include -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/DataTypes.h" -#include "Framework/runDataProcessing.h" #include "CCDB/BasicCCDBManager.h" #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" @@ -27,6 +22,7 @@ #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" #include "Common/DataModel/MftmchMatchingML.h" +#include "Common/DataModel/MatchMFTMuonData.h" #include "Common/Core/trackUtilities.h" #include "PWGDQ/DataModel/ReducedInfoTables.h" #include "PWGDQ/Core/VarManager.h" @@ -65,277 +61,173 @@ #include #include "TDatabasePDG.h" +using namespace std; + using namespace o2; using namespace o2::soa; using namespace o2::aod; using namespace o2::framework; using namespace o2::framework::expressions; -using SMatrix55 = ROOT::Math::SMatrix>; -using SMatrix5 = ROOT::Math::SVector; +#include "Framework/AnalysisTask.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/DataTypes.h" +#include "Framework/runDataProcessing.h" -// using MyEvents = soa::Join; -// using MyEventsWithMults = soa::Join; -// using MyEventsWithFilter = soa::Join; -// using MyEventsWithMultsAndFilter = soa::Join; -// using MyEventsWithCent = soa::Join; -// using MyEventsWithCentAndMults = soa::Join; -using MyMuons = soa::Join; +using MyCollisions = aod::Collisions; +using MyBCs = soa::Join; +using MyMUONs = soa::Join; using MyMFTs = aod::MFTTracks; -// using MyMuonsWithCov = soa::Join; -// using MyMuonsColl = soa::Join; -// using MyMuonsCollWithCov = soa::Join; -using MyBCs = soa::Join; -using ExtBCs = soa::Join; -float mMu = TDatabasePDG::Instance()->GetParticle(13)->Mass(); - -TRandom* rnd = new TRandom(); +using MyCollision = MyCollisions::iterator; +using MyBC = MyBCs::iterator; +using MyMUON = MyMUONs::iterator; +using MyMFT = MyMFTs::iterator; -TLorentzVector muon1LV; -TLorentzVector muon2LV; -TLorentzVector dimuonLV; - -TVector3 V1; -TVector3 V2; +using SMatrix55 = ROOT::Math::SMatrix>; +using SMatrix5 = ROOT::Math::SVector; +float mMu = TDatabasePDG::Instance()->GetParticle(13)->Mass(); +int mRunNumber; +/* + TLorentzVector muon1LV; + TLorentzVector muon2LV; + TLorentzVector dimuonLV; +*/ +unordered_map> map_mfttracks; +unordered_map> map_muontracks; +unordered_map map_collisions; +unordered_map map_has_mfttracks_collisions; +unordered_map map_has_muontracks_collisions; +unordered_map map_vtxz; +unordered_map map_nmfttrack; +/* namespace o2::aod { - -namespace muon_params -{ -DECLARE_SOA_COLUMN(TRACKCHI2, trackChi2, float); -DECLARE_SOA_COLUMN(RABS, rabs, float); -DECLARE_SOA_COLUMN(Q, q, int16_t); -DECLARE_SOA_COLUMN(PT, pt, float); -DECLARE_SOA_COLUMN(ETA, eta, float); -DECLARE_SOA_COLUMN(PHI, phi, float); -DECLARE_SOA_COLUMN(HASMFT, has_mft, bool); -} // namespace muon_params - -DECLARE_SOA_TABLE(MUONParams, "AOD", "MUON", - muon_params::TRACKCHI2, - muon_params::RABS, - muon_params::Q, - muon_params::PT, - muon_params::ETA, - muon_params::PHI, - muon_params::HASMFT); - -namespace mft_params -{ -DECLARE_SOA_COLUMN(NCLUST, nclust, int); -DECLARE_SOA_COLUMN(ISCA, isCA, bool); -DECLARE_SOA_COLUMN(TRACKCHI2, trackChi2, float); -DECLARE_SOA_COLUMN(Q, q, int16_t); -DECLARE_SOA_COLUMN(PT_AT_DCA, pt_dca, float); -DECLARE_SOA_COLUMN(ETA_AT_DCA, eta_dca, float); -DECLARE_SOA_COLUMN(PHI_AT_DCA, phi_dca, float); -DECLARE_SOA_COLUMN(DCAx, dcax, float); -DECLARE_SOA_COLUMN(DCAy, dcay, float); -} // namespace mft_params - -DECLARE_SOA_TABLE(MFTParams, "AOD", "MFT", - mft_params::NCLUST, - mft_params::ISCA, - mft_params::TRACKCHI2, - mft_params::Q, - mft_params::PT_AT_DCA, - mft_params::ETA_AT_DCA, - mft_params::PHI_AT_DCA, - mft_params::DCAx, - mft_params::DCAy); - namespace matching_params { // matching parameters -DECLARE_SOA_COLUMN(NClustMFTTracks, nClustMFT, int); -DECLARE_SOA_COLUMN(Chi2MFTTracks, chi2MFT, float); - -DECLARE_SOA_COLUMN(DeltaP, dp_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPt, dpt_mchplane, float); -DECLARE_SOA_COLUMN(DeltaEta, deta_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPhi, dphi_mchplane, float); -DECLARE_SOA_COLUMN(DeltaX, dx_mchplane, float); -DECLARE_SOA_COLUMN(DeltaY, dy_mchplane, float); - -DECLARE_SOA_COLUMN(MchPt, mchpt, float); -DECLARE_SOA_COLUMN(MchEta, mcheta, float); -DECLARE_SOA_COLUMN(MchPhi, mchphi, float); -DECLARE_SOA_COLUMN(MchQ, mchq, float); - -DECLARE_SOA_COLUMN(MftPt, mftpt, float); -DECLARE_SOA_COLUMN(MftEta, mfteta, float); -DECLARE_SOA_COLUMN(MftPhi, mftphi, float); -DECLARE_SOA_COLUMN(MftQ, mftq, float); - -DECLARE_SOA_COLUMN(MftDCA, mftdca, float); - +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); } // namespace matching_params DECLARE_SOA_TABLE(MatchParams, "AOD", "MATCHING", + matching_params::GMuonPt, + matching_params::GMuonEta, + matching_params::PairQ, matching_params::DeltaPt, - matching_params::DeltaEta, - matching_params::DeltaPhi, matching_params::DeltaX, matching_params::DeltaY, - matching_params::MchPt, - matching_params::MchEta, - matching_params::MchQ, - matching_params::MftEta, - matching_params::MftQ); - -namespace mix_matching_params -{ -// matching parameters -DECLARE_SOA_COLUMN(NClustMFTTracks, nClustMFT, int); -DECLARE_SOA_COLUMN(Chi2MFTTracks, chi2MFT, float); - -DECLARE_SOA_COLUMN(DeltaP, dp_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPt, dpt_mchplane, float); -DECLARE_SOA_COLUMN(DeltaEta, deta_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPhi, dphi_mchplane, float); -DECLARE_SOA_COLUMN(DeltaX, dx_mchplane, float); -DECLARE_SOA_COLUMN(DeltaY, dy_mchplane, float); - -DECLARE_SOA_COLUMN(MchPt, mchpt, float); -DECLARE_SOA_COLUMN(MchEta, mcheta, float); -DECLARE_SOA_COLUMN(MchPhi, mchphi, float); -DECLARE_SOA_COLUMN(MchQ, mchq, float); - -DECLARE_SOA_COLUMN(MftPt, mftpt, float); -DECLARE_SOA_COLUMN(MftEta, mfteta, float); -DECLARE_SOA_COLUMN(MftPhi, mftphi, float); -DECLARE_SOA_COLUMN(MftQ, mftq, float); -DECLARE_SOA_COLUMN(MftDCA, mftdca, float); -} // namespace mix_matching_params - -DECLARE_SOA_TABLE(MixMatchParams, "AOD", "MIXMATCHING", - mix_matching_params::DeltaPt, - mix_matching_params::DeltaEta, - mix_matching_params::DeltaPhi, - mix_matching_params::DeltaX, - mix_matching_params::DeltaY, - mix_matching_params::MchPt, - mix_matching_params::MchEta, - mix_matching_params::MchQ, - mix_matching_params::MftEta, - mix_matching_params::MftQ); + matching_params::DeltaEta, + matching_params::DeltaPhi, + matching_params::IsCorrectMatch); namespace tag_matching_params { // matching parameters -DECLARE_SOA_COLUMN(DeltaPt, dpt_mchplane, float); -DECLARE_SOA_COLUMN(DeltaEta, deta_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPhi, dphi_mchplane, float); -DECLARE_SOA_COLUMN(DeltaX, dx_mchplane, float); -DECLARE_SOA_COLUMN(DeltaY, dy_mchplane, float); - -DECLARE_SOA_COLUMN(MchPt, mchpt, float); -DECLARE_SOA_COLUMN(MchEta, mcheta, float); -DECLARE_SOA_COLUMN(MchQ, mchq, float); - -DECLARE_SOA_COLUMN(MftEta, mfteta, float); -DECLARE_SOA_COLUMN(MftQ, mftq, float); -DECLARE_SOA_COLUMN(IsTaged, isTaged, bool); +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); } // namespace tag_matching_params DECLARE_SOA_TABLE(TagMatchParams, "AOD", "TAGMATCHING", + tag_matching_params::GMuonPt, + tag_matching_params::GMuonEta, + tag_matching_params::PairQ, tag_matching_params::DeltaPt, - tag_matching_params::DeltaEta, - tag_matching_params::DeltaPhi, tag_matching_params::DeltaX, tag_matching_params::DeltaY, - tag_matching_params::MchPt, - tag_matching_params::MchEta, - tag_matching_params::MchQ, - tag_matching_params::MftEta, - tag_matching_params::MftQ, - tag_matching_params::IsTaged); + tag_matching_params::DeltaEta, + tag_matching_params::DeltaPhi, + tag_matching_params::IsCorrectMatch); namespace probe_matching_params { // matching parameters -DECLARE_SOA_COLUMN(NMFTCandTagMuon, nTagMFT, int); -DECLARE_SOA_COLUMN(NMFTCandProbeMuon, nProbeMFT, int); -DECLARE_SOA_COLUMN(TagMuonP, tagmuonp, float); - -DECLARE_SOA_COLUMN(NClustMFTTracks, nClustMFT, int); -DECLARE_SOA_COLUMN(Chi2MFTTracks, chi2MFT, float); - -DECLARE_SOA_COLUMN(DeltaP, dp_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPt, dpt_mchplane, float); -DECLARE_SOA_COLUMN(DeltaEta, deta_mchplane, float); -DECLARE_SOA_COLUMN(DeltaPhi, dphi_mchplane, float); -DECLARE_SOA_COLUMN(DeltaX, dx_mchplane, float); -DECLARE_SOA_COLUMN(DeltaY, dy_mchplane, float); - -DECLARE_SOA_COLUMN(MchPt, mchpt, float); -DECLARE_SOA_COLUMN(MchEta, mcheta, float); -DECLARE_SOA_COLUMN(MchPhi, mchphi, float); -DECLARE_SOA_COLUMN(MchQ, mchq, float); - -DECLARE_SOA_COLUMN(MftPt, mftpt, float); -DECLARE_SOA_COLUMN(MftEta, mfteta, float); -DECLARE_SOA_COLUMN(MftPhi, mftphi, float); -DECLARE_SOA_COLUMN(MftQ, mftq, float); -DECLARE_SOA_COLUMN(MftDCA, mftdca, float); +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(TagGMuonPt, mTagGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); } // namespace probe_matching_params DECLARE_SOA_TABLE(ProbeMatchParams, "AOD", "PROBEMATCHING", - probe_matching_params::NMFTCandTagMuon, - probe_matching_params::NMFTCandProbeMuon, - probe_matching_params::TagMuonP, + probe_matching_params::TagGMuonPt, + probe_matching_params::GMuonPt, + probe_matching_params::GMuonEta, + probe_matching_params::PairQ, probe_matching_params::DeltaPt, - probe_matching_params::DeltaEta, - probe_matching_params::DeltaPhi, probe_matching_params::DeltaX, probe_matching_params::DeltaY, - probe_matching_params::MchPt, - probe_matching_params::MchEta, - probe_matching_params::MchQ, - probe_matching_params::MftEta, - probe_matching_params::MftQ); + probe_matching_params::DeltaEta, + probe_matching_params::DeltaPhi, + probe_matching_params::IsCorrectMatch); + +namespace mix_matching_params +{ +// matching parameters +DECLARE_SOA_COLUMN(DeltaPt, mDeltaPt, float); +DECLARE_SOA_COLUMN(DeltaEta, mDeltaEta, float); +DECLARE_SOA_COLUMN(DeltaPhi, mDeltaPhi, float); +DECLARE_SOA_COLUMN(DeltaX, mDeltaX, float); +DECLARE_SOA_COLUMN(DeltaY, mDeltaY, float); +DECLARE_SOA_COLUMN(GMuonPt, mGMuonPt, float); +DECLARE_SOA_COLUMN(GMuonEta, mGMuonEta, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); +DECLARE_SOA_COLUMN(IsCorrectMatch, mIsCorrectMatch, bool); +} // namespace mix_matching_params + +DECLARE_SOA_TABLE(MixMatchParams, "AOD", "MIXMATCHING", + mix_matching_params::GMuonPt, + mix_matching_params::GMuonEta, + mix_matching_params::PairQ, + mix_matching_params::DeltaPt, + mix_matching_params::DeltaX, + mix_matching_params::DeltaY, + mix_matching_params::DeltaEta, + mix_matching_params::DeltaPhi, + mix_matching_params::IsCorrectMatch); namespace muon_pair { -DECLARE_SOA_COLUMN(NMFT, nMft, int); -DECLARE_SOA_COLUMN(Q, q, int16_t); -DECLARE_SOA_COLUMN(M, m, float); -DECLARE_SOA_COLUMN(Pt, pt, float); -DECLARE_SOA_COLUMN(Rap, rap, float); +// matching parameters +DECLARE_SOA_COLUMN(Mass, mMass, float); +DECLARE_SOA_COLUMN(Pt, mPt, float); +DECLARE_SOA_COLUMN(Rap, mRap, float); +DECLARE_SOA_COLUMN(PairQ, mPairQ, int16_t); } // namespace muon_pair -DECLARE_SOA_TABLE(MuonPair, "AOD", "DIMUON", muon_pair::NMFT, muon_pair::Q, muon_pair::M, muon_pair::Pt, muon_pair::Rap); +DECLARE_SOA_TABLE(MuonPair, "AOD", "MUONPAIR", + muon_pair::PairQ, + muon_pair::Mass, + muon_pair::Pt, + muon_pair::Rap); -namespace tag_muon_pair -{ -DECLARE_SOA_COLUMN(NMFT, nMft, int); -DECLARE_SOA_COLUMN(Q, q, int16_t); -DECLARE_SOA_COLUMN(M, m, float); -DECLARE_SOA_COLUMN(Pt, pt, float); -DECLARE_SOA_COLUMN(Rap, rap, float); -} // namespace tag_muon_pair - -DECLARE_SOA_TABLE(TagMuonPair, "AOD", "TAGDIMUON", tag_muon_pair::NMFT, tag_muon_pair::Q, tag_muon_pair::M, tag_muon_pair::Pt, tag_muon_pair::Rap); } // namespace o2::aod - +*/ struct match_mft_mch_data { - Produces matchingParams; - Produces tagmatchingParams; - Produces probematchingParams; - Produces mixmatchingParams; - Produces muonPairs; - Produces tagmuonPairs; - Produces muonParams; - Produces mftParams; - - HistogramRegistry registry{ - "registry", - {{"hMchP", "MCH track total momentum (at the first station); p [GeV/c]; Counts", {HistType::kTH1F, {{2000, 0, 200}}}}, - {"hMchCorrP", "MCH track total momentum (propagated to PV); p [GeV/c]; Counts", {HistType::kTH1F, {{2000, 0, 200}}}}, - {"hMassCorrMchPair", "Corrected MCH track pair mass (propagated to PV); m [GeV/c^{2}]; Counts", {HistType::kTH1F, {{1000, 0, 10}}}}}}; + //// Variables for matching method + Configurable fMatchingMethod{"cfgMatchingMethod", 0, ""}; //// Variables for selecting muon tracks Configurable fEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; @@ -356,8 +248,8 @@ struct match_mft_mch_data { Configurable fTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; /// Variables to add preselection for the matching table - Configurable fPreselectMatchingX{"cfgPreselectMatchingX", 999.f, ""}; - Configurable fPreselectMatchingY{"cfgPreselectMatchingY", 999.f, ""}; + Configurable fPreselectMatchingX{"cfgPreselectMatchingX", 15.f, ""}; + Configurable fPreselectMatchingY{"cfgPreselectMatchingY", 15.f, ""}; /// Variables to event mixing criteria Configurable fSaveMixedMatchingParamsRate{"cfgSaveMixedMatchingParamsRate", 0.002f, ""}; @@ -367,504 +259,704 @@ struct match_mft_mch_data { //// Variables for selecting tag muon Configurable fTagMassWindowMin{"cfgTagMassWindowMin", 2.8f, ""}; Configurable fTagMassWindowMax{"cfgTagMassWindowMax", 3.3f, ""}; - Configurable fTagXWindowLow{"cfgTagXWindowLow", -0.80f, ""}; - Configurable fTagXWindowUp{"cfgTagXWindowUp", 0.84f, ""}; - Configurable fTagYWindowLow{"cfgTagYWindowLow", -0.64f, ""}; - Configurable fTagYWindowUp{"cfgTagYWindowUp", 0.95f, ""}; - Configurable fTagPhiWindowLow{"cfgTagPhiWindowLow", -0.041f, ""}; - Configurable fTagPhiWindowUp{"cfgTagPhiWindowUp", 0.039f, ""}; - Configurable fTagEtaWindowLow{"cfgTagEtaWindowLow", -0.045f, ""}; - Configurable fTagEtaWindowUp{"cfgTagEtaWindowUp", 0.033f, ""}; - - //// Variables for selecting probe muon - Configurable fProbeXWindowLow{"cfgProbeXWindowLow", -10.f, ""}; - Configurable fProbeXWindowUp{"cfgProbeXWindowUp", 10.f, ""}; - Configurable fProbeYWindowLow{"cfgProbeYWindowLow", -10.f, ""}; - Configurable fProbeYWindowUp{"cfgProbeYWindowUp", 10.f, ""}; - Configurable fProbePhiWindowLow{"cfgProbePhiWindowLow", -1.0f, ""}; - Configurable fProbePhiWindowUp{"cfgProbePhiWindowUp", 1.0f, ""}; - Configurable fProbeEtaWindowLow{"cfgProbeEtaWindowLow", -1.0f, ""}; - Configurable fProbeEtaWindowUp{"cfgProbeEtaWindowUp", 1.0f, ""}; - - Service ccdb; + + //// Variables for ccdb Configurable ccdburl{"ccdb-url", "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"}; - // o2::parameters::GRPMagField* grpmag = nullptr; - o2::globaltracking::MatchGlobalFwd mMatching; - o2::field::MagneticField* fieldB; + //// Variables for Tag matching criteria + Configurable fSigmaXTagMuonCut{"cfgSigmaXTagMuonCut", 1.f, ""}; + Configurable fMeanXTagMuonCut{"cfgMeanXTagMuonCut", 0.f, ""}; + Configurable fSigmaYTagMuonCut{"cfgSigmaYTagMuonCut", 1.f, ""}; + Configurable fMeanYTagMuonCut{"cfgMeanYTagMuonCut", 1.f, ""}; - o2::ccdb::CcdbApi ccdbApi; - int mRunNumber; + Configurable fSigmaEtaTagMuonCut{"cfgSigmaEtaTagMuonCut", 0.2f, ""}; + Configurable fMeanEtaTagMuonCut{"cfgMeanEtaTagMuonCut", 0.f, ""}; + Configurable fSigmaPhiTagMuonCut{"cfgSigmaPhiTagMuonCut", 0.2f, ""}; + Configurable fMeanPhiTagMuonCut{"cfgMeanPhiTagMuonCut", 0.f, ""}; - void init(o2::framework::InitContext&) + template + class FindTagAndProbe { - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - ccdbApi.init(ccdburl); - mRunNumber = 0; - } + private: + o2::dataformats::GlobalFwdTrack muontrack_at_pv[2]; + + TLorentzVector mDimuon; + MUON muontrack1; + MUON muontrack2; + Collision collision; + int tagIdx, probeIdx; + + int16_t mQ; + + inline void fillCovarianceArray(MUON const& muontrack, float cov[15]) const + { + cov[0] = muontrack.cXX(); + cov[1] = muontrack.cXY(); + cov[2] = muontrack.cYY(); + cov[3] = muontrack.cPhiX(); + cov[4] = muontrack.cPhiY(); + cov[5] = muontrack.cPhiPhi(); + cov[6] = muontrack.cTglX(); + cov[7] = muontrack.cTglY(); + cov[8] = muontrack.cTglPhi(); + cov[9] = muontrack.cTglTgl(); + cov[10] = muontrack.c1PtX(); + cov[11] = muontrack.c1PtY(); + cov[12] = muontrack.c1PtPhi(); + cov[13] = muontrack.c1PtTgl(); + cov[14] = muontrack.c1Pt21Pt2(); + } - void initCCDB(ExtBCs::iterator const& bc) - { - if (mRunNumber == bc.runNumber()) { - return; + inline o2::dataformats::GlobalFwdTrack propagateMUONtoPV(MUON const& muontrack) const + { + const double mz = muontrack.z(); + const double mchi2 = muontrack.chi2(); + const float mx = muontrack.x(); + const float my = muontrack.y(); + const float mphi = muontrack.phi(); + const float mtgl = muontrack.tgl(); + const float m1pt = muontrack.signed1Pt(); + + float cov[15]; + fillCovarianceArray(muontrack, cov); + SMatrix5 tpars(mx, my, mphi, mtgl, m1pt); + SMatrix55 tcovs(cov, cov + 15); + + o2::track::TrackParCovFwd parcovmuontrack{mz, tpars, tcovs, mchi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + + o2::globaltracking::MatchGlobalFwd mMatching; + auto mchtrack = mMatching.FwdtoMCH(gtrack); + + o2::mch::TrackExtrap::extrapToVertex(mchtrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + o2::dataformats::GlobalFwdTrack extrap_muontrack; + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + + return extrap_muontrack; } - mRunNumber = bc.runNumber(); - std::map metadata; - auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, mRunNumber); - auto ts = soreor.first; - auto grpmag = ccdbApi.retrieveFromTFileAny(grpmagPath, metadata, ts); - o2::base::Propagator::initFieldFromGRP(grpmag); - if (!o2::base::GeometryManager::isGeometryLoaded()) { - ccdb->get(geoPath); + + inline void setTagAndProbe() + { + if (muontrack1.pt() > muontrack2.pt()) { + tagIdx = 0; + probeIdx = 1; + } else { + tagIdx = 1; + probeIdx = 0; + } } - o2::mch::TrackExtrap::setField(); - fieldB = static_cast(TGeoGlobalMagField::Instance()->GetField()); - } - enum ProagationPoint { ToVtx, - ToDCA }; + public: + inline FindTagAndProbe(const MUON& muon1, const MUON& muon2, const Collision& coll) + : muontrack_at_pv(), mDimuon(), muontrack1(muon1), muontrack2(muon2), collision(coll), tagIdx(-1), probeIdx(-1), mQ(0) + { + mQ = muontrack1.sign() + muontrack2.sign(); + setTagAndProbe(); + } - template - o2::dataformats::GlobalFwdTrack PropagateMUONTrack(T const& muon, int PropType) + void calcMuonPairAtPV() + { + muontrack_at_pv[0] = propagateMUONtoPV(muontrack1); + muontrack_at_pv[1] = propagateMUONtoPV(muontrack2); + TLorentzVector vMuon1, vMuon2; + vMuon1.SetPtEtaPhiM(muontrack_at_pv[0].getPt(), muontrack_at_pv[0].getEta(), muontrack_at_pv[0].getPhi(), mMu); + vMuon2.SetPtEtaPhiM(muontrack_at_pv[1].getPt(), muontrack_at_pv[1].getEta(), muontrack_at_pv[1].getPhi(), mMu); + mDimuon = vMuon1 + vMuon2; + } + inline int getTagMuonIndex() const { return tagIdx; } + inline int getProbeMuonIndex() const { return probeIdx; } + inline float getMass() const { return mDimuon.M(); } + inline float getPt() const { return mDimuon.Pt(); } + inline float getRap() const { return mDimuon.Rapidity(); } + inline int16_t getCharge() const { return mQ; } + inline const o2::dataformats::GlobalFwdTrack& getMuonAtPV(int idx) const { return muontrack_at_pv[idx]; } + }; // end of class FindTagAndProbe + + template + class MatchingParamsML { + private: + MUON muontrack; + MFT mfttrack; + Collision collision; + + float mDX, mDY, mDPt, mDPhi, mDEta; + float mGlobalMuonPtAtDCA, mGlobalMuonEtaAtDCA, mGlobalMuonPhiAtDCA, mGlobalMuonDCAx, mGlobalMuonDCAy, mGlobalMuonQ; + int mMatchingType; + + o2::field::MagneticField* fieldB; + o2::globaltracking::MatchGlobalFwd mMatching; + + inline o2::track::TrackParCovFwd propagateMFTtoMatchingPlane() + { + double covArr[15]{0.0}; + SMatrix55 tmftcovs(covArr, covArr + 15); + + SMatrix5 tmftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); + o2::track::TrackParCovFwd extrap_mfttrack{mfttrack.z(), tmftpars, tmftcovs, mfttrack.chi2()}; + + double propVec[3] = {0.}; + float zPlane = 0.f; + if (mMatchingType == MCH_FIRST_CLUSTER) { + propVec[0] = muontrack.x() - mfttrack.x(); + propVec[1] = muontrack.y() - mfttrack.y(); + propVec[2] = muontrack.z() - mfttrack.z(); + zPlane = muontrack.z(); + } else if (mMatchingType == END_OF_ABSORBER || mMatchingType == BEGINING_OF_ABSORBER) { + auto extrap_muontrack = propagateMUONtoMatchingPlane(); + propVec[0] = extrap_muontrack.getX() - mfttrack.x(); + propVec[1] = extrap_muontrack.getY() - mfttrack.y(); + propVec[2] = extrap_muontrack.getZ() - mfttrack.z(); + zPlane = (mMatchingType == END_OF_ABSORBER) ? -505.f : -90.f; + } else { + zPlane = mfttrack.z(); + } - auto collision = muon.collision(); - - o2::dataformats::GlobalFwdTrack propmuon; - - double chi2 = muon.chi2(); - - SMatrix5 tpars(muon.x(), muon.y(), muon.phi(), muon.tgl(), muon.signed1Pt()); - std::vector v1{muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), - muon.cPhiPhi(), muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), - muon.c1PtX(), muon.c1PtY(), muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2()}; - if (isGoodMUONTrack(muon)) { - SMatrix55 tcovs(v1.begin(), v1.end()); - o2::track::TrackParCovFwd fwdtrack{muon.z(), tpars, tcovs, chi2}; - - o2::dataformats::GlobalFwdTrack track; - track.setParameters(tpars); - track.setZ(fwdtrack.getZ()); - track.setCovariances(tcovs); - auto mchTrack = mMatching.FwdtoMCH(track); - if (PropType == ProagationPoint::ToVtx) - o2::mch::TrackExtrap::extrapToVertex(mchTrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); - else if (PropType == ProagationPoint::ToDCA) - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, collision.posZ()); - - auto proptrack = mMatching.MCHtoFwd(mchTrack); - propmuon.setParameters(proptrack.getParameters()); - propmuon.setZ(proptrack.getZ()); - propmuon.setCovariances(proptrack.getCovariances()); + double centerZ[3] = {mfttrack.x() + propVec[0] / 2., mfttrack.y() + propVec[1] / 2., mfttrack.z() + propVec[2] / 2.}; + float Bz = fieldB->getBz(centerZ); + extrap_mfttrack.propagateToZ(zPlane, Bz); // z in cm + return extrap_mfttrack; } - v1.clear(); - v1.shrink_to_fit(); + inline o2::dataformats::GlobalFwdTrack propagateMUONtoMatchingPlane() + { + float cov[15] = { + muontrack.cXX(), muontrack.cXY(), muontrack.cYY(), + muontrack.cPhiX(), muontrack.cPhiY(), muontrack.cPhiPhi(), + muontrack.cTglX(), muontrack.cTglY(), muontrack.cTglPhi(), + muontrack.cTglTgl(), muontrack.c1PtX(), muontrack.c1PtY(), + muontrack.c1PtPhi(), muontrack.c1PtTgl(), muontrack.c1Pt21Pt2()}; + + SMatrix5 tpars(muontrack.x(), muontrack.y(), muontrack.phi(), muontrack.tgl(), muontrack.signed1Pt()); + SMatrix55 tcovs(cov, cov + 15); + double chi2 = muontrack.chi2(); + + o2::track::TrackParCovFwd parcovmuontrack{muontrack.z(), tpars, tcovs, chi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + + auto mchtrack = mMatching.FwdtoMCH(gtrack); + + if (mMatchingType == MFT_LAST_CLUSTR) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchtrack, mfttrack.z()); + } else if (mMatchingType == END_OF_ABSORBER) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchtrack, -505.); + } else if (mMatchingType == BEGINING_OF_ABSORBER) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchtrack, -90.); + } - return propmuon; - } + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); - template - o2::track::TrackParCovFwd PropagateMFT(T const& mfttrack, int PropType) - { - std::vector mftv1; - SMatrix55 mftcovs{mftv1.begin(), mftv1.end()}; - SMatrix5 mftpars = {mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()}; - o2::track::TrackParCovFwd mftpartrack = {mfttrack.z(), mftpars, mftcovs, mfttrack.chi2()}; - if (PropType == ProagationPoint::ToDCA) { - auto collision = mfttrack.collision(); - double propVec[3] = {fabs(mfttrack.x() - collision.posX()), fabs(mfttrack.y() - collision.posY()), fabs(mfttrack.z() - collision.posZ())}; - double centerZ[3] = {mfttrack.x() - propVec[0] / 2., mfttrack.y() - propVec[1] / 2., mfttrack.z() - propVec[2] / 2.}; + o2::dataformats::GlobalFwdTrack extrap_muontrack; + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + return extrap_muontrack; + } + + inline o2::track::TrackParCovFwd propagateMFTtoDCA() + { + double covArr[15]{0.0}; + SMatrix55 tmftcovs(covArr, covArr + 15); + + SMatrix5 tmftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); + o2::track::TrackParCovFwd extrap_mfttrack{mfttrack.z(), tmftpars, tmftcovs, mfttrack.chi2()}; + + double propVec[3] = {}; + propVec[0] = collision.posX() - mfttrack.x(); + propVec[1] = collision.posY() - mfttrack.y(); + propVec[2] = collision.posZ() - mfttrack.z(); + + double centerZ[3] = {mfttrack.x() + propVec[0] / 2., mfttrack.y() + propVec[1] / 2., mfttrack.z() + propVec[2] / 2.}; float Bz = fieldB->getBz(centerZ); - mftpartrack.propagateToZ(collision.posZ(), Bz); + extrap_mfttrack.propagateToZ(collision.posZ(), Bz); // z in cm + return extrap_mfttrack; } - return mftpartrack; + + inline o2::dataformats::GlobalFwdTrack propagateMUONtoPV() + { + float cov[15] = { + muontrack.cXX(), muontrack.cXY(), muontrack.cYY(), + muontrack.cPhiX(), muontrack.cPhiY(), muontrack.cPhiPhi(), + muontrack.cTglX(), muontrack.cTglY(), muontrack.cTglPhi(), + muontrack.cTglTgl(), muontrack.c1PtX(), muontrack.c1PtY(), + muontrack.c1PtPhi(), muontrack.c1PtTgl(), muontrack.c1Pt21Pt2()}; + + SMatrix5 tpars(muontrack.x(), muontrack.y(), muontrack.phi(), muontrack.tgl(), muontrack.signed1Pt()); + SMatrix55 tcovs(cov, cov + 15); + double chi2 = muontrack.chi2(); + + o2::track::TrackParCovFwd parcovmuontrack{muontrack.z(), tpars, tcovs, chi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + + auto mchtrack = mMatching.FwdtoMCH(gtrack); + o2::mch::TrackExtrap::extrapToVertex(mchtrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + o2::dataformats::GlobalFwdTrack extrap_muontrack; + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + + return extrap_muontrack; + } + + public: + enum MATCHING_TYPE { MCH_FIRST_CLUSTER, + MFT_LAST_CLUSTR, + END_OF_ABSORBER, + BEGINING_OF_ABSORBER }; + + MatchingParamsML(MUON const& muon, MFT const& mft, Collision const& coll, int MType, o2::field::MagneticField* field) : muontrack(muon), mfttrack(mft), collision(coll), mDX(0.f), mDY(0.f), mDPt(0.f), mDPhi(0.f), mDEta(0.f), mGlobalMuonPtAtDCA(0.f), mGlobalMuonEtaAtDCA(0.f), mGlobalMuonPhiAtDCA(0.f), mGlobalMuonDCAx(0.f), mGlobalMuonDCAy(0.f), mGlobalMuonQ(0.f), mMatchingType(MType), fieldB(field) {} + void calcMatchingParams() + { + auto mfttrack_on_matchingP = propagateMFTtoMatchingPlane(); + auto muontrack_on_matchingP = propagateMUONtoMatchingPlane(); + + float dphiRaw = mfttrack_on_matchingP.getPhi() - muontrack_on_matchingP.getPhi(); + float dphi = TVector2::Phi_mpi_pi(dphiRaw); + float deta = mfttrack_on_matchingP.getEta() - muontrack_on_matchingP.getEta(); + + mDX = mfttrack_on_matchingP.getX() - muontrack_on_matchingP.getX(); + mDY = mfttrack_on_matchingP.getY() - muontrack_on_matchingP.getY(); + mDPt = mfttrack_on_matchingP.getPt() - muontrack_on_matchingP.getPt(); + mDPhi = dphi; + mDEta = deta; + } + + void calcGlobalMuonParams() + { + auto mfttrack_at_dca = propagateMFTtoDCA(); + auto muontrack_at_pv = propagateMUONtoPV(); + + float momentum = muontrack_at_pv.getP(); + float theta = mfttrack_at_dca.getTheta(); + float phiTrack = mfttrack_at_dca.getPhi(); + float px = momentum * std::sin(theta) * std::cos(phiTrack); + float py = momentum * std::sin(theta) * std::sin(phiTrack); + + mGlobalMuonQ = muontrack.sign() + mfttrack.sign(); + mGlobalMuonPtAtDCA = std::sqrt(px * px + py * py); + mGlobalMuonEtaAtDCA = mfttrack_at_dca.getEta(); + mGlobalMuonPhiAtDCA = mfttrack_at_dca.getPhi(); + mGlobalMuonDCAx = mfttrack_at_dca.getX() - collision.posX(); + mGlobalMuonDCAy = mfttrack_at_dca.getY() - collision.posY(); + } + + inline float getDx() const { return mDX; } + inline float getDy() const { return mDY; } + inline float getDphi() const { return mDPhi; } + inline float getDeta() const { return mDEta; } + inline float getDpt() const { return mDPt; } + inline float getGMPtAtDCA() const { return mGlobalMuonPtAtDCA; } + inline float getGMEtaAtDCA() const { return mGlobalMuonEtaAtDCA; } + inline float getGMPhiAtDCA() const { return mGlobalMuonPhiAtDCA; } + inline float getGMDcaX() const { return mGlobalMuonDCAx; } + inline float getGMDcaY() const { return mGlobalMuonDCAy; } + inline float getGMDcaXY() const { return std::sqrt(mGlobalMuonDCAx * mGlobalMuonDCAx + mGlobalMuonDCAy * mGlobalMuonDCAy); } + inline int16_t getGMQ() const { return static_cast(mGlobalMuonQ); } + + }; // end of class MatchingParamsML + + template + o2::dataformats::GlobalFwdTrack propagateMUONtoPV(MUON const& muontrack, Collisions const& collisions) + { + auto collision = collisions.rawIteratorAt(muontrack.collisionId()); + o2::globaltracking::MatchGlobalFwd mMatching; + o2::dataformats::GlobalFwdTrack extrap_muontrack; + + SMatrix5 tpars(muontrack.x(), muontrack.y(), muontrack.phi(), muontrack.tgl(), muontrack.signed1Pt()); + std::vector v1{muontrack.cXX(), muontrack.cXY(), muontrack.cYY(), + muontrack.cPhiX(), muontrack.cPhiY(), muontrack.cPhiPhi(), + muontrack.cTglX(), muontrack.cTglY(), muontrack.cTglPhi(), + muontrack.cTglTgl(), muontrack.c1PtX(), muontrack.c1PtY(), + muontrack.c1PtPhi(), muontrack.c1PtTgl(), muontrack.c1Pt21Pt2()}; + SMatrix55 tcovs(v1.begin(), v1.end()); + double chi2 = muontrack.chi2(); + o2::track::TrackParCovFwd parcovmuontrack{muontrack.z(), tpars, tcovs, chi2}; + + o2::dataformats::GlobalFwdTrack gtrack; + gtrack.setParameters(tpars); + gtrack.setZ(parcovmuontrack.getZ()); + gtrack.setCovariances(tcovs); + auto mchtrack = mMatching.FwdtoMCH(gtrack); + + o2::mch::TrackExtrap::extrapToVertex(mchtrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + auto fwdtrack = mMatching.MCHtoFwd(mchtrack); + extrap_muontrack.setParameters(fwdtrack.getParameters()); + extrap_muontrack.setZ(fwdtrack.getZ()); + extrap_muontrack.setCovariances(fwdtrack.getCovariances()); + + return extrap_muontrack; } - template - o2::track::TrackParCovFwd PropagateMFTtoMatchingPlane(MFT const& mfttrack, FWD const& fwdtrack) + inline bool isGoodTagDimuon(float M) { - std::vector v1; // Temporary null vector for the computation of the covariance matrix - double propVec[3] = {fwdtrack.x() - mfttrack.x(), fwdtrack.y() - mfttrack.y(), fwdtrack.z() - mfttrack.z()}; - double centerZ[3] = {mfttrack.x() + propVec[0] / 2., mfttrack.y() + propVec[1] / 2., mfttrack.z() + propVec[2] / 2.}; - float Bz = fieldB->getBz(centerZ); // gives error if the propagator is not initFielded - SMatrix55 tmftcovs(v1.begin(), v1.end()); - SMatrix5 tmftpars(mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()); - o2::track::TrackParCovFwd extrap_mfttrack{mfttrack.z(), tmftpars, tmftcovs, mfttrack.chi2()}; - extrap_mfttrack.propagateToZ(fwdtrack.z(), Bz); // z in cm - return extrap_mfttrack; + return !(M < fTagMassWindowMin || M > fTagMassWindowMax); } - template - bool isGoodMUONTrack(T track) + inline bool isGoodTagMatching(float mDX, float mDY, float mDEta, float mDPhi) + { + float dxNorm = (mDX - fMeanXTagMuonCut) / (fSigmaXTagMuonCut * 3); + float dyNorm = (mDY - fMeanYTagMuonCut) / (fSigmaYTagMuonCut * 3); + float detaNorm = (mDEta - fMeanEtaTagMuonCut) / (fSigmaEtaTagMuonCut * 3); + float dphiNorm = (mDPhi - fMeanPhiTagMuonCut) / (fSigmaPhiTagMuonCut * 3); + + float rTagXY = dxNorm * dxNorm + dyNorm * dyNorm; + float rTagEtaPhi = detaNorm * detaNorm + dphiNorm * dphiNorm; + + return (rTagXY < 1.f && rTagEtaPhi > 0.f); + } + + template + bool isCorrectMatching(MUON const& muontrack, MFT const& mfttrack) + { + + int idmuon = muontrack.mcParticleId(); + int idmft = mfttrack.mcParticleId(); + + if (idmuon == -1 || idmft == -1) + return false; + if (idmuon != idmft) + return false; + else + return true; + }; + + template + bool isGoodMuonQuality(MUON muontrack) { - if (!track.has_collision()) + if (!muontrack.has_collision()) return false; - if (track.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) + if (muontrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) return false; - if (track.chi2() > fTrackChi2MchUp) + if (muontrack.chi2() > fTrackChi2MchUp) return false; - if (fRabsLow1 > track.rAtAbsorberEnd() || track.rAtAbsorberEnd() > fRabsUp2) + if (fRabsLow1 > muontrack.rAtAbsorberEnd() || muontrack.rAtAbsorberEnd() > fRabsUp2) return false; - if (track.rAtAbsorberEnd() < fRabsUp1 && fPdcaUp1 < track.pDca()) + if (muontrack.rAtAbsorberEnd() < fRabsUp1 && fPdcaUp1 < muontrack.pDca()) return false; - if (track.rAtAbsorberEnd() > fRabsLow2 && fPdcaUp2 < track.pDca()) + if (muontrack.rAtAbsorberEnd() > fRabsLow2 && fPdcaUp2 < muontrack.pDca()) return false; return true; } - template - int selectTagMUON(T track1, T track2) + template + bool isGoodMuonKine(MUON muontrack) { - if (track1.pt() > track2.pt()) { - return track1.globalIndex(); - } else { - return track2.globalIndex(); - } - } - - template - int selectProbeMUONTrack(T track1, T track2) - { - if (track1.pt() < track2.pt()) { - return track1.globalIndex(); - } else { - return track2.globalIndex(); - } - } - - bool isGoodKenematicMUONTrack(o2::dataformats::GlobalFwdTrack track) - { - if (fEtaMchLow > track.getEta() || track.getEta() > fEtaMchUp) + if (fEtaMchLow > muontrack.getEta() || muontrack.getEta() > fEtaMchUp) return false; return true; } - template - bool isGoodMFTTrack(T track) + template + bool isGoodMFTQuality(MFT mfttrack) { - if (!track.has_collision()) + if (!mfttrack.has_collision()) return false; - if (track.chi2() > fTrackChi2MftUp) + if (mfttrack.chi2() > fTrackChi2MftUp) return false; - if (track.nClusters() < fTrackNClustMftLow) + if (mfttrack.nClusters() < fTrackNClustMftLow) return false; return true; } - bool isGoodKenematicMFTTrack(o2::track::TrackParCovFwd track) + template + bool isGoodMFTKine(MFT mfttrack) { - if (fEtaMftLow > track.getEta() || track.getEta() > fEtaMftUp) { + if (fEtaMftLow > mfttrack.getEta() || mfttrack.getEta() > fEtaMftUp) return false; - } return true; } - void process(aod::Collisions const& collisions, ExtBCs const& ebcs, - MyMuons const& fwdtracks, MyMFTs const& mfttracks) + inline bool isPassMatchingPreselection(float Dx, float Dy) { - initCCDB(ebcs.begin()); - - std::unordered_set bcs_mfttrack; - std::unordered_map map_vtxZ; - std::unordered_map nmfttracks; - std::unordered_map> map_mfttraks; - - for (const auto& mfttrack : mfttracks) { + return !(std::abs(Dx) > fPreselectMatchingX || std::abs(Dy) > fPreselectMatchingY); + } - if (!isGoodMFTTrack(mfttrack)) { + template + void setMUONs(MUONs const& muontracks, Collisions const& collisions) + { + for (auto muontrack : muontracks) { + if (!isGoodMuonQuality(muontrack)) continue; - } - - bcs_mfttrack.insert(mfttrack.collisionId()); - std::vector& tracks = map_mfttraks[mfttrack.collisionId()]; - tracks.push_back(mfttrack.globalIndex()); - - o2::track::TrackParCovFwd mftpartrack = PropagateMFT(mfttrack, ProagationPoint::ToDCA); - - if (!isGoodKenematicMFTTrack(mftpartrack)) { + o2::dataformats::GlobalFwdTrack muontrack_at_pv = propagateMUONtoPV(muontrack, collisions); + if (!isGoodMuonKine(muontrack_at_pv)) continue; - } - auto collision = mfttrack.collision(); + auto collision = collisions.rawIteratorAt(muontrack.collisionId()); - map_vtxZ[mfttrack.collisionId()] = collision.posZ(); + bool& has = map_has_muontracks_collisions[muontrack.collisionId()]; + has = true; - float dx = mftpartrack.getX() - collision.posX(); - float dy = mftpartrack.getY() - collision.posY(); - - mftParams(mfttrack.nClusters(), mfttrack.isCA(), mfttrack.chi2(), mfttrack.sign(), mftpartrack.getPt(), mftpartrack.getEta(), mftpartrack.getPhi(), dx, dy); - nmfttracks[mfttrack.collisionId()]++; + vector& arr_muontracks = map_muontracks[collision.globalIndex()]; + arr_muontracks.push_back(muontrack.globalIndex()); } + } - std::unordered_map nfwdtracks; - std::unordered_map> map_fwdtraks; - - for (auto fwdtrack : fwdtracks) { + template + o2::track::TrackParCovFwd PropagateMFTtoDCA(MFT const& mfttrack, Collisions const& collisions, o2::field::MagneticField* field) + { + auto collision = collisions.rawIteratorAt(mfttrack.collisionId()); + std::vector mftv1; + SMatrix55 mftcovs{mftv1.begin(), mftv1.end()}; + SMatrix5 mftpars = {mfttrack.x(), mfttrack.y(), mfttrack.phi(), mfttrack.tgl(), mfttrack.signed1Pt()}; + o2::track::TrackParCovFwd mftpartrack = {mfttrack.z(), mftpars, mftcovs, mfttrack.chi2()}; + double propVec[3] = {fabs(mfttrack.x() - collision.posX()), + fabs(mfttrack.y() - collision.posY()), + fabs(mfttrack.z() - collision.posZ())}; + double centerZ[3] = {mfttrack.x() - propVec[0] / 2., + mfttrack.y() - propVec[1] / 2., + mfttrack.z() - propVec[2] / 2.}; + float Bz = field->getBz(centerZ); + mftpartrack.propagateToZ(collision.posZ(), Bz); + return mftpartrack; + } - if (!isGoodMUONTrack(fwdtrack)) { + template + void setMFTs(MFTs const& mfttracks, Collisions const& collisions, o2::field::MagneticField* field) + { + for (auto mfttrack : mfttracks) { + if (!isGoodMFTQuality(mfttrack)) continue; - } - o2::dataformats::GlobalFwdTrack propmuonAtPV = PropagateMUONTrack(fwdtrack, ProagationPoint::ToVtx); - if (!isGoodKenematicMUONTrack(propmuonAtPV)) { + o2::track::TrackParCovFwd mfttrack_at_dca = PropagateMFTtoDCA(mfttrack, collisions, field); + if (!isGoodMFTKine(mfttrack_at_dca)) continue; - } - - std::vector& tracks = map_fwdtraks[fwdtrack.collisionId()]; - tracks.push_back(fwdtrack.globalIndex()); - bool hasMFT = false; + auto collision = collisions.rawIteratorAt(mfttrack.collisionId()); - std::vector& mfttracks = map_mfttraks[fwdtrack.collisionId()]; + map_vtxz[mfttrack.collisionId()] = collision.posZ(); + map_nmfttrack[mfttrack.collisionId()] += 1; - if (mfttracks.size() > 0) { - hasMFT = true; - } + bool& has = map_has_mfttracks_collisions[mfttrack.collisionId()]; + has = true; - muonParams(fwdtrack.chi2(), fwdtrack.rAtAbsorberEnd(), fwdtrack.sign(), propmuonAtPV.getPt(), propmuonAtPV.getEta(), propmuonAtPV.getPhi(), hasMFT); - nfwdtracks[fwdtrack.collisionId()]++; + vector& arr_mfttracks = map_mfttracks[collision.globalIndex()]; + arr_mfttracks.push_back(mfttrack.globalIndex()); } + } - for (auto fwdtrack1 : fwdtracks) { - - if (!isGoodMUONTrack(fwdtrack1)) { - continue; - } - - int ibc = fwdtrack1.collisionId(); - auto collision = fwdtrack1.collision(); - - o2::dataformats::GlobalFwdTrack fwdtrackAtPV1 = PropagateMUONTrack(fwdtrack1, ProagationPoint::ToVtx); - if (!isGoodKenematicMUONTrack(fwdtrackAtPV1)) { - continue; - } - - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////// MIXED EVENT /////////////// - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - for (auto bc_mfttrack : bcs_mfttrack) { - - if (ibc == bc_mfttrack) - continue; - if (fabs(nmfttracks[ibc] - nmfttracks[bc_mfttrack]) > fEventMaxDeltaNMFT) { - continue; - } - if (fabs(map_vtxZ[bc_mfttrack] - collision.posZ()) > fEventMaxDeltaVtxZ) { - continue; - } - - std::vector& mfttrackGlobalIndex = map_mfttraks[bc_mfttrack]; - - for (int idmfttrack1 = 0; idmfttrack1 < static_cast(mfttrackGlobalIndex.size()); ++idmfttrack1) { - - auto mfttrack1 = mfttracks.rawIteratorAt(mfttrackGlobalIndex[idmfttrack1]); - o2::track::TrackParCovFwd mfttrack_at_matching = PropagateMFTtoMatchingPlane(mfttrack1, fwdtrack1); - - V1.SetPtEtaPhi(mfttrack_at_matching.getPt(), mfttrack_at_matching.getEta(), mfttrack_at_matching.getPhi()); - V2.SetPtEtaPhi(fwdtrack1.pt(), fwdtrack1.eta(), fwdtrack1.phi()); - - double deltaPt = mfttrack_at_matching.getPt() - fwdtrack1.pt(); - double deltaX = mfttrack_at_matching.getX() - fwdtrack1.x(); - double deltaY = mfttrack_at_matching.getY() - fwdtrack1.y(); - double deltaPhi = V1.DeltaPhi(V2); - double deltaEta = mfttrack_at_matching.getEta() - fwdtrack1.eta(); - - if (fabs(deltaX) > fPreselectMatchingX) { - continue; - } - if (fabs(deltaY) > fPreselectMatchingY) { - continue; - } - - o2::track::TrackParCovFwd mfttrack_at_dca = PropagateMFT(mfttrack1, ProagationPoint::ToDCA); - - mixmatchingParams(deltaPt, deltaEta, deltaPhi, deltaX, deltaY, fwdtrackAtPV1.getPt(), fwdtrackAtPV1.getEta(), fwdtrack1.sign(), mfttrack_at_dca.getEta(), mfttrack1.sign()); - } - } // end of loop bc_mfttrack + Produces tableMatchingParams; + Produces tableTagMatchingParams; + Produces tableProbeMatchingParams; + Produces tableMixMatchingParams; + Produces tableMuonPair; - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - /////////////// SAME EVENT /////////////// - ////////////////////////////////////////////////////////////////////////////////////////////////////////////////// + Service ccdbManager; - std::vector& mfttrackGlobalIndex = map_mfttraks[ibc]; + o2::field::MagneticField* fieldB; + o2::ccdb::CcdbApi ccdbApi; - for (int idmfttrack1 = 0; idmfttrack1 < static_cast(mfttrackGlobalIndex.size()); ++idmfttrack1) { + template + void initCCDB(BC const& bc) + { + if (mRunNumber == bc.runNumber()) + return; - auto mfttrack1 = mfttracks.rawIteratorAt(mfttrackGlobalIndex[idmfttrack1]); - o2::track::TrackParCovFwd mfttrack_at_matching = PropagateMFTtoMatchingPlane(mfttrack1, fwdtrack1); - V1.SetPtEtaPhi(mfttrack_at_matching.getPt(), mfttrack_at_matching.getEta(), mfttrack_at_matching.getPhi()); - V2.SetPtEtaPhi(fwdtrack1.pt(), fwdtrack1.eta(), fwdtrack1.phi()); + mRunNumber = bc.runNumber(); + std::map metadata; + auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(ccdbApi, mRunNumber); + auto ts = soreor.first; + auto grpmag = ccdbApi.retrieveFromTFileAny(grpmagPath, metadata, ts); + o2::base::Propagator::initFieldFromGRP(grpmag); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdbManager->get(geoPath); + } + o2::mch::TrackExtrap::setField(); + fieldB = static_cast(TGeoGlobalMagField::Instance()->GetField()); + } - double deltaPt = mfttrack_at_matching.getPt() - fwdtrack1.pt(); - double deltaX = mfttrack_at_matching.getX() - fwdtrack1.x(); - double deltaY = mfttrack_at_matching.getY() - fwdtrack1.y(); - double deltaPhi = V1.DeltaPhi(V2); - double deltaEta = mfttrack_at_matching.getEta() - fwdtrack1.eta(); + void init(o2::framework::InitContext&) + { + ccdbManager->setURL(ccdburl); + ccdbManager->setCaching(true); + ccdbManager->setLocalObjectValidityChecking(); + ccdbManager->setFatalWhenNull(false); + ccdbApi.init(ccdburl); + mRunNumber = 0; + } - if (fabs(deltaX) > fPreselectMatchingX) { - continue; - } - if (fabs(deltaY) > fPreselectMatchingY) { + void process(MyCollisions const& collisions, + MyBCs const& bcs, + MyMUONs const& muontracks, + MyMFTs const& mfttracks) + { + LOG(info) << "Process() "; + map_muontracks.clear(); + map_mfttracks.clear(); + map_collisions.clear(); + map_has_muontracks_collisions.clear(); + map_has_mfttracks_collisions.clear(); + + initCCDB(bcs.begin()); + setMUONs(muontracks, collisions); + setMFTs(mfttracks, collisions, fieldB); + + for (auto map_has_muontracks_collision : map_has_muontracks_collisions) { + auto idmuontrack_collisions = map_has_muontracks_collision.first; + for (auto map_has_mfttracks_collision : map_has_mfttracks_collisions) { + auto idmfttrack_collisions = map_has_mfttracks_collision.first; + if (idmuontrack_collisions != idmfttrack_collisions) continue; - } - - o2::track::TrackParCovFwd mfttrack_at_dca = PropagateMFT(mfttrack1, ProagationPoint::ToDCA); - - matchingParams(deltaPt, deltaEta, deltaPhi, deltaX, deltaY, - fwdtrackAtPV1.getPt(), fwdtrackAtPV1.getEta(), fwdtrack1.sign(), - mfttrack_at_dca.getEta(), mfttrack1.sign()); - - } // end of loop idmfttrack1 + map_collisions[idmfttrack_collisions] = true; + } + } - std::vector& fwdtrackGlobalIndex = map_fwdtraks[ibc]; + for (auto const& map_collision : map_collisions) { + auto const& collision = collisions.rawIteratorAt(map_collision.first); - for (int idfwdtrack2 = 0; idfwdtrack2 < static_cast(fwdtrackGlobalIndex.size()); ++idfwdtrack2) { + for (auto const& imuontrack1 : map_muontracks[map_collision.first]) { + auto const& muontrack1 = muontracks.rawIteratorAt(imuontrack1); - if (fwdtrack1.globalIndex() == fwdtrackGlobalIndex[idfwdtrack2]) { - continue; - } + for (auto const& imfttrack1 : map_mfttracks[map_collision.first]) { + auto const& mfttrack1 = mfttracks.rawIteratorAt(imfttrack1); - auto fwdtrack2 = fwdtracks.rawIteratorAt(fwdtrackGlobalIndex[idfwdtrack2]); + MatchingParamsML matching(muontrack1, mfttrack1, collision, fMatchingMethod, fieldB); + matching.calcMatchingParams(); - if (!isGoodMUONTrack(fwdtrack2)) { - continue; - } - - o2::dataformats::GlobalFwdTrack fwdtrackAtPV2 = PropagateMUONTrack(fwdtrack2, ProagationPoint::ToVtx); - if (!isGoodKenematicMUONTrack(fwdtrackAtPV2)) { - continue; - } + if (!isPassMatchingPreselection(matching.getDx(), matching.getDy())) + continue; - muon1LV.SetPtEtaPhiM(fwdtrackAtPV1.getPt(), fwdtrackAtPV1.getEta(), fwdtrackAtPV1.getPhi(), mMu); - muon2LV.SetPtEtaPhiM(fwdtrackAtPV2.getPt(), fwdtrackAtPV2.getEta(), fwdtrackAtPV2.getPhi(), mMu); - dimuonLV = muon1LV + muon2LV; + matching.calcGlobalMuonParams(); - muonPairs(nmfttracks[ibc], fwdtrack1.sign() + fwdtrack2.sign(), dimuonLV.M(), dimuonLV.Pt(), dimuonLV.Rapidity()); + bool isTrue = false; // isCorrectMatching(muontrack1,mfttrack1); - if (fabs(fwdtrack1.sign() + fwdtrack2.sign()) > 0) { - continue; - } - if (fTagMassWindowMin > dimuonLV.M() || dimuonLV.M() > fTagMassWindowMax) { - continue; - } - if (nmfttracks[ibc] < 1) { - continue; + tableMatchingParams(matching.getGMPtAtDCA(), + matching.getGMEtaAtDCA(), + static_cast(matching.getGMQ()), + matching.getDpt(), + matching.getDx(), + matching.getDy(), + matching.getDeta(), + matching.getDphi(), + isTrue); } - tagmuonPairs(nmfttracks[ibc], fwdtrack1.sign() + fwdtrack2.sign(), dimuonLV.M(), dimuonLV.Pt(), dimuonLV.Rapidity()); - - bool isGoodTag = false; - int nMFTCandsTagMUON = 0; - - auto tagfwdtrack = fwdtracks.rawIteratorAt(selectTagMUON(fwdtrack1, fwdtrack2)); - o2::dataformats::GlobalFwdTrack tagfwdtrackAtPV = PropagateMUONTrack(tagfwdtrack, ProagationPoint::ToVtx); - - for (int idmfttrack1 = 0; idmfttrack1 < static_cast(mfttrackGlobalIndex.size()); ++idmfttrack1) { - - auto mfttrack1 = mfttracks.rawIteratorAt(mfttrackGlobalIndex[idmfttrack1]); - o2::track::TrackParCovFwd mfttrack_at_matching = PropagateMFTtoMatchingPlane(mfttrack1, tagfwdtrack); - - V1.SetPtEtaPhi(mfttrack_at_matching.getPt(), mfttrack_at_matching.getEta(), mfttrack_at_matching.getPhi()); - V2.SetPtEtaPhi(tagfwdtrack.pt(), tagfwdtrack.eta(), tagfwdtrack.phi()); - - double deltaPt = mfttrack_at_matching.getPt() - tagfwdtrack.pt(); - double deltaX = mfttrack_at_matching.getX() - tagfwdtrack.x(); - double deltaY = mfttrack_at_matching.getY() - tagfwdtrack.y(); - double deltaPhi = V1.DeltaPhi(V2); - double deltaEta = mfttrack_at_matching.getEta() - tagfwdtrack.eta(); - - if (fabs(deltaX) > fPreselectMatchingX) { + for (auto const& map_mfttrack : map_mfttracks) { + if (map_mfttrack.first == map_collision.first) continue; - } - if (fabs(deltaY) > fPreselectMatchingY) { + if (fabs(map_vtxz[map_mfttrack.first] - map_vtxz[map_collision.first]) > fEventMaxDeltaVtxZ) + continue; + if (fabs(map_nmfttrack[map_mfttrack.first] - map_nmfttrack[map_collision.first]) > fEventMaxDeltaNMFT) continue; - } - - o2::track::TrackParCovFwd mfttrack_at_dca = PropagateMFT(mfttrack1, ProagationPoint::ToDCA); - - bool dummyTag = false; - if (fTagXWindowLow < deltaX && deltaX < fTagXWindowUp && - fTagXWindowLow < deltaY && deltaY < fTagXWindowUp && - fTagPhiWindowLow < deltaPhi && deltaPhi < fTagPhiWindowUp && - fTagEtaWindowLow < deltaEta && deltaEta < fTagEtaWindowUp) { - isGoodTag = true; - dummyTag = true; - nMFTCandsTagMUON++; + for (auto const& imfttrack1 : map_mfttrack.second) { + auto const& mfttrack1 = mfttracks.rawIteratorAt(imfttrack1); + MatchingParamsML matching(muontrack1, mfttrack1, collision, fMatchingMethod, fieldB); + matching.calcMatchingParams(); + if (!isPassMatchingPreselection(matching.getDx(), matching.getDy())) + continue; + matching.calcGlobalMuonParams(); + + bool isTrue = false; // isCorrectMatching(muontrack1,mfttrack1); + + tableMixMatchingParams(matching.getGMPtAtDCA(), + matching.getGMEtaAtDCA(), + static_cast(matching.getGMQ()), + matching.getDpt(), + matching.getDx(), + matching.getDy(), + matching.getDeta(), + matching.getDphi(), + isTrue); } - - tagmatchingParams(deltaPt, deltaEta, deltaPhi, deltaX, deltaY, - tagfwdtrackAtPV.getPt(), tagfwdtrackAtPV.getEta(), tagfwdtrack.sign(), mfttrack_at_dca.getEta(), mfttrack1.sign(), dummyTag); - } - - if (!isGoodTag) { - continue; } - auto probefwdtrack = fwdtracks.rawIteratorAt(selectProbeMUONTrack(fwdtrack1, fwdtrack2)); - o2::dataformats::GlobalFwdTrack probefwdtrackAtPV = PropagateMUONTrack(probefwdtrack, ProagationPoint::ToVtx); - - int nMFTCandsProbeMUON = 0; - /// Counting the number of MFT tracks in a probe muon track - for (int idmfttrack1 = 0; idmfttrack1 < static_cast(mfttrackGlobalIndex.size()); ++idmfttrack1) { - - auto mfttrack1 = mfttracks.rawIteratorAt(mfttrackGlobalIndex[idmfttrack1]); - o2::track::TrackParCovFwd mfttrack_at_matching = PropagateMFTtoMatchingPlane(mfttrack1, probefwdtrack); - V1.SetPtEtaPhi(mfttrack_at_matching.getPt(), mfttrack_at_matching.getEta(), mfttrack_at_matching.getPhi()); - V2.SetPtEtaPhi(probefwdtrack.pt(), probefwdtrack.eta(), probefwdtrack.phi()); - - double deltaX = mfttrack_at_matching.getX() - probefwdtrack.x(); - double deltaY = mfttrack_at_matching.getY() - probefwdtrack.y(); - double deltaPhi = V1.DeltaPhi(V2); - double deltaEta = mfttrack_at_matching.getEta() - probefwdtrack.eta(); + for (auto const& imuontrack2 : map_muontracks[map_collision.first]) { - if (fProbeXWindowLow < deltaX && deltaX < fProbeXWindowUp && - fProbeXWindowLow < deltaY && deltaY < fProbeXWindowUp && - fProbePhiWindowLow < deltaPhi && deltaPhi < fProbePhiWindowUp && - fProbeEtaWindowLow < deltaEta && deltaEta < fProbeEtaWindowUp) { - nMFTCandsProbeMUON++; - } - } + if (imuontrack1 == imuontrack2) + continue; - for (int idmfttrack1 = 0; idmfttrack1 < static_cast(mfttrackGlobalIndex.size()); ++idmfttrack1) { + auto const& muontrack2 = muontracks.rawIteratorAt(imuontrack2); - auto mfttrack1 = mfttracks.rawIteratorAt(mfttrackGlobalIndex[idmfttrack1]); - o2::track::TrackParCovFwd mfttrack_at_matching = PropagateMFTtoMatchingPlane(mfttrack1, probefwdtrack); + FindTagAndProbe tagdimuon(muontrack1, muontrack2, collision); + tagdimuon.calcMuonPairAtPV(); + tableMuonPair(tagdimuon.getCharge(), tagdimuon.getMass(), tagdimuon.getPt(), tagdimuon.getRap()); - V1.SetPtEtaPhi(mfttrack_at_matching.getPt(), mfttrack_at_matching.getEta(), mfttrack_at_matching.getPhi()); - V2.SetPtEtaPhi(probefwdtrack.pt(), probefwdtrack.eta(), probefwdtrack.phi()); + if (!isGoodTagDimuon(tagdimuon.getMass())) + continue; - double deltaPt = mfttrack_at_matching.getPt() - probefwdtrack.pt(); - double deltaX = mfttrack_at_matching.getX() - probefwdtrack.x(); - double deltaY = mfttrack_at_matching.getY() - probefwdtrack.y(); - double deltaPhi = V1.DeltaPhi(V2); - double deltaEta = mfttrack_at_matching.getEta() - probefwdtrack.eta(); + auto tagmuontrack = muontrack1; + auto probemuontrack = muontrack2; - if (fabs(deltaX) > fPreselectMatchingX) { - continue; - } - if (fabs(deltaY) > fPreselectMatchingY) { - continue; + if (tagdimuon.getTagMuonIndex() == 1) { + tagmuontrack = muontrack2; + probemuontrack = muontrack1; } - o2::track::TrackParCovFwd mfttrack_at_dca = PropagateMFT(mfttrack1, ProagationPoint::ToDCA); - - probematchingParams(nMFTCandsTagMUON, nMFTCandsProbeMUON, tagfwdtrack.p(), deltaPt, deltaEta, deltaPhi, deltaX, deltaY, probefwdtrackAtPV.getPt(), probefwdtrackAtPV.getEta(), probefwdtrack.sign(), mfttrack_at_dca.getEta(), mfttrack1.sign()); - } - } - } - } + int nTagMFTCand = 0; + int nProbeMFTCand = 0; + + unordered_map> map_tagMatchingParams; + unordered_map> map_probeMatchingParams; + + for (auto const& imfttrack1 : map_mfttracks[map_collision.first]) { + auto const& mfttrack1 = mfttracks.rawIteratorAt(imfttrack1); + MatchingParamsML matchingTag(tagmuontrack, mfttrack1, collision, fMatchingMethod, fieldB); + matchingTag.calcMatchingParams(); + matchingTag.calcGlobalMuonParams(); + if (isGoodTagMatching(matchingTag.getDx(), matchingTag.getDy(), matchingTag.getDeta(), matchingTag.getDphi()) && + isPassMatchingPreselection(matchingTag.getDx(), matchingTag.getDy())) { + bool isTrue = false; // isCorrectMatching(tagmuontrack,mfttrack1); + tableTagMatchingParams(matchingTag.getGMPtAtDCA(), + matchingTag.getGMEtaAtDCA(), + matchingTag.getGMQ(), + matchingTag.getDpt(), + matchingTag.getDx(), + matchingTag.getDy(), + matchingTag.getDeta(), + matchingTag.getDphi(), + isTrue); + ++nTagMFTCand; + } + + if (nTagMFTCand == 1) { + MatchingParamsML matchingProbe(probemuontrack, mfttrack1, collision, fMatchingMethod, fieldB); + matchingProbe.calcMatchingParams(); + matchingProbe.calcGlobalMuonParams(); + if (isPassMatchingPreselection(matchingProbe.getDx(), matchingProbe.getDy())) { + bool isTrue = false; // isCorrectMatching(probemuontrack,mfttrack1); + tableProbeMatchingParams(matchingTag.getGMPtAtDCA(), + matchingProbe.getGMPtAtDCA(), + matchingProbe.getGMEtaAtDCA(), + matchingProbe.getGMQ(), + matchingProbe.getDpt(), + matchingProbe.getDx(), + matchingProbe.getDy(), + matchingProbe.getDeta(), + matchingProbe.getDphi(), + isTrue); + } + } // nTagMFTCand + } // end of loop imfttrack1 + } // end of loop imuontrack2 + } // end of loop imuontrack1 + } // end of loop map_collision + + } // end of processMC }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; }