Add track selection MVA plots to tracking DQM (92X)

DQM/TrackingMonitor/interface/TrackBuildingAnalyzer.h

@@ -34,6 +34,9 @@ class DQMStore;
 class TrackBuildingAnalyzer 
 {
     public:
+        using MVACollection = std::vector<float>;
+        using QualityMaskCollection = std::vector<unsigned char>;
+
         TrackBuildingAnalyzer(const edm::ParameterSet&);
         ~TrackBuildingAnalyzer();
         void initHisto(DQMStore::IBooker & ibooker, const edm::ParameterSet&);
@@ -55,6 +58,12 @@ class TrackBuildingAnalyzer
             const edm::ESHandle<MagneticField>& theMF,
             const edm::ESHandle<TransientTrackingRecHitBuilder>& theTTRHBuilder
         );
+        void analyze
+        (
+            const edm::View<reco::Track>& trackCollection,
+            const std::vector<const MVACollection *>& mvaCollections,
+            const std::vector<const QualityMaskCollection *>& qualityMaskCollections
+        );
 
     private:
 
@@ -92,6 +101,13 @@ class TrackBuildingAnalyzer
 	MonitorElement* stoppingSource;
 	MonitorElement* stoppingSourceVSeta;
 	MonitorElement* stoppingSourceVSphi;
+
+	std::vector<MonitorElement *> trackMVAs;
+	std::vector<MonitorElement *> trackMVAsHP;
+	std::vector<MonitorElement *> trackMVAsVsPtProfile;
+	std::vector<MonitorElement *> trackMVAsHPVsPtProfile;
+	std::vector<MonitorElement *> trackMVAsVsEtaProfile;
+	std::vector<MonitorElement *> trackMVAsHPVsEtaProfile;
 
         std::string histname;  //for naming the histograms according to algorithm used
 
@@ -112,5 +128,6 @@ class TrackBuildingAnalyzer
 	bool doProfPHI;
 	bool doProfETA;
 	bool doStopSource;
+	bool doMVAPlots;
 };
 #endif

DQM/TrackingMonitor/interface/TrackingMonitor.h

@@ -51,6 +51,9 @@ class GenericTriggerEventFlag;
 class TrackingMonitor : public DQMEDAnalyzer 
 {
     public:
+        using MVACollection = std::vector<float>;
+        using QualityMaskCollection = std::vector<unsigned char>;
+
         explicit TrackingMonitor(const edm::ParameterSet&);
         ~TrackingMonitor();
         virtual void beginJob(void);
@@ -95,6 +98,9 @@ class TrackingMonitor : public DQMEDAnalyzer
 	edm::EDGetTokenT<edmNew::DetSetVector<SiStripCluster> > stripClustersToken_;
 	edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster> > pixelClustersToken_;
 
+	std::vector<std::tuple<edm::EDGetTokenT<MVACollection>, edm::EDGetTokenT<QualityMaskCollection> > > mvaQualityTokens_;
+	edm::EDGetTokenT<edm::View<reco::Track> > mvaTrackToken_;
+
 	std::string Quality_;
 	std::string AlgoName_;
 
@@ -181,6 +187,7 @@ class TrackingMonitor : public DQMEDAnalyzer
 	bool doGeneralPropertiesPlots_;
 	bool doHitPropertiesPlots_;
 	bool doTkCandPlots;
+	bool doMVAPlots;
 	bool doSeedNumberPlot;
 	bool doSeedLumiAnalysis_;
 	bool doSeedVsClusterPlot;

DQM/TrackingMonitor/python/MonitorTrackSTAMuons_cfi.py

@@ -12,10 +12,8 @@
 MonitorTrackSTAMuons.ClusterLabels = cms.vstring('Tot')
 
 # output parameters
-MonitorTrackSTAMuons.OutputMEsInRootFile = cms.bool(False)
 MonitorTrackSTAMuons.AlgoName            = cms.string('sta')
 MonitorTrackSTAMuons.Quality             = cms.string('')
-MonitorTrackSTAMuons.OutputFileName      = cms.string('monitortrackparameters_stamuons.root')
 MonitorTrackSTAMuons.FolderName          = cms.string('Muons/standAloneMuonsUpdatedAtVtx')
 MonitorTrackSTAMuons.BSFolderName        = cms.string('Muons/standAloneMuonsUpdatedAtVtx/BeamSpotParameters')
 

DQM/TrackingMonitor/python/TrackingMonitorAllTrackingSequences_cff.py

@@ -6,8 +6,6 @@
 from DQM.TrackingMonitor.TrackingMonitor_cfi import *
 
 # properties
-TrackMon.OutputMEsInRootFile    = cms.bool(False)
-TrackMon.OutputFileName         = cms.string('TrackingMonitorAllSequences.root')
 TrackMon.MeasurementState       = cms.string('ImpactPoint')
 
 # which plots to do

DQM/TrackingMonitor/python/TrackingMonitorSeed_cfi.py

@@ -4,8 +4,6 @@
 
 TrackMonSeed = DQM.TrackingMonitor.TrackingMonitor_cfi.TrackMon.clone()
 
-TrackMonSeed.OutputMEsInRootFile        = cms.bool(False)
-TrackMonSeed.OutputFileName             = cms.string('TrackingMonitorSeedMultiplicity.root')
 TrackMonSeed.MeasurementState           = cms.string('ImpactPoint')
 TrackMonSeed.FolderName                 = cms.string('Tracking/TrackParameters')
 TrackMonSeed.BSFolderName               = cms.string('Tracking/TrackParameters/BeamSpotParameters')

DQM/TrackingMonitor/python/TrackingMonitor_cfi.py

@@ -11,6 +11,8 @@
     TrackProducer    = cms.InputTag("generalTracks"),
     SeedProducer     = cms.InputTag("initialStepSeeds"),
     TCProducer       = cms.InputTag("initialStepTrackCandidates"),
+    MVAProducers     = cms.vstring("initialStepClassifier1", "initialStepClassifier2"),
+    TrackProducerForMVA = cms.InputTag("initialStepTracks"),
     ClusterLabels    = cms.vstring('Tot'), # to decide which Seeds-Clusters correlation plots to have default is Total other options 'Strip', 'Pix'
     beamSpot         = cms.InputTag("offlineBeamSpot"),
     primaryVertex    = cms.InputTag('offlinePrimaryVertices'),
@@ -31,10 +33,8 @@
     pvLabels = cms.vstring(),
 
     # output parameters
-    OutputMEsInRootFile = cms.bool(False),
     AlgoName            = cms.string('GenTk'),
     Quality             = cms.string(''),
-    OutputFileName      = cms.string('MonitorTrack.root'),
     FolderName          = cms.string('Tracking/GlobalParameters'),
     BSFolderName        = cms.string('Tracking/ParametersVsBeamSpot'),
     PVFolderName        = cms.string('Tracking/PrimaryVertices'),
@@ -89,6 +89,7 @@
     doPlotsVsLUMI                       = cms.bool(False),
     doPlotsVsBX                         = cms.bool(False),
     doHIPlots                           = cms.bool(False),                              
+    doMVAPlots                          = cms.bool(False),
     qualityString = cms.string("highPurity"),                      
     #which seed plots to do
     doSeedNumberHisto = cms.bool(False),
@@ -361,6 +362,11 @@
     TCDzMax = cms.double(400.0),
     TCDzMin = cms.double(-400.0),                                                
 
+    # Track selection MVA
+    MVABin  = cms.int32(100),
+    MVAMin  = cms.double(-1),
+    MVAMax  = cms.double(1),
+
 #######################################
 ## needed for tracksVScluster and seedVScluster
 

DQM/TrackingMonitor/src/TrackBuildingAnalyzer.cc

@@ -96,10 +96,15 @@ void TrackBuildingAnalyzer::initHisto(DQMStore::IBooker & ibooker, const edm::Pa
   int    TCHitBin   = iConfig.getParameter<int>(   "TCHitBin");
   double TCHitMin   = iConfig.getParameter<double>("TCHitMin");
   double TCHitMax   = iConfig.getParameter<double>("TCHitMax");
+
+  int MVABin        = iConfig.getParameter<int>(   "MVABin");
+  double MVAMin     = iConfig.getParameter<double>("MVAMin");
+  double MVAMax     = iConfig.getParameter<double>("MVAMax");
 
 
   edm::InputTag seedProducer   = iConfig.getParameter<edm::InputTag>("SeedProducer");
   edm::InputTag tcProducer     = iConfig.getParameter<edm::InputTag>("TCProducer");
+  std::vector<std::string> mvaProducers = iConfig.getParameter<std::vector<std::string> >("MVAProducers");
 
   doAllPlots     = iConfig.getParameter<bool>("doAllPlots");
   doAllSeedPlots = iConfig.getParameter<bool>("doSeedParameterHistos");
@@ -117,6 +122,7 @@ void TrackBuildingAnalyzer::initHisto(DQMStore::IBooker & ibooker, const edm::Pa
   doProfPHI      = iConfig.getParameter<bool>("doSeedNVsPhiProf");
   doProfETA      = iConfig.getParameter<bool>("doSeedNVsEtaProf");
   doStopSource   = iConfig.getParameter<bool>("doStopSource");
+  doMVAPlots     = iConfig.getParameter<bool>("doMVAPlots");
 
   //    if (doAllPlots){doAllSeedPlots=true; doTCPlots=true;}
 
@@ -319,6 +325,52 @@ void TrackBuildingAnalyzer::initHisto(DQMStore::IBooker & ibooker, const edm::Pa
     TrackCandPhiVsEta = ibooker.book2D(histname+CatagoryName, histname+CatagoryName, EtaBin, EtaMin, EtaMax, PhiBin, PhiMin, PhiMax);
     TrackCandPhiVsEta->setAxisTitle("Track Candidate #eta", 1);
     TrackCandPhiVsEta->setAxisTitle("Track Candidate #phi", 2);
+
+    if(doAllTCPlots || doMVAPlots) {
+      for(size_t i=1, end=mvaProducers.size(); i<=end; ++i) {
+        auto num = std::to_string(i);
+        std::string pfix;
+
+        if(i == 1) {
+          trackMVAsHP.push_back(nullptr);
+          trackMVAsHPVsPtProfile.push_back(nullptr);
+          trackMVAsHPVsEtaProfile.push_back(nullptr);
+        }
+        else {
+          pfix = " (not loose-selected)";
+          std::string pfix2 = " (not HP-selected)";
+          histname = "TrackMVA"+num+"HP_"+tcProducer.label() + "_";
+          trackMVAsHP.push_back(ibooker.book1D(histname+CatagoryName, histname+CatagoryName+pfix2, MVABin, MVAMin, MVAMax));
+          trackMVAsHP.back()->setAxisTitle("Track selection MVA"+num, 1);
+          trackMVAsHP.back()->setAxisTitle("Number of tracks", 2);
+
+          histname = "TrackMVA"+num+"HPVsPtProfile_"+tcProducer.label() + "_";
+          trackMVAsHPVsPtProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix2, TrackPtBin, TrackPtMin, TrackPtMax, MVABin, MVAMin, MVAMax));
+          trackMVAsHPVsPtProfile.back()->setAxisTitle("Track p_{T} (GeV/c)", 1);
+          trackMVAsHPVsPtProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
+
+          histname = "TrackMVA"+num+"HPVsEtaProfile_"+tcProducer.label() + "_";
+          trackMVAsHPVsEtaProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix2, EtaBin, EtaMin, EtaMax, MVABin, MVAMin, MVAMax));
+          trackMVAsHPVsEtaProfile.back()->setAxisTitle("Track #eta", 1);
+          trackMVAsHPVsEtaProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
+        }
+
+        histname = "TrackMVA"+num+"_"+tcProducer.label() + "_";
+        trackMVAs.push_back(ibooker.book1D(histname+CatagoryName, histname+CatagoryName+pfix, MVABin, MVAMin, MVAMax));
+        trackMVAs.back()->setAxisTitle("Track selection MVA"+num, 1);
+        trackMVAs.back()->setAxisTitle("Number of tracks", 2);
+
+        histname = "TrackMVA"+num+"VsPtProfile_"+tcProducer.label() + "_";
+        trackMVAsVsPtProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix, TrackPtBin, TrackPtMin, TrackPtMax, MVABin, MVAMin, MVAMax));
+        trackMVAsVsPtProfile.back()->setAxisTitle("Track p_{T} (GeV/c)", 1);
+        trackMVAsVsPtProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
+
+        histname = "TrackMVA"+num+"VsEtaProfile_"+tcProducer.label() + "_";
+        trackMVAsVsEtaProfile.push_back(ibooker.bookProfile(histname+CatagoryName, histname+CatagoryName+pfix, EtaBin, EtaMin, EtaMax, MVABin, MVAMin, MVAMax));
+        trackMVAsVsEtaProfile.back()->setAxisTitle("Track #eta", 1);
+        trackMVAsVsEtaProfile.back()->setAxisTitle("Track selection MVA"+num, 2);
+      }
+    }
   }
 
 }
@@ -443,3 +495,64 @@ void TrackBuildingAnalyzer::analyze
     if (doAllTCPlots) NumberOfRecHitsPerTrackCandVsPhiProfile->Fill( phi, numberOfHits );
   }
 }
+
+namespace {
+  bool trackSelected(unsigned char mask, unsigned char qual) {
+    return mask & 1<<qual;
+  }
+}
+void TrackBuildingAnalyzer::analyze(const edm::View<reco::Track>& trackCollection,
+                                    const std::vector<const MVACollection *>& mvaCollections,
+                                    const std::vector<const QualityMaskCollection *>& qualityMaskCollections) {
+  if(!(doAllTCPlots || doMVAPlots))
+    return;
+  if(trackCollection.empty())
+    return;
+
+  const auto ntracks = trackCollection.size();
+  const auto nmva = mvaCollections.size();
+  for(const auto mva: mvaCollections) {
+    if(mva->size() != ntracks) {
+      edm::LogError("LogicError") << "TrackBuildingAnalyzer: Incompatible size of MVACollection, " << mva->size() << " differs from the size of the track collection " << ntracks;
+      return;
+    }
+  }
+  for(const auto qual: qualityMaskCollections) {
+    if(qual->size() != ntracks) {
+      edm::LogError("LogicError") << "TrackBuildingAnalyzer: Incompatible size of QualityMaskCollection, " << qual->size() << " differs from the size of the track collection " << ntracks;
+      return;
+    }
+  }
+
+
+  for(size_t iTrack=0; iTrack<ntracks; ++iTrack) {
+    // Fill MVA1 histos with all tracks, MVA2 histos only with tracks
+    // not selected by MVA1 etc
+    bool selectedLoose = false;
+    bool selectedHP = false;
+
+    const auto pt = trackCollection[iTrack].pt();
+    const auto eta = trackCollection[iTrack].eta();
+
+    for(size_t iMVA=0; iMVA<nmva; ++iMVA) {
+      const auto mva = (*(mvaCollections[iMVA]))[iTrack];
+      if(!selectedLoose) {
+        trackMVAs[iMVA]->Fill(mva);
+        trackMVAsVsPtProfile[iMVA]->Fill(pt, mva);
+        trackMVAsVsEtaProfile[iMVA]->Fill(eta, mva);
+      }
+      if(iMVA >= 1 && !selectedHP) {
+        trackMVAsHP[iMVA]->Fill(mva);
+        trackMVAsHPVsPtProfile[iMVA]->Fill(pt, mva);
+        trackMVAsHPVsEtaProfile[iMVA]->Fill(eta, mva);
+      }
+
+      const auto qual = (*(qualityMaskCollections)[iMVA])[iTrack];
+      selectedLoose |= trackSelected(qual, reco::TrackBase::loose);
+      selectedHP |= trackSelected(qual, reco::TrackBase::highPurity);
+
+      if(selectedLoose && selectedHP)
+        break;
+    }
+  }
+}

DQM/TrackingMonitor/src/TrackingMonitor.cc

@@ -13,6 +13,7 @@
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/Utilities/interface/InputTag.h"
+#include "FWCore/Utilities/interface/transform.h"
 #include "DQMServices/Core/interface/DQMStore.h"
 #include "DQM/TrackingMonitor/interface/TrackBuildingAnalyzer.h"
 #include "DQM/TrackingMonitor/interface/TrackAnalyzer.h"
@@ -119,6 +120,16 @@ TrackingMonitor::TrackingMonitor(const edm::ParameterSet& iConfig)
   trackCandidateToken_ = consumes<TrackCandidateCollection>(tcProducer); 
   seedToken_           = consumes<edm::View<TrajectorySeed> >(seedProducer);
 
+  doMVAPlots = iConfig.getParameter<bool>("doMVAPlots");
+  if(doMVAPlots) {
+    mvaQualityTokens_ = edm::vector_transform(iConfig.getParameter<std::vector<std::string> >("MVAProducers"),
+                                              [&](const std::string& tag) {
+                                                return std::make_tuple(consumes<MVACollection>(edm::InputTag(tag, "MVAValues")),
+                                                                       consumes<QualityMaskCollection>(edm::InputTag(tag, "QualityMasks")));
+                                              });
+    mvaTrackToken_ = consumes<edm::View<reco::Track> >(iConfig.getParameter<edm::InputTag>("TrackProducerForMVA"));
+  }
+
   edm::InputTag stripClusterInputTag_ = iConfig.getParameter<edm::InputTag>("stripCluster");
   edm::InputTag pixelClusterInputTag_ = iConfig.getParameter<edm::InputTag>("pixelCluster");
   stripClustersToken_ = mayConsume<edmNew::DetSetVector<SiStripCluster> > (stripClusterInputTag_);
@@ -852,6 +863,26 @@ void TrackingMonitor::analyze(const edm::Event& iEvent, const edm::EventSetup& i
 	} else {
 	  edm::LogWarning("TrackingMonitor") << "No Track Candidates in the event.  Not filling associated histograms";
 	}
+
+	if(doMVAPlots) {
+	  // Get MVA and quality mask collections
+	  std::vector<const MVACollection *> mvaCollections;
+	  std::vector<const QualityMaskCollection *> qualityMaskCollections;
+
+	  edm::Handle<edm::View<reco::Track> > htracks;
+	  iEvent.getByToken(mvaTrackToken_, htracks);
+
+	  edm::Handle<MVACollection> hmva;
+	  edm::Handle<QualityMaskCollection> hqual;
+	  for(const auto& tokenTpl: mvaQualityTokens_) {
+	    iEvent.getByToken(std::get<0>(tokenTpl), hmva);
+	    iEvent.getByToken(std::get<1>(tokenTpl), hqual);
+
+	    mvaCollections.push_back(hmva.product());
+	    qualityMaskCollections.push_back(hqual.product());
+	  }
+	  theTrackBuildingAnalyzer->analyze(*htracks, mvaCollections, qualityMaskCollections);
+	}
       }
 
       //plots for trajectory seeds

DQM/TrackingMonitor/test/TrackingMonitor_StandAlone_cfg.py

@@ -76,8 +76,6 @@
 
 # properties
 process.TrackMon.AlgoName               = cms.string('GenTk')
-process.TrackMon.OutputFileName 	= cms.string('TrackingMonitor.root')
-process.TrackMon.OutputMEsInRootFile 	= cms.bool(True)
 
 process.TrackMon.FolderName             = cms.string('Track/GlobalParameters')
 process.TrackMon.BSFolderName           = cms.string('Track/BeamSpotParameters')

DQM/TrackingMonitorSource/python/TrackingSourceConfig_Tier0_cff.py

@@ -1,5 +1,6 @@
 import FWCore.ParameterSet.Config as cms
 import RecoTracker.IterativeTracking.iterativeTkConfig as _cfg
+import RecoTracker.IterativeTracking.iterativeTkUtils as _utils
 
 ### load which are the tracks collection 2 be monitored
 from DQM.TrackingMonitorSource.TrackCollections2monitor_cff import *
@@ -310,8 +311,14 @@
     TrackingDQMSourceTier0 += cms.ignore(locals()[label])
 # seeding monitoring
 for _eraName, _postfix, _era in _cfg.allEras():
+    mvaSel = _utils.getMVASelectors(_postfix)
     _seq = cms.Sequence()
     for step in locals()["selectedIterTrackingStep"+_postfix]:
+        if step in mvaSel:
+            locals()["TrackSeedMon"+step].doMVAPlots = True
+            locals()["TrackSeedMon"+step].TrackProducerForMVA = mvaSel[step][0]
+            locals()["TrackSeedMon"+step].MVAProducers = mvaSel[step][1]
+
         _seq += locals()["TrackSeedMon"+step]
     if _eraName == "":
         locals()["TrackSeedMonSequence"] = _seq

RecoTracker/IterativeTracking/python/iterativeTkUtils.py

@@ -0,0 +1,32 @@
+# This file provides additional helpers for getting information of
+# iterations in automated way.
+import RecoTracker.IterativeTracking.iterativeTkConfig as _cfg
+import RecoTracker.IterativeTracking.iterativeTk_cff as _iterativeTk_cff
+
+def getMVASelectors(postfix):
+    # assume naming convention that the iteration name (when first
+    # letter in lower case) is the selector name
+
+    ret = {}
+
+    for iterName, seqName in _cfg.iterationAlgos(postfix, includeSequenceName=True):
+        if hasattr(_iterativeTk_cff, iterName):
+            mod = getattr(_iterativeTk_cff, iterName)
+            seq = getattr(_iterativeTk_cff, seqName)
+
+            # Ignore iteration if the MVA selector module is not in the sequence
+            try:
+                seq.index(mod)
+            except:
+                continue
+
+            typeName = mod._TypedParameterizable__type
+            classifiers = []
+            if typeName == "ClassifierMerger":
+                classifiers = mod.inputClassifiers.value()
+            elif "TrackMVAClassifier" in typeName:
+                classifiers = [iterName]
+            if len(classifiers) > 0:
+                ret[iterName] = (iterName+"Tracks", classifiers)
+
+    return ret