Merge pull request #29771 from CTPPS/pps_re_miniaod_10_6

PPS: updates for re-miniAOD (backport of #29770)

Configuration/StandardSequences/python/Reconstruction_cff.py

@@ -90,8 +90,8 @@
 ###########################################
 _fastSim_localreco = localreco.copyAndExclude([
     castorreco,
-    totemRPLocalReconstruction,totemTimingLocalReconstruction,ctppsDiamondLocalReconstruction,
-      ctppsLocalTrackLiteProducer,ctppsPixelLocalReconstruction,ctppsProtons,
+    totemRPLocalReconstructionTask,totemTimingLocalReconstructionTask,ctppsDiamondLocalReconstructionTask,
+      ctppsLocalTrackLiteProducer,ctppsPixelLocalReconstructionTask,ctppsProtons,
     trackerlocalreco
 ])
 fastSim.toReplaceWith(localreco, _fastSim_localreco)

PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py

@@ -548,13 +548,12 @@ def miniAOD_customizeOutput(out):
 def miniAOD_customizeData(process):
     from PhysicsTools.PatAlgos.tools.coreTools import runOnData
     runOnData( process, outputModules = [] )
-    process.load("RecoCTPPS.TotemRPLocal.ctppsLocalTrackLiteProducer_cff")
-    process.load("RecoCTPPS.ProtonReconstruction.ctppsProtons_cff")
-    process.load("Geometry.VeryForwardGeometry.geometryRPFromDB_cfi")
+    process.load("RecoCTPPS.Configuration.recoCTPPS_cff")
     task = getPatAlgosToolsTask(process)
     from Configuration.Eras.Modifier_ctpps_2016_cff import ctpps_2016
-    ctpps_2016.toModify(task, func=lambda t: t.add(process.ctppsLocalTrackLiteProducer))
-    ctpps_2016.toModify(task, func=lambda t: t.add(process.ctppsProtons))
+    from Configuration.Eras.Modifier_run2_miniAOD_devel_cff import run2_miniAOD_devel
+    (ctpps_2016 & ~run2_miniAOD_devel).toModify(task, func=lambda t: t.add(process.ctppsLocalTrackLiteProducer, process.ctppsProtons))
+    (ctpps_2016 & run2_miniAOD_devel).toModify(task, func=lambda t: t.add(process.recoCTPPSTask))
 
 def miniAOD_customizeAllData(process):
     miniAOD_customizeCommon(process)

RecoCTPPS/Configuration/python/recoCTPPS_cff.py

@@ -11,11 +11,12 @@
 
 from Geometry.VeryForwardGeometry.geometryRPFromDB_cfi import *
 
-recoCTPPS = cms.Sequence(
-    totemRPLocalReconstruction *
-    ctppsDiamondLocalReconstruction *
-    totemTimingLocalReconstruction *
-    ctppsPixelLocalReconstruction *
-    ctppsLocalTrackLiteProducer *
+recoCTPPSTask = cms.Task(
+    totemRPLocalReconstructionTask ,
+    ctppsDiamondLocalReconstructionTask ,
+    totemTimingLocalReconstructionTask ,
+    ctppsPixelLocalReconstructionTask ,
+    ctppsLocalTrackLiteProducer ,
     ctppsProtons
 )
+recoCTPPS = cms.Sequence(recoCTPPSTask)

RecoCTPPS/PixelLocal/python/ctppsPixelLocalReconstruction_cff.py

@@ -9,6 +9,7 @@
 # local track producer
 from RecoCTPPS.PixelLocal.ctppsPixelLocalTracks_cfi import ctppsPixelLocalTracks
 
-ctppsPixelLocalReconstruction = cms.Sequence(
-    ctppsPixelClusters*ctppsPixelRecHits*ctppsPixelLocalTracks
+ctppsPixelLocalReconstructionTask = cms.Task(
+    ctppsPixelClusters,ctppsPixelRecHits,ctppsPixelLocalTracks
 )
+ctppsPixelLocalReconstruction = cms.Sequence(ctppsPixelLocalReconstructionTask)

RecoCTPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc

@@ -45,6 +45,8 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<>
 
     edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
 
+    bool pixelDiscardBXShiftedTracks_;
+
     std::string lhcInfoLabel_;
     std::string opticsLabel_;
 
@@ -61,26 +63,33 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<>
     struct AssociationCuts
     {
       bool x_cut_apply;
-      double x_cut_value;
+      double x_cut_mean, x_cut_value;
       bool y_cut_apply;
-      double y_cut_value;
+      double y_cut_mean, y_cut_value;
       bool xi_cut_apply;
-      double xi_cut_value;
+      double xi_cut_mean, xi_cut_value;
       bool th_y_cut_apply;
-      double th_y_cut_value;
+      double th_y_cut_mean, th_y_cut_value;
 
       double ti_tr_min;
       double ti_tr_max;
 
       void load(const edm::ParameterSet &ps)
       {
         x_cut_apply    = ps.getParameter<bool>  ("x_cut_apply");
+        x_cut_mean    = ps.getParameter<double>("x_cut_mean");
         x_cut_value    = ps.getParameter<double>("x_cut_value");
+
         y_cut_apply    = ps.getParameter<bool>  ("y_cut_apply");
+        y_cut_mean    = ps.getParameter<double>("y_cut_mean");
         y_cut_value    = ps.getParameter<double>("y_cut_value");
+
         xi_cut_apply   = ps.getParameter<bool>  ("xi_cut_apply");
+        xi_cut_mean   = ps.getParameter<double>("xi_cut_mean");
         xi_cut_value   = ps.getParameter<double>("xi_cut_value");
+
         th_y_cut_apply = ps.getParameter<bool>  ("th_y_cut_apply");
+        th_y_cut_mean = ps.getParameter<double>("th_y_cut_mean");
         th_y_cut_value = ps.getParameter<double>("th_y_cut_value");
 
         ti_tr_min = ps.getParameter<double>("ti_tr_min");
@@ -92,12 +101,19 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<>
         edm::ParameterSetDescription desc;
 
         desc.add<bool>("x_cut_apply", false)->setComment("whether to apply track-association cut in x");
+        desc.add<double>("x_cut_mean", 0E-6)->setComment("mean of track-association cut in x, mm");
         desc.add<double>("x_cut_value", 800E-6)->setComment("threshold of track-association cut in x, mm");
+
         desc.add<bool>("y_cut_apply", false)->setComment("whether to apply track-association cut in y");
+        desc.add<double>("y_cut_mean", 0E-6)->setComment("mean of track-association cut in y, mm");
         desc.add<double>("y_cut_value", 600E-6)->setComment("threshold of track-association cut in y, mm");
+
         desc.add<bool>("xi_cut_apply", true)->setComment("whether to apply track-association cut in xi");
+        desc.add<double>("xi_cut_mean", 0.)->setComment("mean of track-association cut in xi");
         desc.add<double>("xi_cut_value", 0.013)->setComment("threshold of track-association cut in xi");
+
         desc.add<bool>("th_y_cut_apply", true)->setComment("whether to apply track-association cut in th_y");
+        desc.add<double>("th_y_cut_mean", 0E-6)->setComment("mean of track-association cut in th_y, rad");
         desc.add<double>("th_y_cut_value", 20E-6)->setComment("threshold of track-association cut in th_y, rad");
 
         desc.add<double>("ti_tr_min", -1.)->setComment("minimum value for timing-tracking association cut");
@@ -110,6 +126,7 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<>
     std::map<unsigned int, AssociationCuts> association_cuts_;  // map: arm -> AssociationCuts
 
     unsigned int max_n_timing_tracks_;
+    double default_time_;
 
     ProtonReconstructionAlgorithm algorithm_;
 
@@ -121,6 +138,9 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<>
 
 CTPPSProtonProducer::CTPPSProtonProducer(const edm::ParameterSet& iConfig) :
   tracksToken_                (consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagLocalTrackLite"))),
+
+  pixelDiscardBXShiftedTracks_(iConfig.getParameter<bool>("pixelDiscardBXShiftedTracks")),
+
   lhcInfoLabel_               (iConfig.getParameter<std::string>("lhcInfoLabel")),
   opticsLabel_                (iConfig.getParameter<std::string>("opticsLabel")),
   verbosity_                  (iConfig.getUntrackedParameter<unsigned int>("verbosity", 0)),
@@ -135,6 +155,7 @@ CTPPSProtonProducer::CTPPSProtonProducer(const edm::ParameterSet& iConfig) :
   localAngleYMax_             (iConfig.getParameter<double>("localAngleYMax")),
 
   max_n_timing_tracks_        (iConfig.getParameter<unsigned int>("max_n_timing_tracks")),
+  default_time_               (iConfig.getParameter<double>("default_time")),
 
   algorithm_                  (iConfig.getParameter<bool>("fitVtxY"), iConfig.getParameter<bool>("useImprovedInitialEstimate"), verbosity_),
   opticsValid_(false)
@@ -161,6 +182,9 @@ void CTPPSProtonProducer::fillDescriptions(edm::ConfigurationDescriptions& descr
   desc.add<edm::InputTag>("tagLocalTrackLite", edm::InputTag("ctppsLocalTrackLiteProducer"))
     ->setComment("specification of the input lite-track collection");
 
+  desc.add<bool>("pixelDiscardBXShiftedTracks", false)
+      ->setComment("whether to discard pixel tracks built from BX-shifted planes");
+
   desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
   desc.add<std::string>("opticsLabel", "")->setComment("label of the optics record");
 
@@ -193,6 +217,8 @@ void CTPPSProtonProducer::fillDescriptions(edm::ConfigurationDescriptions& descr
 
   desc.add<unsigned int>("max_n_timing_tracks", 5)->setComment("maximum number of timing tracks per RP");
 
+  desc.add<double>("default_time", 0.)->setComment("proton time to be used when no timing information available");
+
   desc.add<bool>("fitVtxY", true)
     ->setComment("for multi-RP reconstruction, flag whether y* should be free fit parameter");
 
@@ -261,6 +287,13 @@ void CTPPSProtonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe
             || tr.getTy() < localAngleYMin_ || tr.getTy() > localAngleYMax_)
           continue;
 
+        if (pixelDiscardBXShiftedTracks_)
+        {
+          if (tr.getPixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes ||
+                tr.getPixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes)
+            continue;
+        }
+
         const CTPPSDetId rpId(tr.getRPId());
 
         if (verbosity_)
@@ -308,34 +341,37 @@ void CTPPSProtonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe
         // do multi-RP reco if chosen
         if (doMultiRPReconstruction_ && rpIds.size() == 2)
         {
-          // find matching track pairs from different tracking RPs
+          // find matching track pairs from different tracking RPs, ordered: i=near, j=far RP
           std::vector<std::pair<unsigned int, unsigned int>> idx_pairs;
           std::map<unsigned int, unsigned int> idx_pair_multiplicity;
           for (const auto &i : indices)
           {
             for (const auto &j : indices)
             {
-              if (j <= i)
-                continue;
-
               const auto &tr_i = hTracks->at(i);
               const auto &tr_j = hTracks->at(j);
 
+              const double z_i = hGeometry->getRPTranslation(tr_i.getRPId()).z();
+              const double z_j = hGeometry->getRPTranslation(tr_j.getRPId()).z();
+
               const auto &pr_i = singleRPResultsIndexed[i];
               const auto &pr_j = singleRPResultsIndexed[j];
 
               if (tr_i.getRPId() == tr_j.getRPId())
                 continue;
 
+              if (std::abs(z_i) >= std::abs(z_j))
+                continue;
+
               bool matching = true;
 
-              if (ac.x_cut_apply && std::abs(tr_i.getX() - tr_j.getX()) > ac.x_cut_value)
+              if (ac.x_cut_apply && std::abs(tr_i.getX() - tr_j.getX() - ac.x_cut_mean) > ac.x_cut_value)
                 matching = false;
-              if (ac.y_cut_apply && std::abs(tr_i.getY() - tr_j.getY()) > ac.y_cut_value)
+              else if (ac.y_cut_apply && std::abs(tr_i.getY() - tr_j.getY() - ac.y_cut_mean) > ac.y_cut_value)
                 matching = false;
-              if (ac.xi_cut_apply && std::abs(pr_i.xi() - pr_j.xi()) > ac.xi_cut_value)
+              else if (ac.xi_cut_apply && std::abs(pr_i.xi() - pr_j.xi() - ac.xi_cut_mean) > ac.xi_cut_value)
                 matching = false;
-              if (ac.th_y_cut_apply && std::abs(pr_i.thetaY() - pr_j.thetaY()) > ac.th_y_cut_value)
+              else if (ac.th_y_cut_apply && std::abs(pr_i.thetaY() - pr_j.thetaY() - ac.th_y_cut_mean) > ac.th_y_cut_value)
                 matching = false;
 
               if (!matching)
@@ -447,7 +483,7 @@ void CTPPSProtonProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSe
               swt += w * tr.getTime();
             }
 
-            float time = 0., time_unc = 0.;
+            float time = default_time_, time_unc = 0.;
             if (sw > 0.)
             {
               time = swt / sw;

RecoCTPPS/ProtonReconstruction/python/ctppsProtons_cff.py

@@ -14,6 +14,8 @@
 from Configuration.Eras.Modifier_ctpps_2017_cff import ctpps_2017
 from Configuration.Eras.Modifier_ctpps_2018_cff import ctpps_2018
 
+from Configuration.Eras.Modifier_run2_miniAOD_devel_cff import run2_miniAOD_devel
+
 def applyDefaultSettings(ctppsProtons):
   ctppsProtons.association_cuts_45.x_cut_apply = False
   ctppsProtons.association_cuts_45.y_cut_apply = False
@@ -27,12 +29,48 @@ def applyDefaultSettings(ctppsProtons):
   ctppsProtons.association_cuts_56.xi_cut_value = 0.015
   ctppsProtons.association_cuts_56.th_y_cut_apply = False
 
+  # update for re-miniAOD
+  run2_miniAOD_devel.toModify(ctppsProtons,
+    pixelDiscardBXShiftedTracks = True,
+    association_cuts_45 = dict(ti_tr_min = -1.5, ti_tr_max = 2.0),
+    association_cuts_56 = dict(ti_tr_min = -1.5, ti_tr_max = 2.0),
+    default_time = -999.
+  )
+
 ctpps_2016.toModify(ctppsProtons, applyDefaultSettings) # applied for all Run2 years (2016, 2017 and 2018)
 
 def apply2017Settings(ctppsProtons):
   ctppsProtons.association_cuts_45.xi_cut_value = 0.010
   ctppsProtons.association_cuts_56.xi_cut_value = 0.015
 
+  # update for re-miniAOD
+  run2_miniAOD_devel.toModify(ctppsProtons,
+    association_cuts_45 = dict(
+      x_cut_apply = False,
+      y_cut_apply = False,
+
+      xi_cut_apply = True,
+      xi_cut_value = 5. * 0.00121,
+      xi_cut_mean = +6.0695e-5,
+
+      th_y_cut_apply = False
+    ),
+
+    association_cuts_56 = dict(
+      x_cut_apply = False,
+
+      y_cut_apply = True,
+      y_cut_value = 5. * 0.14777,
+      y_cut_mean = -0.022612,
+
+      xi_cut_apply = True,
+      xi_cut_value = 5. * 0.0020627,
+      xi_cut_mean = +8.012857e-5,
+
+      th_y_cut_apply = False
+    )
+  )
+
 ctpps_2017.toModify(ctppsProtons, apply2017Settings)
 
 def apply2018Settings(ctppsProtons):
@@ -44,4 +82,39 @@ def apply2018Settings(ctppsProtons):
   ctppsProtons.association_cuts_56.th_y_cut_apply = True
   ctppsProtons.association_cuts_56.th_y_cut_value = 20E-6
 
+  # update for re-miniAOD
+  run2_miniAOD_devel.toModify(ctppsProtons,
+    association_cuts_45 = dict(
+      x_cut_apply = True,
+      x_cut_value = 4. * 0.16008188,
+      x_cut_mean = -0.065194856,
+
+      y_cut_apply = True,
+      y_cut_value = 4. * 0.1407986,
+      y_cut_mean = +0.10973631,
+
+      xi_cut_apply = True,
+      xi_cut_value = 4. * 0.0012403586,
+      xi_cut_mean = +3.113062e-5,
+
+      th_y_cut_apply = False
+    ),
+
+    association_cuts_56 = dict(
+      x_cut_apply = True,
+      x_cut_value = 5. * 0.18126434,
+      x_cut_mean = +0.073016431,
+
+      y_cut_apply = True,
+      y_cut_value = 5. * 0.14990802,
+      y_cut_mean = +0.064261029,
+
+      xi_cut_apply = True,
+      xi_cut_value = 5. * 0.002046409,
+      xi_cut_mean = -1.1852528e-5,
+
+      th_y_cut_apply = False
+    )
+  )
+
 ctpps_2018.toModify(ctppsProtons, apply2018Settings)

RecoCTPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc

@@ -251,6 +251,8 @@ reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromMultiRP(const
 
   if (verbosity_)
     os
+      << "    fit valid=" << result.IsValid()
+      << std::endl
       << "    xi=" << params[0] << " +- " << result.Error(0)
       << ", th_x=" << params[1] << " +-" << result.Error(1)
       << ", th_y=" << params[2] << " +-" << result.Error(2)

RecoCTPPS/TotemRPLocal/python/ctppsDiamondLocalReconstruction_cff.py

@@ -6,7 +6,8 @@
 # local track fitting
 from RecoCTPPS.TotemRPLocal.ctppsDiamondLocalTracks_cfi import ctppsDiamondLocalTracks
 
-ctppsDiamondLocalReconstruction = cms.Sequence(
-    ctppsDiamondRecHits *
+ctppsDiamondLocalReconstructionTask = cms.Task(
+    ctppsDiamondRecHits,
     ctppsDiamondLocalTracks
 )
+ctppsDiamondLocalReconstruction = cms.Sequence(ctppsDiamondLocalReconstructionTask)

RecoCTPPS/TotemRPLocal/python/totemRPLocalReconstruction_cff.py

@@ -12,9 +12,10 @@
 # local track fitting
 from RecoCTPPS.TotemRPLocal.totemRPLocalTrackFitter_cfi import *
 
-totemRPLocalReconstruction = cms.Sequence(
-    totemRPClusterProducer *
-    totemRPRecHitProducer *
-    totemRPUVPatternFinder *
+totemRPLocalReconstructionTask = cms.Task(
+    totemRPClusterProducer ,
+    totemRPRecHitProducer ,
+    totemRPUVPatternFinder ,
     totemRPLocalTrackFitter
 )
+totemRPLocalReconstruction = cms.Sequence(totemRPLocalReconstructionTask)

RecoCTPPS/TotemRPLocal/python/totemTimingLocalReconstruction_cff.py

@@ -6,7 +6,8 @@
 # local track fitting
 from RecoCTPPS.TotemRPLocal.totemTimingLocalTracks_cfi import totemTimingLocalTracks
 
-totemTimingLocalReconstruction = cms.Sequence(
-    totemTimingRecHits *
+totemTimingLocalReconstructionTask = cms.Task(
+    totemTimingRecHits ,
     totemTimingLocalTracks
 )
+totemTimingLocalReconstruction = cms.Sequence(totemTimingLocalReconstructionTask)

Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc

@@ -192,7 +192,7 @@ void CTPPSDirectProtonSimulation::produce(edm::Event &iEvent, const edm::EventSe
   std::unique_ptr<std::vector<CTPPSLocalTrackLite>> pTracks(new std::vector<CTPPSLocalTrackLite>());
 
   // loop over event vertices
-  auto evt = new HepMC::GenEvent(*hepmc_prod->GetEvent());
+  auto evt = hepmc_prod->GetEvent();
   for (auto it_vtx = evt->vertices_begin(); it_vtx != evt->vertices_end(); ++it_vtx) {
     auto vtx = *(it_vtx);