diff --git a/RecoMTD/TimingIDTools/interface/MTDTrackQualityMVA.h b/RecoMTD/TimingIDTools/interface/MTDTrackQualityMVA.h index 4d31e1acdeff0..38f3329ff79a5 100644 --- a/RecoMTD/TimingIDTools/interface/MTDTrackQualityMVA.h +++ b/RecoMTD/TimingIDTools/interface/MTDTrackQualityMVA.h @@ -10,20 +10,23 @@ #include "CommonTools/MVAUtils/interface/TMVAEvaluator.h" #define MTDTRACKQUALITYMVA_VARS(MTDBDTVAR) \ - MTDBDTVAR(pt) \ - MTDBDTVAR(eta) \ - MTDBDTVAR(phi) \ - MTDBDTVAR(chi2) \ - MTDBDTVAR(ndof) \ - MTDBDTVAR(numberOfValidHits) \ - MTDBDTVAR(numberOfValidPixelBarrelHits) \ - MTDBDTVAR(numberOfValidPixelEndcapHits) \ - MTDBDTVAR(btlMatchChi2) \ - MTDBDTVAR(btlMatchTimeChi2) \ - MTDBDTVAR(etlMatchChi2) \ - MTDBDTVAR(etlMatchTimeChi2) \ - MTDBDTVAR(mtdt) \ - MTDBDTVAR(path_len) + MTDBDTVAR(Track_pt) \ + MTDBDTVAR(Track_eta) \ + MTDBDTVAR(Track_phi) \ + MTDBDTVAR(Track_dz) \ + MTDBDTVAR(Track_dxy) \ + MTDBDTVAR(Track_chi2) \ + MTDBDTVAR(Track_ndof) \ + MTDBDTVAR(Track_npixBarrelValidHits) \ + MTDBDTVAR(Track_npixEndcapValidHits) \ + MTDBDTVAR(Track_BTLchi2) \ + MTDBDTVAR(Track_BTLtime_chi2) \ + MTDBDTVAR(Track_ETLchi2) \ + MTDBDTVAR(Track_ETLtime_chi2) \ + MTDBDTVAR(Track_Tmtd) \ + MTDBDTVAR(Track_sigmaTmtd) \ + MTDBDTVAR(Track_length) \ + MTDBDTVAR(Track_lHitPos) #define MTDBDTVAR_ENUM(ENUM) ENUM, #define MTDBDTVAR_STRING(STRING) #STRING, @@ -38,6 +41,7 @@ class MTDTrackQualityMVA { //---getters--- // 4D float operator()(const reco::TrackRef& trk, + const reco::BeamSpot& beamspot, const edm::ValueMap& npixBarrels, const edm::ValueMap& npixEndcaps, const edm::ValueMap& btl_chi2s, @@ -45,7 +49,9 @@ class MTDTrackQualityMVA { const edm::ValueMap& etl_chi2s, const edm::ValueMap& etl_time_chi2s, const edm::ValueMap& tmtds, - const edm::ValueMap& trk_lengths) const; + const edm::ValueMap& sigmatmtds, + const edm::ValueMap& trk_lengths, + const edm::ValueMap& trk_lhitpos) const; private: std::vector vars_, spec_vars_; diff --git a/RecoMTD/TimingIDTools/plugins/MTDTrackQualityMVAProducer.cc b/RecoMTD/TimingIDTools/plugins/MTDTrackQualityMVAProducer.cc index 54fd444fe946d..cfd0bf1f218ea 100644 --- a/RecoMTD/TimingIDTools/plugins/MTDTrackQualityMVAProducer.cc +++ b/RecoMTD/TimingIDTools/plugins/MTDTrackQualityMVAProducer.cc @@ -37,15 +37,17 @@ class MTDTrackQualityMVAProducer : public edm::stream::EDProducer<> { edm::EDGetTokenT tracksToken_; edm::EDGetTokenT tracksMTDToken_; - edm::EDGetTokenT> btlMatchChi2Token_; + edm::EDGetTokenT RecBeamSpotToken_; edm::EDGetTokenT> btlMatchTimeChi2Token_; edm::EDGetTokenT> etlMatchChi2Token_; edm::EDGetTokenT> etlMatchTimeChi2Token_; edm::EDGetTokenT> mtdTimeToken_; + edm::EDGetTokenT> sigmamtdTimeToken_; edm::EDGetTokenT> pathLengthToken_; edm::EDGetTokenT> npixBarrelToken_; edm::EDGetTokenT> npixEndcapToken_; + edm::EDGetTokenT> outermostHitPositionToken_; MTDTrackQualityMVA mva_; }; @@ -53,15 +55,19 @@ class MTDTrackQualityMVAProducer : public edm::stream::EDProducer<> { MTDTrackQualityMVAProducer::MTDTrackQualityMVAProducer(const ParameterSet& iConfig) : tracksToken_(consumes(iConfig.getParameter("tracksSrc"))), btlMatchChi2Token_(consumes>(iConfig.getParameter("btlMatchChi2Src"))), + RecBeamSpotToken_(consumes(iConfig.getParameter("offlineBS"))), btlMatchTimeChi2Token_( consumes>(iConfig.getParameter("btlMatchTimeChi2Src"))), etlMatchChi2Token_(consumes>(iConfig.getParameter("etlMatchChi2Src"))), etlMatchTimeChi2Token_( consumes>(iConfig.getParameter("etlMatchTimeChi2Src"))), mtdTimeToken_(consumes>(iConfig.getParameter("mtdTimeSrc"))), + sigmamtdTimeToken_(consumes>(iConfig.getParameter("sigmamtdTimeSrc"))), pathLengthToken_(consumes>(iConfig.getParameter("pathLengthSrc"))), npixBarrelToken_(consumes>(iConfig.getParameter("npixBarrelSrc"))), npixEndcapToken_(consumes>(iConfig.getParameter("npixEndcapSrc"))), + outermostHitPositionToken_( + consumes>(iConfig.getParameter("outermostHitPositionSrc"))), mva_(iConfig.getParameter("qualityBDT_weights_file").fullPath()) { produces>(mvaName); } @@ -79,15 +85,20 @@ void MTDTrackQualityMVAProducer::fillDescriptions(edm::ConfigurationDescriptions desc.add("etlMatchTimeChi2Src", edm::InputTag("trackExtenderWithMTD", "etlMatchTimeChi2")) ->setComment("ETL Chi2 Matching value Map"); desc.add("mtdTimeSrc", edm::InputTag("trackExtenderWithMTD", "generalTracktmtd")) - ->setComment("MTD TIme value Map"); + ->setComment("MTD Time value Map"); + desc.add("sigmamtdTimeSrc", edm::InputTag("trackExtenderWithMTD", "generalTracksigmatmtd")) + ->setComment("sigma MTD Time value Map"); desc.add("pathLengthSrc", edm::InputTag("trackExtenderWithMTD", "generalTrackPathLength")) ->setComment("MTD PathLength value Map"); desc.add("npixBarrelSrc", edm::InputTag("trackExtenderWithMTD", "npixBarrel")) ->setComment("# of Barrel pixel associated to refitted tracks"); desc.add("npixEndcapSrc", edm::InputTag("trackExtenderWithMTD", "npixEndcap")) ->setComment("# of Endcap pixel associated to refitted tracks"); + desc.add("outermostHitPositionSrc", + edm::InputTag("trackExtenderWithMTD", "generalTrackOutermostHitPosition")); + desc.add("offlineBS", edm::InputTag("offlineBeamSpot")); desc.add("qualityBDT_weights_file", - edm::FileInPath("RecoMTD/TimingIDTools/data/clf4D_MTDquality_bo.xml")) + edm::FileInPath("RecoMTD/TimingIDTools/data/BDT_nvars_17_d7.xml")) ->setComment("Track MTD quality BDT weights"); descriptions.add("mtdTrackQualityMVAProducer", desc); } @@ -109,6 +120,11 @@ void MTDTrackQualityMVAProducer::produce(edm::Event& ev, const edm::EventSetup& ev.getByToken(tracksToken_, tracksH); const auto& tracks = *tracksH; + reco::BeamSpot beamSpot; + edm::Handle BeamSpotH; + ev.getByToken(RecBeamSpotToken_, BeamSpotH); + beamSpot = *BeamSpotH; + const auto& btlMatchChi2 = ev.get(btlMatchChi2Token_); const auto& btlMatchTimeChi2 = ev.get(btlMatchTimeChi2Token_); const auto& etlMatchChi2 = ev.get(etlMatchChi2Token_); @@ -117,6 +133,8 @@ void MTDTrackQualityMVAProducer::produce(edm::Event& ev, const edm::EventSetup& const auto& npixBarrel = ev.get(npixBarrelToken_); const auto& npixEndcap = ev.get(npixEndcapToken_); const auto& mtdTime = ev.get(mtdTimeToken_); + const auto& sigmamtdTime = ev.get(sigmamtdTimeToken_); + const auto& lHitPos = ev.get(outermostHitPositionToken_); std::vector mvaOutRaw; @@ -127,6 +145,7 @@ void MTDTrackQualityMVAProducer::produce(edm::Event& ev, const edm::EventSetup& mvaOutRaw.push_back(-1.); else { mvaOutRaw.push_back(mva_(trackref, + beamSpot, npixBarrel, npixEndcap, btlMatchChi2, @@ -134,7 +153,9 @@ void MTDTrackQualityMVAProducer::produce(edm::Event& ev, const edm::EventSetup& etlMatchChi2, etlMatchTimeChi2, mtdTime, - pathLength)); + sigmamtdTime, + pathLength, + lHitPos)); } } fillValueMap(ev, tracksH, mvaOutRaw, mvaName); diff --git a/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc b/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc index 0ee37dde7ca88..ccfbde216417b 100644 --- a/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc +++ b/RecoMTD/TimingIDTools/plugins/MVATrainingNtuple.cc @@ -11,6 +11,8 @@ #include "DataFormats/VertexReco/interface/Vertex.h" #include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h" +#include "DataFormats/ForwardDetId/interface/ETLDetId.h" +// #include "DataFormats/ForwardDetId/interface/BTLDetId.h" #include "SimDataFormats/Associations/interface/TrackToTrackingParticleAssociator.h" #include "SimDataFormats/Associations/interface/MtdSimLayerClusterToTPAssociatorBaseImpl.h" @@ -46,11 +48,12 @@ class MVATrainingNtuple : public edm::one::EDAnalyzer private: void analyze(const edm::Event&, const edm::EventSetup&) override; + const bool trkTPSelAll(const TrackingParticle&); + const bool trkRecSel(const reco::TrackBase&); + const edm::Ref>* getAnyMatchedTP(const reco::TrackBaseRef&); double timeFromTrueMass(double, double, double, double); - bool isSameCluster(const FTLCluster&, const FTLCluster&); - std::vector getSimPVs(const edm::Handle&); edm::Service fs_; @@ -61,18 +64,19 @@ class MVATrainingNtuple : public edm::one::EDAnalyzer static constexpr double simUnit_ = 1e9; //sim time in s while reco time in ns static constexpr double c_ = 2.99792458e1; //c in cm/ns - static constexpr double BTL_eta_cut = 1.5; + static constexpr double etacutGEN_ = 4.; // |eta| < 4; + static constexpr double etacutREC_ = 3.; // |eta| < 3; + static constexpr double etacutBTL_ = 1.5; // |eta| < 1.5; + static constexpr double pTcutBTL_ = 0.7; // PT > 0.7 GeV + static constexpr double pTcutETL_ = 0.2; // PT > 0.2 GeV + static constexpr double rBTL_ = 110.0; + static constexpr double zETL_ = 290.0; std::string fileName_; bool saveNtupleforBDT_; bool saveNtupleforGNN_; - // cuts for BDT training input - static constexpr double BDT_track_eta_cut = 3.0; - static constexpr double BDT_track_pt_cut = 0.5; - const reco::RecoToSimCollection* r2s_; - const reco::SimToRecoCollection* s2r_; // GNN input variables std::vector gnn_pt, gnn_eta, gnn_phi, gnn_z_pca, gnn_dz, gnn_t_Pi, gnn_t_K, gnn_t_P, gnn_t0safe, gnn_t0pid, @@ -87,10 +91,10 @@ class MVATrainingNtuple : public edm::one::EDAnalyzer // BDT input variables std::vector Ttrack_pt, Ttrack_eta, Ttrack_phi, Ttrack_dz, Ttrack_dxy, Ttrack_chi2, Ttrack_BTLchi2, Ttrack_BTLtime_chi2, Ttrack_ETLchi2, Ttrack_ETLtime_chi2, Ttrack_t0, Ttrack_sigmat0, Ttrack_Tmtd, - Ttrack_sigmaTmtd, Ttrack_lenght, Ttrack_MtdMVA, Ttrack_lHitPos, TtrackTP_pt, TtrackTP_eta, TtrackTP_phi; - std::vector Ttrack_ndof, Ttrack_nValidHits, Ttrack_npixBarrelValidHits, Ttrack_npixEndcapValidHits, - TtrackTP_nValidHits; - std::vector Ttrack_Signal, Ttrack_Associated; + Ttrack_sigmaTmtd, Ttrack_length, Ttrack_MtdMVA, Ttrack_lHitPos; + std::vector Ttrack_ndof, Ttrack_nValidHits, Ttrack_npixBarrelValidHits, Ttrack_npixEndcapValidHits; + std::vector Ttrack_Signal, Ttrack_Associated, Ttrack_TPDirClu, Ttrack_TPOtherClu, Ttrack_isBTL, Ttrack_isETL, + Ttrack_withMTD; edm::EDGetTokenT> btlMatchChi2Token_; edm::EDGetTokenT> btlMatchTimeChi2Token_; @@ -118,8 +122,9 @@ class MVATrainingNtuple : public edm::one::EDAnalyzer edm::EDGetTokenT> momentumToken_; edm::EDGetTokenT> sigmatimeToken_; edm::EDGetTokenT> t0SrcToken_; - edm::EDGetTokenT> Sigmat0SrcToken_; + edm::EDGetTokenT> sigmat0SrcToken_; edm::EDGetTokenT> t0PidToken_; + edm::EDGetTokenT> sigmat0PidToken_; edm::EDGetTokenT> t0SafePidToken_; edm::EDGetTokenT> sigmat0SafePidToken_; edm::EDGetTokenT> trackMVAQualToken_; @@ -160,8 +165,9 @@ MVATrainingNtuple::MVATrainingNtuple(const edm::ParameterSet& iConfig) momentumToken_ = consumes>(iConfig.getParameter("momentumSrc")); sigmatimeToken_ = consumes>(iConfig.getParameter("sigmaSrc")); t0SrcToken_ = consumes>(iConfig.getParameter("t0Src")); - Sigmat0SrcToken_ = consumes>(iConfig.getParameter("sigmat0Src")); + sigmat0SrcToken_ = consumes>(iConfig.getParameter("sigmat0Src")); t0PidToken_ = consumes>(iConfig.getParameter("t0PID")); + sigmat0PidToken_ = consumes>(iConfig.getParameter("sigmat0PID")); t0SafePidToken_ = consumes>(iConfig.getParameter("t0SafePID")); sigmat0SafePidToken_ = consumes>(iConfig.getParameter("sigmat0SafePID")); trackMVAQualToken_ = consumes>(iConfig.getParameter("trackMVAQual")); @@ -230,11 +236,6 @@ double MVATrainingNtuple::timeFromTrueMass(double mass, double pathlength, doubl } } -bool MVATrainingNtuple::isSameCluster(const FTLCluster& clu1, const FTLCluster& clu2) { - return clu1.id() == clu2.id() && clu1.size() == clu2.size() && clu1.x() == clu2.x() && clu1.y() == clu2.y() && - clu1.time() == clu2.time(); -} - std::vector MVATrainingNtuple::getSimPVs( const edm::Handle& tVC) { std::vector simpv; @@ -331,14 +332,12 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& const auto& trackAssoc = iEvent.get(trackAssocToken_); - std::vector vertices; - edm::Handle> RecVertexHandle; - iEvent.getByToken(RecVertexToken_, RecVertexHandle); - vertices = *RecVertexHandle; - const auto& tp2SimAssociationMap = iEvent.get(tp2SimAssociationMapToken_); const auto& r2sAssociationMap = iEvent.get(r2sAssociationMapToken_); + // auto RecTrackHandle = makeValid(iEvent.getHandle(RecTrackToken_)); + auto RecTrackHandle = iEvent.getHandle(RecTrackToken_); + const auto& btlRecCluHandle = iEvent.getHandle(btlRecCluToken_); const auto& etlRecCluHandle = iEvent.getHandle(etlRecCluToken_); @@ -346,8 +345,9 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& const auto& momentum = iEvent.get(momentumToken_); const auto& sigmatimemtd = iEvent.get(sigmatimeToken_); const auto& t0Src = iEvent.get(t0SrcToken_); - const auto& Sigmat0Src = iEvent.get(Sigmat0SrcToken_); + const auto& sigmat0Src = iEvent.get(sigmat0SrcToken_); const auto& t0Pid = iEvent.get(t0PidToken_); + const auto& sigmat0Pid = iEvent.get(sigmat0PidToken_); const auto& t0Safe = iEvent.get(t0SafePidToken_); const auto& sigmat0Safe = iEvent.get(sigmat0SafePidToken_); const auto& mtdQualMVA = iEvent.get(trackMVAQualToken_); @@ -560,10 +560,6 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& BDTtree->Branch("Track_chi2", &Ttrack_chi2); BDTtree->Branch("Track_ndof", &Ttrack_ndof); BDTtree->Branch("Track_nValidHits", &Ttrack_nValidHits); - BDTtree->Branch("TrackTP_pt", &TtrackTP_pt); - BDTtree->Branch("TrackTP_eta", &TtrackTP_eta); - BDTtree->Branch("TrackTP_phi", &TtrackTP_phi); - BDTtree->Branch("TrackTP_nValidHits", &TtrackTP_nValidHits); BDTtree->Branch("Track_npixBarrelValidHits", &Ttrack_npixBarrelValidHits); BDTtree->Branch("Track_npixEndcapValidHits", &Ttrack_npixEndcapValidHits); BDTtree->Branch("Track_Signal", &Ttrack_Signal); @@ -578,7 +574,12 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& BDTtree->Branch("Track_MtdMVA", &Ttrack_MtdMVA); BDTtree->Branch("Track_lHitPos", &Ttrack_lHitPos); BDTtree->Branch("Track_sigmaTmtd", &Ttrack_sigmaTmtd); - BDTtree->Branch("Track_lenght", &Ttrack_lenght); + BDTtree->Branch("Track_length", &Ttrack_length); + BDTtree->Branch("Track_DirClu", &Ttrack_TPDirClu); + BDTtree->Branch("Track_OtherClu", &Ttrack_TPOtherClu); + BDTtree->Branch("Track_isBTL", &Ttrack_isBTL); + BDTtree->Branch("Track_isETL", &Ttrack_isETL); + BDTtree->Branch("Track_withMTD", &Ttrack_withMTD); Ttrack_pt.clear(); Ttrack_eta.clear(); @@ -588,10 +589,6 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& Ttrack_chi2.clear(); Ttrack_ndof.clear(); Ttrack_nValidHits.clear(); - TtrackTP_pt.clear(); - TtrackTP_eta.clear(); - TtrackTP_phi.clear(); - TtrackTP_nValidHits.clear(); Ttrack_npixBarrelValidHits.clear(); Ttrack_npixEndcapValidHits.clear(); Ttrack_Signal.clear(); @@ -606,10 +603,17 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& Ttrack_MtdMVA.clear(); Ttrack_lHitPos.clear(); Ttrack_sigmaTmtd.clear(); - Ttrack_lenght.clear(); + Ttrack_length.clear(); + Ttrack_TPDirClu.clear(); + Ttrack_TPOtherClu.clear(); + Ttrack_isBTL.clear(); + Ttrack_isETL.clear(); + Ttrack_withMTD.clear(); unsigned int index = 0; - for (const auto& trackGen : *tracksH) { + + // --- Loop over all RECO tracks --- + for (const auto& trackGen : *RecTrackHandle) { const reco::TrackRef trackref(iEvent.getHandle(RecTrackToken_), index); index++; @@ -621,160 +625,185 @@ void MVATrainingNtuple::analyze(const edm::Event& iEvent, const edm::EventSetup& const reco::TrackRef mtdTrackref = reco::TrackRef(iEvent.getHandle(RecMTDTrackToken_), trackAssoc[trackref]); const reco::Track& track = *mtdTrackref; - if (std::abs(trackGen.eta()) < BDT_track_eta_cut && trackGen.pt() > BDT_track_pt_cut) { - const reco::TrackBaseRef trkrefb(trackref); - auto found = r2s_->find(trkrefb); // Find TP! - if (found != r2s_->end()) { - bool good_association = false; - - Ttrack_pt.push_back(trackGen.pt()); - Ttrack_phi.push_back(trackGen.phi()); - Ttrack_eta.push_back(trackGen.eta()); - Ttrack_dz.push_back(std::abs(trackGen.dz(beamSpot.position()))); - Ttrack_dxy.push_back(std::abs(trackGen.dxy(beamSpot.position()))); - Ttrack_chi2.push_back(trackGen.chi2()); - Ttrack_ndof.push_back(trackGen.ndof()); - Ttrack_nValidHits.push_back(trackGen.numberOfValidHits()); - - Ttrack_npixBarrelValidHits.push_back(npixBarrel[trackref]); - Ttrack_npixEndcapValidHits.push_back(npixEndcap[trackref]); - Ttrack_BTLchi2.push_back(btlMatchChi2[trackref]); - Ttrack_BTLtime_chi2.push_back(btlMatchTimeChi2[trackref]); - Ttrack_ETLchi2.push_back(etlMatchChi2[trackref]); - Ttrack_ETLtime_chi2.push_back(etlMatchTimeChi2[trackref]); - - Ttrack_t0.push_back(t0Src[trackref]); - Ttrack_sigmat0.push_back(Sigmat0Src[trackref]); - Ttrack_Tmtd.push_back(tMtd[trackref]); - Ttrack_sigmaTmtd.push_back(sigmatimemtd[trackref]); - Ttrack_lenght.push_back(pathLength[trackref]); - Ttrack_MtdMVA.push_back(mtdQualMVA[trackref]); - Ttrack_lHitPos.push_back(outermostHitPosition[trackref]); - - const auto& tp = (found->val) - [0]; // almost all tracks have just one TP, a few have 2. (can scan through with "for(const auto& tp : found->val)") - - TtrackTP_pt.push_back(tp.first->pt()); - TtrackTP_eta.push_back(tp.first->eta()); - TtrackTP_phi.push_back(tp.first->phi()); - TtrackTP_nValidHits.push_back(tp.first->numberOfHits()); - - auto simClustersRefs = tp2SimAssociationMap.find(tp.first); // finds a simClusterReference!! - const bool withMTD = (simClustersRefs != tp2SimAssociationMap.end()); - - // 1) Link track RecHit to MTdTrackingRecHit (I know which RecHits, hit MTD) - // 2) Get the MTD Reco Cluster from MTDTrackingRecHit info - // 3) Find the MTD sim cluster that is linked to MTD reco cluster in the previous step - // 4) Check if the MTD sim cluster found in previous step is the same as MTD Sim cluster that is linked to TP. - - if (withMTD) { // TP link to MTDsimCluster - - // In test file, all TPs had only 1 simCluster linked to them - const auto& SimCluRefs = (simClustersRefs->val)[0]; - if ((*SimCluRefs).trackIdOffset() == 0) { // SimCluster linked to TP is from DirectHit!!! - - for (const auto& hit : track.recHits()) { // Extended track with MTD - if (good_association) - continue; // if goodd assoc found, do not go through all the following checks. - if (hit->isValid() == false) - continue; - - MTDDetId Hit = hit->geographicalId(); - - if ((Hit.det() == 6) && (Hit.subdetId() == 1) && - (Hit.mtdSubDetector() == 1 || Hit.mtdSubDetector() == 2)) { // trackingRecHit is a hit in MTD - - const MTDTrackingRecHit* mtdhit1 = static_cast( - hit); // Why I can't I access the mtdcluster info directly from TrackingRecHit? - const FTLCluster& hit_cluster_check = mtdhit1->mtdCluster(); - - if (abs(track.eta()) < BTL_eta_cut) { // Should be a BTL cluster - for (const auto& DetSetCluBTL : *btlRecCluHandle) { // BTL check - if (good_association) - break; - for (const auto& clusterBTL : DetSetCluBTL) { // Scan throguh btl reco clusters to find a match - if (good_association) - break; - if (isSameCluster(hit_cluster_check, - clusterBTL)) { // find the reco Cluster inside the recoCluster collections - - edm::Ref, FTLCluster> clusterRefBTL = edmNew::makeRefTo( - btlRecCluHandle, - &clusterBTL); // get the reference to reco cluster inside the collections - auto itp = r2sAssociationMap.equal_range(clusterRefBTL); // find the linked simCluster - if (itp.first != itp.second) { // find the linked simCluster - std::vector simClustersRefs_RecoMatchBTL = - (*itp.first).second; // the range of itp.first, itp.second should be always 1 - - for (unsigned int i = 0; i < simClustersRefs_RecoMatchBTL.size(); i++) { - auto simClusterRef_RecoMatchBTL = simClustersRefs_RecoMatchBTL[i]; - - if ((*simClusterRef_RecoMatchBTL).simLCTime() == - (*SimCluRefs) - .simLCTime()) { // check if the sim cluster linked to reco cluster is the same as the one linked to TP. - good_association = true; - break; - } - } - } - } else { - continue; - } // mtd hit matched to btl reco cluster - } // loop through BTL reco clusters - } // loop thorugh set of BTL reco clusters - } else { // Should be an ETL cluster - for (const auto& DetSetCluETL : *etlRecCluHandle) { // ETL check - if (good_association) - break; - for (const auto& clusterETL : DetSetCluETL) { // Scan throguh etl reco clusters to find a match - if (good_association) - break; - if (isSameCluster(hit_cluster_check, clusterETL)) { - edm::Ref, FTLCluster> clusterRefETL = - edmNew::makeRefTo(etlRecCluHandle, &clusterETL); - auto itp = r2sAssociationMap.equal_range(clusterRefETL); - if (itp.first != itp.second) { - std::vector simClustersRefs_RecoMatchETL = - (*itp.first).second; // the range of itp.first, itp.second should be always 1 - - for (unsigned int i = 0; i < simClustersRefs_RecoMatchETL.size(); i++) { - auto simClusterRef_RecoMatchETL = simClustersRefs_RecoMatchETL[i]; - - if ((*simClusterRef_RecoMatchETL).simLCTime() == (*SimCluRefs).simLCTime()) { - good_association = true; - break; - } - } - } - } else { - continue; - } // mtd hit matched to etl reco cluster - } // loop through ETL reco clusters - } // loop thorugh set of ETL reco clusters - } // BTL/ETL cluster search split - - } else { // trackingRecHit is a hit in MTD - continue; - } // Hits in MTD - } // Loop through trackHits - } - } // TP link to MTDsimCluster + bool withMTD = false; + bool isBTL = false; + bool isETL = false; + bool twoETLdiscs = false; + bool isTPDirClu = false; + bool isTPOtherClu = false; + + bool isTPmtdDirectBTL = false, isTPmtdOtherBTL = false, isTPmtdDirectETL = false, isTPmtdOtherETL = false, + isTPmtdCorrect = false; + + if (trkRecSel(trackGen)) { + if (std::round(sigmatimemtd[trackref] - sigmat0Pid[trackref]) != 0) { + LogWarning("mtdTracks") + << "TimeError associated to refitted track is different from TimeError stored in tofPID " + "sigmat0 ValueMap: this should not happen"; + } - if (tp.first->eventId().bunchCrossing() == 0 && - tp.first->eventId().event() == 0) { // Signal vs PU seperation - Ttrack_Signal.push_back(true); // Signal track - } else { - Ttrack_Signal.push_back(false); // PU track? + std::vector, FTLCluster>> recoClustersRefs; + + if (std::abs(trackGen.eta()) < etacutBTL_) { + // --- all BTL tracks (with and without hit in MTD) --- + + for (const auto hit : track.recHits()) { + if (hit->isValid() == false) + continue; + MTDDetId Hit = hit->geographicalId(); + if ((Hit.det() == 6) && (Hit.subdetId() == 1) && (Hit.mtdSubDetector() == 1)) { + isBTL = true; + const auto* mtdhit = static_cast(hit); + const auto& hitCluster = mtdhit->mtdCluster(); + if (hitCluster.size() != 0) { + auto recoClusterRef = edmNew::makeRefTo(btlRecCluHandle, &hitCluster); + recoClustersRefs.push_back(recoClusterRef); + } + } } + } //loop over (geometrical) BTL tracks + else { + // --- all ETL tracks (with and without hit in MTD) --- + for (const auto hit : track.recHits()) { + if (hit->isValid() == false) + continue; + MTDDetId Hit = hit->geographicalId(); + if ((Hit.det() == 6) && (Hit.subdetId() == 1) && (Hit.mtdSubDetector() == 2)) { + isETL = true; + + const auto* mtdhit = static_cast(hit); + const auto& hitCluster = mtdhit->mtdCluster(); + if (hitCluster.size() != 0) { + auto recoClusterRef = edmNew::makeRefTo(etlRecCluHandle, &hitCluster); + recoClustersRefs.push_back(recoClusterRef); + } + } + } + } - Ttrack_Associated.push_back(good_association); + LogDebug("MVATrainingNtuple") << "Track p/pt = " << trackGen.p() << " " << trackGen.pt() << " eta " + << trackGen.eta() << " BTL " << isBTL << " ETL " << isETL << " 2disks " + << twoETLdiscs; - // Found TP that is matched to the GTrack + // == TrackingParticle based matching + const reco::TrackBaseRef trkrefb(trackref); + auto tp_info = getAnyMatchedTP(trkrefb); + if (tp_info != nullptr && trkTPSelAll(**tp_info)) { + // == MC truth matching + + auto simClustersRefsIt = tp2SimAssociationMap.find(*tp_info); + withMTD = (simClustersRefsIt != tp2SimAssociationMap.end()); + + // If there is a mtdSimLayerCluster from the tracking particle + if (withMTD) { + // -- Get the refs to MtdSimLayerClusters associated to the TP + std::vector> simClustersRefs; + for (const auto& ref : simClustersRefsIt->val) { + simClustersRefs.push_back(ref); + } + // === ETL + // -- Sort ETL sim clusters by time + if (std::abs(trackGen.eta()) >= etacutBTL_ && !simClustersRefs.empty()) { + std::sort(simClustersRefs.begin(), simClustersRefs.end(), [](const auto& a, const auto& b) { + return a->simLCTime() < b->simLCTime(); + }); + // Check if TP has direct or other sim cluster for BTL + for (const auto& simClusterRef : simClustersRefs) { + if (simClusterRef->trackIdOffset() == 0) { + isTPmtdDirectETL = true; + } else if (simClusterRef->trackIdOffset() != 0) { + isTPmtdOtherETL = true; + } + } + } + // === BTL + // -- Sort BTL sim clusters by time + if (std::abs(trackGen.eta()) < etacutBTL_ && !simClustersRefs.empty()) { + std::sort(simClustersRefs.begin(), simClustersRefs.end(), [](const auto& a, const auto& b) { + return a->simLCTime() < b->simLCTime(); + }); + // Check if TP has direct or other sim cluster for BTL + for (const auto& simClusterRef : simClustersRefs) { + if (simClusterRef->trackIdOffset() == 0) { + isTPmtdDirectBTL = true; + } else if (simClusterRef->trackIdOffset() != 0) { + isTPmtdOtherBTL = true; + } + } + } + + // == Check if the track-cluster association is correct: Track->RecoClus->SimClus == Track->TP->SimClus + for (const auto& recClusterRef : recoClustersRefs) { + if (recClusterRef.isNonnull()) { + auto itp = r2sAssociationMap.equal_range(recClusterRef); + if (itp.first != itp.second) { + auto& simClustersRefs_RecoMatch = (*itp.first).second; + for (const auto& simClusterRef_RecoMatch : simClustersRefs_RecoMatch) { + // Check if simClusterRef_RecoMatch exists in SimClusters + auto simClusterIt = + std::find(simClustersRefs.begin(), simClustersRefs.end(), simClusterRef_RecoMatch); + // SimCluster found in SimClusters + if (simClusterIt != simClustersRefs.end()) { + isTPmtdCorrect = true; + } + } + } + } + } /// end loop over reco clusters associated to this track. + } // --- end "withMTD" + } // TP matching + Ttrack_pt.push_back(trackGen.pt()); + Ttrack_phi.push_back(trackGen.phi()); + Ttrack_eta.push_back(trackGen.eta()); + Ttrack_dz.push_back(std::abs(trackGen.dz(beamSpot.position()))); + Ttrack_dxy.push_back(std::abs(trackGen.dxy(beamSpot.position()))); + Ttrack_chi2.push_back(trackGen.chi2()); + Ttrack_ndof.push_back(trackGen.ndof()); + Ttrack_nValidHits.push_back(trackGen.numberOfValidHits()); + + Ttrack_npixBarrelValidHits.push_back(npixBarrel[trackref]); + Ttrack_npixEndcapValidHits.push_back(npixEndcap[trackref]); + Ttrack_BTLchi2.push_back(btlMatchChi2[trackref]); + Ttrack_BTLtime_chi2.push_back(btlMatchTimeChi2[trackref]); + Ttrack_ETLchi2.push_back(etlMatchChi2[trackref]); + Ttrack_ETLtime_chi2.push_back(etlMatchTimeChi2[trackref]); + + Ttrack_t0.push_back(t0Src[trackref]); + Ttrack_sigmat0.push_back(sigmat0Src[trackref]); + Ttrack_Tmtd.push_back(tMtd[trackref]); + Ttrack_sigmaTmtd.push_back(sigmatimemtd[trackref]); + Ttrack_length.push_back(pathLength[trackref]); + Ttrack_MtdMVA.push_back(mtdQualMVA[trackref]); + Ttrack_lHitPos.push_back(outermostHitPosition[trackref]); + + Ttrack_isBTL.push_back(isBTL); + Ttrack_isETL.push_back(isETL); + Ttrack_withMTD.push_back(withMTD); + + if ((tp_info != nullptr) && (*tp_info)->eventId().bunchCrossing() == 0 && + (*tp_info)->eventId().event() == 0) { // Signal vs PU seperation + Ttrack_Signal.push_back(true); // Signal track + } else { + Ttrack_Signal.push_back(false); // PU track? + } + + if (isTPmtdCorrect) { + Ttrack_Associated.push_back(true); // Associated track + } else { + Ttrack_Associated.push_back(false); // Not associated track + } + if (isTPmtdDirectBTL || isTPmtdDirectETL) { + isTPDirClu = true; // Direct cluster + } + if (isTPmtdOtherBTL || isTPmtdOtherETL) { + isTPOtherClu = true; // Not direct cluster } - } // basic track eta/pT cuts + Ttrack_TPDirClu.push_back(isTPDirClu); + Ttrack_TPOtherClu.push_back(isTPOtherClu); + } // trkRecSel - } // Loop on reco tracks + } // RECO tracks loop BDTtree->Fill(); @@ -825,8 +854,8 @@ void MVATrainingNtuple::fillDescriptions(edm::ConfigurationDescriptions& descrip desc.add("outermostHitPositionSrc", edm::InputTag("trackExtenderWithMTD", "generalTrackOutermostHitPosition")); desc.addUntracked("fileName", "file.root"); - desc.add("ntupleforBDT", false); - desc.add("ntupleforGNN", true); + desc.add("ntupleforBDT", true); + desc.add("ntupleforGNN", false); { edm::ParameterSetDescription psd0; HITrackFilterForPVFinding::fillPSetDescription(psd0); // extension of TrackFilterForPVFinding @@ -836,5 +865,28 @@ void MVATrainingNtuple::fillDescriptions(edm::ConfigurationDescriptions& descrip descriptions.add("mvaTrainingNtuple", desc); } +const bool MVATrainingNtuple::trkTPSelAll(const TrackingParticle& tp) { + bool match = false; + + auto x_pv = tp.parentVertex()->position().x(); + auto y_pv = tp.parentVertex()->position().y(); + auto z_pv = tp.parentVertex()->position().z(); + + auto r_pv = std::sqrt(x_pv * x_pv + y_pv * y_pv); + + match = tp.charge() != 0 && std::abs(tp.eta()) < etacutGEN_ && + ((std::abs(tp.eta()) < etacutBTL_ && tp.pt() > pTcutBTL_) || + (std::abs(tp.eta()) > etacutBTL_ && tp.pt() > pTcutETL_)) && + r_pv < rBTL_ && std::abs(z_pv) < zETL_; + return match; +} + +const bool MVATrainingNtuple::trkRecSel(const reco::TrackBase& trk) { + bool match = false; + match = std::abs(trk.eta()) <= etacutREC_ && ((std::abs(trk.eta()) <= etacutBTL_ && trk.pt() > pTcutBTL_) || + (std::abs(trk.eta()) > etacutBTL_ && trk.pt() > pTcutETL_)); + return match; +} + //define this as a plug-in DEFINE_FWK_MODULE(MVATrainingNtuple); diff --git a/RecoMTD/TimingIDTools/src/MTDTrackQualityMVA.cc b/RecoMTD/TimingIDTools/src/MTDTrackQualityMVA.cc index f3a5cc0db37f8..d69d86ac834dd 100644 --- a/RecoMTD/TimingIDTools/src/MTDTrackQualityMVA.cc +++ b/RecoMTD/TimingIDTools/src/MTDTrackQualityMVA.cc @@ -10,9 +10,11 @@ MTDTrackQualityMVA::MTDTrackQualityMVA(std::string weights_file) { mva_ = std::make_unique(); mva_->initialize(options, method, weights_file, vars_, spec_vars_, true, false); //use GBR, GradBoost -} + //mva_->initialize(options, method, weights_file, vars_, spec_vars_, false, false); //use TMVA +}; float MTDTrackQualityMVA::operator()(const reco::TrackRef& trk, + const reco::BeamSpot& beamspot, const edm::ValueMap& npixBarrels, const edm::ValueMap& npixEndcaps, const edm::ValueMap& btl_chi2s, @@ -20,30 +22,44 @@ float MTDTrackQualityMVA::operator()(const reco::TrackRef& trk, const edm::ValueMap& etl_chi2s, const edm::ValueMap& etl_time_chi2s, const edm::ValueMap& tmtds, - const edm::ValueMap& trk_lengths) const { + const edm::ValueMap& sigmatmtds, + const edm::ValueMap& trk_lengths, + const edm::ValueMap& trk_lhitpos) const { std::map vars; - //---training performed only above 0.5 GeV - constexpr float minPtForMVA = 0.5; - if (trk->pt() < minPtForMVA) + static constexpr double etacutREC_ = 3.; // |eta| < 3 + static constexpr double etacutBTL_ = 1.5; // |eta| < 1.5 for BTL + static constexpr double pTcutBTL_ = 0.7; // PT > 0.7 GeV + static constexpr double pTcutETL_ = 0.2; // PT > 0.2 GeV + + //---training performed only for the specified cuts + if (std::abs(trk->eta()) > etacutREC_) // max eta cut + return -1; + if (std::abs(trk->eta()) <= etacutBTL_ && trk->pt() < pTcutBTL_) // min PT cut for BTL + return -1; + if (std::abs(trk->eta()) > etacutBTL_ && trk->pt() < pTcutETL_) // min PT cut for ETL return -1; //---training performed only for tracks with MTD hits if (tmtds[trk] > 0) { - vars.emplace(vars_[int(VarID::pt)], trk->pt()); - vars.emplace(vars_[int(VarID::eta)], trk->eta()); - vars.emplace(vars_[int(VarID::phi)], trk->phi()); - vars.emplace(vars_[int(VarID::chi2)], trk->chi2()); - vars.emplace(vars_[int(VarID::ndof)], trk->ndof()); - vars.emplace(vars_[int(VarID::numberOfValidHits)], trk->numberOfValidHits()); - vars.emplace(vars_[int(VarID::numberOfValidPixelBarrelHits)], npixBarrels[trk]); - vars.emplace(vars_[int(VarID::numberOfValidPixelEndcapHits)], npixEndcaps[trk]); - vars.emplace(vars_[int(VarID::btlMatchChi2)], btl_chi2s.contains(trk.id()) ? btl_chi2s[trk] : -1); - vars.emplace(vars_[int(VarID::btlMatchTimeChi2)], btl_time_chi2s.contains(trk.id()) ? btl_time_chi2s[trk] : -1); - vars.emplace(vars_[int(VarID::etlMatchChi2)], etl_chi2s.contains(trk.id()) ? etl_chi2s[trk] : -1); - vars.emplace(vars_[int(VarID::etlMatchTimeChi2)], etl_time_chi2s.contains(trk.id()) ? etl_time_chi2s[trk] : -1); - vars.emplace(vars_[int(VarID::mtdt)], tmtds[trk]); - vars.emplace(vars_[int(VarID::path_len)], trk_lengths[trk]); + vars.emplace(vars_[int(VarID::Track_pt)], trk->pt()); + vars.emplace(vars_[int(VarID::Track_eta)], trk->eta()); + vars.emplace(vars_[int(VarID::Track_phi)], trk->phi()); + vars.emplace(vars_[int(VarID::Track_dz)], trk->dz(beamspot.position())); + vars.emplace(vars_[int(VarID::Track_dxy)], trk->dxy(beamspot.position())); + vars.emplace(vars_[int(VarID::Track_chi2)], trk->chi2()); + vars.emplace(vars_[int(VarID::Track_ndof)], trk->ndof()); + vars.emplace(vars_[int(VarID::Track_npixBarrelValidHits)], npixBarrels[trk]); + vars.emplace(vars_[int(VarID::Track_npixEndcapValidHits)], npixEndcaps[trk]); + vars.emplace(vars_[int(VarID::Track_BTLchi2)], btl_chi2s.contains(trk.id()) ? btl_chi2s[trk] : -1); + vars.emplace(vars_[int(VarID::Track_BTLtime_chi2)], btl_time_chi2s.contains(trk.id()) ? btl_time_chi2s[trk] : -1); + vars.emplace(vars_[int(VarID::Track_ETLchi2)], etl_chi2s.contains(trk.id()) ? etl_chi2s[trk] : -1); + vars.emplace(vars_[int(VarID::Track_ETLtime_chi2)], etl_time_chi2s.contains(trk.id()) ? etl_time_chi2s[trk] : -1); + vars.emplace(vars_[int(VarID::Track_Tmtd)], tmtds[trk]); + vars.emplace(vars_[int(VarID::Track_sigmaTmtd)], sigmatmtds[trk]); + vars.emplace(vars_[int(VarID::Track_length)], trk_lengths[trk]); + vars.emplace(vars_[int(VarID::Track_lHitPos)], trk_lhitpos[trk]); + return 1. / (1 + sqrt(2 / (1 + mva_->evaluate(vars, false)) - 1)); //return values between 0-1 (probability) } else return -1;