RecoVertex & RecoMTD: remove first iteration of 4D and use 3D with ti…

…ming, update tofPIDProducer, update 4D time calculation

Configuration/ProcessModifiers/python/useMVASel_cff.py

@@ -0,0 +1,6 @@
+import FWCore.ParameterSet.Config as cms
+
+# Modifier to enable the use of the MVA selection on 
+# tracks for the 4D vertex reco
+
+useMVASel =  cms.Modifier()

RecoMTD/TimingIDTools/plugins/TOFPIDProducer.cc

@@ -49,11 +49,18 @@ class TOFPIDProducer : public edm::stream::EDProducer<> {
   edm::EDGetTokenT<edm::ValueMap<float>> tofkToken_;
   edm::EDGetTokenT<edm::ValueMap<float>> tofpToken_;
   edm::EDGetTokenT<reco::VertexCollection> vtxsToken_;
+  edm::EDGetTokenT<edm::ValueMap<float>> trackMTDTimeQualityToken_;
   double vtxMaxSigmaT_;
   double maxDz_;
   double maxDtSignificance_;
   double minProbHeavy_;
   double fixedT0Error_;
+  double probPion_;
+  double probKaon_;
+  double probProton_;
+  double minTrackTimeQuality_;
+  bool MVASel_;
+  bool vertexReassignment_;
 };
 
 TOFPIDProducer::TOFPIDProducer(const ParameterSet& iConfig)
@@ -65,11 +72,19 @@ TOFPIDProducer::TOFPIDProducer(const ParameterSet& iConfig)
       tofkToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofkSrc"))),
       tofpToken_(consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("tofpSrc"))),
       vtxsToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vtxsSrc"))),
+      trackMTDTimeQualityToken_(
+          consumes<edm::ValueMap<float>>(iConfig.getParameter<edm::InputTag>("trackMTDTimeQualityVMapTag"))),
       vtxMaxSigmaT_(iConfig.getParameter<double>("vtxMaxSigmaT")),
       maxDz_(iConfig.getParameter<double>("maxDz")),
       maxDtSignificance_(iConfig.getParameter<double>("maxDtSignificance")),
       minProbHeavy_(iConfig.getParameter<double>("minProbHeavy")),
-      fixedT0Error_(iConfig.getParameter<double>("fixedT0Error")) {
+      fixedT0Error_(iConfig.getParameter<double>("fixedT0Error")),
+      probPion_(iConfig.getParameter<double>("probPion")),
+      probKaon_(iConfig.getParameter<double>("probKaon")),
+      probProton_(iConfig.getParameter<double>("probProton")),
+      minTrackTimeQuality_(iConfig.getParameter<double>("minTrackTimeQuality")),
+      MVASel_(iConfig.getParameter<bool>("MVASel")),
+      vertexReassignment_(iConfig.getParameter<bool>("vertexReassignment")) {
   produces<edm::ValueMap<float>>(t0Name);
   produces<edm::ValueMap<float>>(sigmat0Name);
   produces<edm::ValueMap<float>>(t0safeName);
@@ -97,6 +112,8 @@ void TOFPIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio
       ->setComment("Input ValueMap for track tof as proton");
   desc.add<edm::InputTag>("vtxsSrc", edm::InputTag("unsortedOfflinePrimaryVertices4DwithPID"))
       ->setComment("Input primary vertex collection");
+  desc.add<edm::InputTag>("trackMTDTimeQualityVMapTag", edm::InputTag("mtdTrackQualityMVA:mtdQualMVA"))
+      ->setComment("Track MVA quality value");
   desc.add<double>("vtxMaxSigmaT", 0.025)
       ->setComment("Maximum primary vertex time uncertainty for use in particle id [ns]");
   desc.add<double>("maxDz", 0.1)
@@ -107,6 +124,12 @@ void TOFPIDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptio
   desc.add<double>("minProbHeavy", 0.75)
       ->setComment("Minimum probability for a particle to be a kaon or proton before reassigning the timestamp");
   desc.add<double>("fixedT0Error", 0.)->setComment("Use a fixed T0 uncertainty [ns]");
+  desc.add<double>("probPion", 1.)->setComment("A priori probability pions");
+  desc.add<double>("probKaon", 1.)->setComment("A priori probability kaons");
+  desc.add<double>("probProton", 1.)->setComment("A priori probability for protons");
+  desc.add<double>("minTrackTimeQuality", 0.8)->setComment("Minimum MVA Quality selection on tracks");
+  desc.add<bool>("MVASel", false)->setComment("Use MVA Quality selection");
+  desc.add<bool>("vertexReassignment", true)->setComment("Track-vertex reassignment");
 
   descriptions.add("tofPIDProducer", desc);
 }
@@ -142,6 +165,8 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
 
   const auto& vtxs = ev.get(vtxsToken_);
 
+  const auto& trackMVAQualIn = ev.get(trackMTDTimeQualityToken_);
+
   //output value maps (PID probabilities and recalculated time at beamline)
   std::vector<float> t0OutRaw;
   std::vector<float> sigmat0OutRaw;
@@ -165,7 +190,9 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
     float prob_k = -1.;
     float prob_p = -1.;
 
-    if (sigmat0 > 0.) {
+    float trackMVAQual = trackMVAQualIn[trackref];
+
+    if (sigmat0 > 0. && (!MVASel_ || (MVASel_ && trackMVAQual >= minTrackTimeQuality_))) {
       double rsigmazsq = 1. / track.dzError() / track.dzError();
       double rsigmat = 1. / sigmatmtd;
 
@@ -239,7 +266,12 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
         double chisqmin_k = std::numeric_limits<double>::max();
         double chisqmin_p = std::numeric_limits<double>::max();
         //loop through vertices and check for better matches
-        for (const reco::Vertex& vtx : vtxs) {
+        for (unsigned int ivtx = 0; ivtx < vtxs.size(); ++ivtx) {
+          const reco::Vertex& vtx = vtxs[ivtx];
+          if (!vertexReassignment_) {
+            if (ivtx != (unsigned int)vtxidx)
+              continue;
+          }
           if (!(vtx.tError() > 0. && vtx.tError() < vtxMaxSigmaT_)) {
             continue;
           }
@@ -283,9 +315,9 @@ void TOFPIDProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
 
         //compute PID probabilities
         //*TODO* deal with heavier nucleons and/or BSM case here?
-        double rawprob_pi = exp(-0.5 * chisqmin_pi);
-        double rawprob_k = exp(-0.5 * chisqmin_k);
-        double rawprob_p = exp(-0.5 * chisqmin_p);
+        double rawprob_pi = probPion_ * exp(-0.5 * chisqmin_pi);
+        double rawprob_k = probKaon_ * exp(-0.5 * chisqmin_k);
+        double rawprob_p = probProton_ * exp(-0.5 * chisqmin_p);
 
         double normprob = 1. / (rawprob_pi + rawprob_k + rawprob_p);
 

RecoVertex/Configuration/python/RecoVertex_EventContent_cff.py

@@ -17,9 +17,7 @@
                                             'keep *_offlinePrimaryVertices4DWithBS__*',
                                             'keep *_trackTimeValueMapProducer_*_*' ]
 
-_phase2_tktiming_layer_RecoVertexEventContent = [ 'keep *_offlinePrimaryVertices4DnoPID__*',
-                                                  'keep *_offlinePrimaryVertices4DnoPIDWithBS__*',
-                                                  'keep *_tofPID_*_*']
+_phase2_tktiming_layer_RecoVertexEventContent = [ 'keep *_tofPID_*_*']
 phase2_timing.toModify( RecoVertexAOD,
      outputCommands = RecoVertexAOD.outputCommands + _phase2_tktiming_RecoVertexEventContent)
 phase2_timing_layer.toModify( RecoVertexAOD,

RecoVertex/Configuration/python/RecoVertex_cff.py

@@ -45,15 +45,11 @@
 from RecoVertex.Configuration.RecoVertex_phase2_timing_cff import (tpClusterProducer ,
                                                                   quickTrackAssociatorByHits ,
                                                                   trackTimeValueMapProducer ,
-                                                                  unsortedOfflinePrimaryVertices4DnoPID ,
-                                                                  trackWithVertexRefSelectorBeforeSorting4DnoPID ,
-                                                                  trackRefsForJetsBeforeSorting4DnoPID ,
-                                                                  offlinePrimaryVertices4DnoPID ,
-                                                                  offlinePrimaryVertices4DnoPIDWithBS,
                                                                   unsortedOfflinePrimaryVertices4DwithPID ,
                                                                   offlinePrimaryVertices4DwithPID ,
                                                                   offlinePrimaryVertices4DwithPIDWithBS,
                                                                   tofPID,
+                                                                  tofPID3D,
                                                                   tofPID4DnoPID,
                                                                   unsortedOfflinePrimaryVertices4D,
                                                                   trackWithVertexRefSelectorBeforeSorting4D,
@@ -73,13 +69,9 @@
                                             )
 
 _phase2_tktiming_layer_vertexrecoTask = cms.Task( _phase2_tktiming_vertexrecoTask.copy() ,
-                                            unsortedOfflinePrimaryVertices4DnoPID ,
-                                            trackWithVertexRefSelectorBeforeSorting4DnoPID ,
-                                            trackRefsForJetsBeforeSorting4DnoPID ,
-                                            offlinePrimaryVertices4DnoPID ,
-                                            offlinePrimaryVertices4DnoPIDWithBS,
+                                            tofPID3D,
                                             tofPID,
-                                            tofPID4DnoPID,
+                                            tofPID4DnoPID
                                             )
 
 from Configuration.Eras.Modifier_phase2_timing_cff import phase2_timing
@@ -93,3 +85,6 @@
 phase2_timing_layer.toModify(offlinePrimaryVertices4D, vertices = "unsortedOfflinePrimaryVertices4D", particles = "trackRefsForJetsBeforeSorting4D")
 phase2_timing_layer.toModify(offlinePrimaryVertices4DWithBS, vertices = "unsortedOfflinePrimaryVertices4D:WithBS", particles = "trackRefsForJetsBeforeSorting4D")
 
+from Configuration.ProcessModifiers.useMVASel_cff import useMVASel
+useMVASel.toModify(unsortedOfflinePrimaryVertices4D, useMVACut = cms.bool(True))
+useMVASel.toModify(unsortedOfflinePrimaryVertices4DwithPID, useMVACut = cms.bool(True))

RecoVertex/Configuration/python/RecoVertex_phase2_timing_cff.py

@@ -7,6 +7,30 @@
     TkClusParameters = DA2D_vectParameters,
     TrackTimesLabel = cms.InputTag("trackTimeValueMapProducer","generalTracksConfigurableFlatResolutionModel"),
     TrackTimeResosLabel = cms.InputTag("trackTimeValueMapProducer","generalTracksConfigurableFlatResolutionModelResolution"),
+    vertexCollections = cms.VPSet(
+        cms.PSet(
+            algorithm = cms.string('AdaptiveVertexFitter'),
+            chi2cutoff = cms.double(2.5),
+            label = cms.string(''),
+            maxDistanceToBeam = cms.double(1.0),
+            minNdof = cms.double(0.0),
+            useBeamConstraint = cms.bool(False),
+            vertexTimeParameters = cms.PSet(
+                algorithm = cms.string('legacy4D')
+            )
+        ),
+        cms.PSet(
+            algorithm = cms.string('AdaptiveVertexFitter'),
+            chi2cutoff = cms.double(2.5),
+            label = cms.string('WithBS'),
+            maxDistanceToBeam = cms.double(1.0),
+            minNdof = cms.double(2.0),
+            useBeamConstraint = cms.bool(True),
+            vertexTimeParameters = cms.PSet(
+                algorithm = cms.string('legacy4D')
+            )
+        )
+    )
 )
 trackWithVertexRefSelectorBeforeSorting4D = trackWithVertexRefSelector.clone(
     vertexTag = "unsortedOfflinePrimaryVertices4D",
@@ -27,28 +51,6 @@
     vertices = "unsortedOfflinePrimaryVertices4D:WithBS"
 )
 
-unsortedOfflinePrimaryVertices4DnoPID = unsortedOfflinePrimaryVertices4D.clone(
-    TrackTimesLabel = "trackExtenderWithMTD:generalTrackt0",
-    TrackTimeResosLabel = "trackExtenderWithMTD:generalTracksigmat0"
-)
-trackWithVertexRefSelectorBeforeSorting4DnoPID = trackWithVertexRefSelector.clone(
-    vertexTag = "unsortedOfflinePrimaryVertices4DnoPID",
-    ptMax = 9e99,
-    ptErrorCut = 9e99
-)
-trackRefsForJetsBeforeSorting4DnoPID = trackRefsForJets.clone(
-    src = "trackWithVertexRefSelectorBeforeSorting4DnoPID"
-)
-offlinePrimaryVertices4DnoPID = offlinePrimaryVertices4D.clone(
-    vertices = "unsortedOfflinePrimaryVertices4DnoPID",
-    particles = "trackRefsForJetsBeforeSorting4DnoPID",
-    trackTimeTag = "trackExtenderWithMTD:generalTrackt0",
-    trackTimeResoTag = "trackExtenderWithMTD:generalTracksigmat0"
-)
-offlinePrimaryVertices4DnoPIDWithBS=offlinePrimaryVertices4DnoPID.clone(
-    vertices = "unsortedOfflinePrimaryVertices4DnoPID:WithBS"
-)
-
 unsortedOfflinePrimaryVertices4DwithPID = unsortedOfflinePrimaryVertices4D.clone(
     TrackTimesLabel = "tofPID4DnoPID:t0safe",
     TrackTimeResosLabel = "tofPID4DnoPID:sigmat0safe"
@@ -76,9 +78,37 @@
 from SimTracker.TrackAssociation.trackTimeValueMapProducer_cfi import trackTimeValueMapProducer
 from RecoMTD.TimingIDTools.tofPIDProducer_cfi import tofPIDProducer
 
-tofPID4DnoPID=tofPIDProducer.clone(vtxsSrc='unsortedOfflinePrimaryVertices4DnoPID')
+tofPID4DnoPID=tofPIDProducer.clone(vtxsSrc='unsortedOfflinePrimaryVertices')
 tofPID=tofPIDProducer.clone()
+tofPID3D=tofPIDProducer.clone(vtxsSrc='unsortedOfflinePrimaryVertices')
 
 from Configuration.Eras.Modifier_phase2_timing_layer_cff import phase2_timing_layer
-phase2_timing_layer.toModify(tofPID, vtxsSrc='unsortedOfflinePrimaryVertices4D')
+phase2_timing_layer.toModify(tofPID, vtxsSrc='unsortedOfflinePrimaryVertices4D', vertexReassignment=False)
+phase2_timing_layer.toModify(tofPID3D, vertexReassignment=False)
+phase2_timing_layer.toModify(unsortedOfflinePrimaryVertices, 
+    vertexCollections = cms.VPSet(
+        cms.PSet(
+            algorithm = cms.string('AdaptiveVertexFitter'),
+            chi2cutoff = cms.double(2.5),
+            label = cms.string(''),
+            maxDistanceToBeam = cms.double(1.0),
+            minNdof = cms.double(0.0),
+            useBeamConstraint = cms.bool(False),
+            vertexTimeParameters = cms.PSet(
+                algorithm = cms.string('fromTracksPID')
+            )
+        ),
+        cms.PSet(
+            algorithm = cms.string('AdaptiveVertexFitter'),
+            chi2cutoff = cms.double(2.5),
+            label = cms.string('WithBS'),
+            maxDistanceToBeam = cms.double(1.0),
+            minNdof = cms.double(2.0),
+            useBeamConstraint = cms.bool(True),
+            vertexTimeParameters = cms.PSet(
+                algorithm = cms.string('fromTracksPID')
+            )
+        )
+    )
+)
 

RecoVertex/PrimaryVertexProducer/interface/PrimaryVertexProducer.h

@@ -106,6 +106,11 @@ class PrimaryVertexProducer : public edm::stream::EDProducer<> {
   edm::EDGetTokenT<reco::TrackCollection> trkToken;
   edm::EDGetTokenT<edm::ValueMap<float> > trkTimesToken;
   edm::EDGetTokenT<edm::ValueMap<float> > trkTimeResosToken;
+  edm::EDGetTokenT<edm::ValueMap<float> > trackMTDTimeQualityToken;
 
   bool useTransientTrackTime;
+  bool useMVASelection;
+  edm::ValueMap<float> trackMTDTimeQualities_;
+  edm::ValueMap<float> trackTimes_;
+  double minTrackTimeQuality;
 };