diff --git a/Configuration/StandardSequences/python/Reconstruction_cff.py b/Configuration/StandardSequences/python/Reconstruction_cff.py
index 13f475ceb8b7f..11dfffe90809e 100644
--- a/Configuration/StandardSequences/python/Reconstruction_cff.py
+++ b/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)
diff --git a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py b/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
index 5dcd650d031b4..2c837a10c0e0d 100644
--- a/PhysicsTools/PatAlgos/python/slimming/miniAOD_tools.py
+++ b/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)
diff --git a/RecoCTPPS/Configuration/python/recoCTPPS_cff.py b/RecoCTPPS/Configuration/python/recoCTPPS_cff.py
index 735a13129458d..c3fdf7f4563c8 100644
--- a/RecoCTPPS/Configuration/python/recoCTPPS_cff.py
+++ b/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)
diff --git a/RecoCTPPS/PixelLocal/python/ctppsPixelLocalReconstruction_cff.py b/RecoCTPPS/PixelLocal/python/ctppsPixelLocalReconstruction_cff.py
index 5a8dc459605fc..348e98db107aa 100644
--- a/RecoCTPPS/PixelLocal/python/ctppsPixelLocalReconstruction_cff.py
+++ b/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)
diff --git a/RecoCTPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc b/RecoCTPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc
index 9579ae1321dca..fbe372115771e 100644
--- a/RecoCTPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc
+++ b/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,13 +63,13 @@ 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;
@@ -75,12 +77,19 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<>
       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;
diff --git a/RecoCTPPS/ProtonReconstruction/python/ctppsProtons_cff.py b/RecoCTPPS/ProtonReconstruction/python/ctppsProtons_cff.py
index c67cd87471496..1016061033dfc 100644
--- a/RecoCTPPS/ProtonReconstruction/python/ctppsProtons_cff.py
+++ b/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)
diff --git a/RecoCTPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc b/RecoCTPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc
index d6c4dd0e3d590..95aec35565e0f 100644
--- a/RecoCTPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc
+++ b/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)
diff --git a/RecoCTPPS/TotemRPLocal/python/ctppsDiamondLocalReconstruction_cff.py b/RecoCTPPS/TotemRPLocal/python/ctppsDiamondLocalReconstruction_cff.py
index 1243e3a0fd2aa..b17cd18f76260 100644
--- a/RecoCTPPS/TotemRPLocal/python/ctppsDiamondLocalReconstruction_cff.py
+++ b/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)
diff --git a/RecoCTPPS/TotemRPLocal/python/totemRPLocalReconstruction_cff.py b/RecoCTPPS/TotemRPLocal/python/totemRPLocalReconstruction_cff.py
index 854f8b6747fad..f8491569c951b 100644
--- a/RecoCTPPS/TotemRPLocal/python/totemRPLocalReconstruction_cff.py
+++ b/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)
diff --git a/RecoCTPPS/TotemRPLocal/python/totemTimingLocalReconstruction_cff.py b/RecoCTPPS/TotemRPLocal/python/totemTimingLocalReconstruction_cff.py
index 99d8051cf6668..4536fa219c67f 100644
--- a/RecoCTPPS/TotemRPLocal/python/totemTimingLocalReconstruction_cff.py
+++ b/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)
diff --git a/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc b/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc
index 3862ffce8a33e..c4aca9b419588 100644
--- a/Validation/CTPPS/plugins/CTPPSDirectProtonSimulation.cc
+++ b/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);
 
diff --git a/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc b/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc
index 99379d6564e0e..977b47235c240 100644
--- a/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc
+++ b/Validation/CTPPS/plugins/CTPPSLHCInfoPlotter.cc
@@ -34,6 +34,9 @@ class CTPPSLHCInfoPlotter : public edm::one::EDAnalyzer<> {
   TH1D *h_beamEnergy_;
   TH1D *h_xangle_;
   TH1D *h_betaStar_;
+
+  TH1D *h_fill_;
+  TH1D *h_run_;
 };
 
 //----------------------------------------------------------------------------------------------------
@@ -49,7 +52,10 @@ CTPPSLHCInfoPlotter::CTPPSLHCInfoPlotter(const edm::ParameterSet &iConfig)
 
       h_beamEnergy_(new TH1D("h_beamEnergy", ";beam energy   (GeV)", 81, -50., 8050.)),
       h_xangle_(new TH1D("h_xangle", ";(half) crossing angle   (#murad)", 201, -0.5, 200.5)),
-      h_betaStar_(new TH1D("h_betaStar", ";#beta^{*}   (m)", 101, -0.005, 1.005)) {}
+      h_betaStar_(new TH1D("h_betaStar", ";#beta^{*}   (m)", 101, -0.005, 1.005)),
+
+      h_fill_(new TH1D("h_fill", ";fill", 4001, 3999.5, 8000.5)),
+      h_run_(new TH1D("h_run", ";run", 6000, 270E3, 330E3)) {}
 
 //----------------------------------------------------------------------------------------------------
 
@@ -60,6 +66,9 @@ void CTPPSLHCInfoPlotter::analyze(const edm::Event &iEvent, const edm::EventSetu
   h_beamEnergy_->Fill(hLHCInfo->energy());
   h_xangle_->Fill(hLHCInfo->crossingAngle());
   h_betaStar_->Fill(hLHCInfo->betaStar());
+
+  h_fill_->Fill(hLHCInfo->fillNumber());
+  h_run_->Fill(iEvent.id().run());
 }
 
 //----------------------------------------------------------------------------------------------------
@@ -70,6 +79,9 @@ void CTPPSLHCInfoPlotter::endJob() {
   h_beamEnergy_->Write();
   h_xangle_->Write();
   h_betaStar_->Write();
+
+  h_fill_->Write();
+  h_run_->Write();
 }
 
 //----------------------------------------------------------------------------------------------------
diff --git a/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc b/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc
index 0cb26773a3968..189a20eda6604 100644
--- a/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc
+++ b/Validation/CTPPS/plugins/CTPPSOpticsPlotter.cc
@@ -34,6 +34,9 @@ class CTPPSOpticsPlotter : public edm::one::EDAnalyzer<> {
 
   std::string opticsLabel_;
 
+  unsigned int rpId_45_N_, rpId_45_F_;
+  unsigned int rpId_56_N_, rpId_56_F_;
+
   std::string outputFile_;
 
   struct RPPlots {
@@ -75,6 +78,24 @@ class CTPPSOpticsPlotter : public edm::one::EDAnalyzer<> {
   };
 
   std::map<unsigned int, RPPlots> rp_plots_;
+
+  struct ArmPlots {
+    unsigned int id_N, id_F;
+
+    std::unique_ptr<TGraph> g_de_x_vs_x_disp, g_de_y_vs_x_disp;
+
+    ArmPlots() : g_de_x_vs_x_disp(new TGraph), g_de_y_vs_x_disp(new TGraph) {}
+
+    void write() const {
+      g_de_x_vs_x_disp->SetTitle(";x_N   (cm);x_F - x_N   (cm)");
+      g_de_x_vs_x_disp->Write("g_de_x_vs_x_disp");
+
+      g_de_y_vs_x_disp->SetTitle(";x_N   (cm);y_F - y_N   (cm)");
+      g_de_y_vs_x_disp->Write("g_de_y_vs_x_disp");
+    }
+  };
+
+  std::map<unsigned int, ArmPlots> arm_plots_;
 };
 
 //----------------------------------------------------------------------------------------------------
@@ -86,7 +107,19 @@ using namespace edm;
 
 CTPPSOpticsPlotter::CTPPSOpticsPlotter(const edm::ParameterSet& iConfig)
     : opticsLabel_(iConfig.getParameter<std::string>("opticsLabel")),
-      outputFile_(iConfig.getParameter<string>("outputFile")) {}
+
+      rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
+      rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
+      rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
+      rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
+
+      outputFile_(iConfig.getParameter<string>("outputFile")) {
+  arm_plots_[0].id_N = rpId_45_N_;
+  arm_plots_[0].id_F = rpId_45_F_;
+
+  arm_plots_[1].id_N = rpId_56_N_;
+  arm_plots_[1].id_F = rpId_56_F_;
+}
 
 //----------------------------------------------------------------------------------------------------
 
@@ -103,7 +136,7 @@ void CTPPSOpticsPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup
   if (hOpticalFunctions->empty())
     return;
 
-  // make plots
+  // make per-RP plots
   for (const auto& it : *hOpticalFunctions) {
     CTPPSDetId rpId(it.first);
     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
@@ -151,6 +184,48 @@ void CTPPSOpticsPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup
       pl.h_y_vs_x_disp->SetPoint(idx, k_out_xi.x - k_out_beam.x, k_out_xi.y - k_out_beam.y);
     }
   }
+
+  // make per-arm plots
+  for (const auto& ap : arm_plots_) {
+    // find optics objects
+    const LHCInterpolatedOpticalFunctionsSet *opt_N = NULL, *opt_F = NULL;
+
+    for (const auto& it : *hOpticalFunctions) {
+      CTPPSDetId rpId(it.first);
+      unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
+
+      if (rpDecId == ap.second.id_N)
+        opt_N = &it.second;
+      if (rpDecId == ap.second.id_F)
+        opt_F = &it.second;
+    }
+
+    if (!opt_N || !opt_F) {
+      edm::LogError("CTPPSOpticsPlotter::analyze") << "Cannot find optics objects for arm " << ap.first;
+      continue;
+    }
+
+    LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
+
+    LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam_N, k_out_beam_F;
+    opt_N->transport(k_in_beam, k_out_beam_N);
+    opt_F->transport(k_in_beam, k_out_beam_F);
+
+    for (double xi = 0.; xi < 0.30001; xi += 0.001) {
+      LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
+
+      LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_N, k_out_xi_F;
+      opt_N->transport(k_in_xi, k_out_xi_N);
+      opt_F->transport(k_in_xi, k_out_xi_F);
+
+      int idx = ap.second.g_de_x_vs_x_disp->GetN();
+
+      ap.second.g_de_x_vs_x_disp->SetPoint(
+          idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.x - k_out_beam_F.x) - (k_out_xi_N.x - k_out_beam_N.x));
+      ap.second.g_de_y_vs_x_disp->SetPoint(
+          idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.y - k_out_beam_F.y) - (k_out_xi_N.y - k_out_beam_N.y));
+    }
+  }
 }
 
 //----------------------------------------------------------------------------------------------------
@@ -162,6 +237,11 @@ void CTPPSOpticsPlotter::endJob() {
     gDirectory = f_out->mkdir(Form("%u", p.first));
     p.second.write();
   }
+
+  for (const auto& p : arm_plots_) {
+    gDirectory = f_out->mkdir(Form("arm %u", p.first));
+    p.second.write();
+  }
 }
 
 //----------------------------------------------------------------------------------------------------
diff --git a/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorData.cc b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorData.cc
new file mode 100644
index 0000000000000..f423035d129b5
--- /dev/null
+++ b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorData.cc
@@ -0,0 +1,681 @@
+/****************************************************************************
+ * Authors:
+ *   Jan Kašpar
+ *   Mauricio Thiel
+ ****************************************************************************/
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/one/EDAnalyzer.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/ESWatcher.h"
+
+#include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
+#include "CondFormats/CTPPSReadoutObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
+
+#include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
+
+#include "DataFormats/ProtonReco/interface/ForwardProton.h"
+#include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
+#include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
+#include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
+
+#include "TFile.h"
+#include "TH1D.h"
+#include "TH2D.h"
+#include "TProfile.h"
+#include "TF1.h"
+#include "TGraph.h"
+
+#include <map>
+#include <string>
+#include <sstream>
+
+//----------------------------------------------------------------------------------------------------
+
+class CTPPSProtonReconstructionEfficiencyEstimatorData : public edm::one::EDAnalyzer<> {
+public:
+  explicit CTPPSProtonReconstructionEfficiencyEstimatorData(const edm::ParameterSet &);
+
+  static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
+
+private:
+  void analyze(const edm::Event &, const edm::EventSetup &) override;
+  void endJob() override;
+
+  edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tokenTracks_;
+
+  edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsMultiRP_;
+
+  bool pixelDiscardBXShiftedTracks_;
+
+  double localAngleXMin_, localAngleXMax_, localAngleYMin_, localAngleYMax_;
+
+  std::string opticsLabel_;
+
+  unsigned int n_prep_events_;
+
+  unsigned int n_exp_prot_max_;
+
+  std::vector<double> n_sigmas_;
+
+  std::string outputFile_;
+
+  unsigned int verbosity_;
+
+  edm::ESWatcher<CTPPSInterpolatedOpticsRcd> opticsWatcher_;
+
+  struct ArmData {
+    unsigned int rpId_N, rpId_F;
+
+    std::shared_ptr<const TSpline3> s_x_to_xi_N;
+
+    unsigned int n_events_with_tracks;
+
+    std::unique_ptr<TH1D> h_de_x, h_de_y;
+    std::unique_ptr<TH2D> h2_de_y_vs_de_x;
+
+    double de_x_mean, de_x_sigma;
+    double de_y_mean, de_y_sigma;
+
+    std::vector<double> n_sigmas;
+
+    unsigned int n_exp_prot_max;
+
+    struct EffPlots {
+      std::unique_ptr<TProfile> p_eff1_vs_x_N;
+      std::unique_ptr<TProfile> p_eff1_vs_xi_N;
+
+      std::unique_ptr<TProfile> p_eff2_vs_x_N;
+      std::unique_ptr<TProfile> p_eff2_vs_xi_N;
+
+      EffPlots()
+          : p_eff1_vs_x_N(new TProfile("", ";x_{N}   (mm);efficiency", 50, 0., 25.)),
+            p_eff1_vs_xi_N(new TProfile("", ";#xi_{si,N};efficiency", 50, 0., 0.25)),
+
+            p_eff2_vs_x_N(new TProfile("", ";x_{N}   (mm);efficiency", 50, 0., 25.)),
+            p_eff2_vs_xi_N(new TProfile("", ";#xi_{si,N};efficiency", 50, 0., 0.25)) {}
+
+      void write() const {
+        p_eff1_vs_x_N->Write("p_eff1_vs_x_N");
+        p_eff1_vs_xi_N->Write("p_eff1_vs_xi_N");
+
+        p_eff2_vs_x_N->Write("p_eff2_vs_x_N");
+        p_eff2_vs_xi_N->Write("p_eff2_vs_xi_N");
+      }
+    };
+
+    // (n exp protons, index in n_sigmas) --> plots
+    // n exp protons = 0 --> no condition (summary case)
+    std::map<unsigned int, std::map<unsigned int, EffPlots>> effPlots;
+
+    ArmData()
+        : n_events_with_tracks(0),
+
+          h_de_x(new TH1D("", ";x_{F} - x_{N}   (mm)", 200, -1., +1.)),
+          h_de_y(new TH1D("", ";y_{F} - y_{N}   (mm)", 200, -1., +1.)),
+          h2_de_y_vs_de_x(new TH2D("", ";x_{F} - x_{N}   (mm);y_{F} - y_{N}   (mm)", 100, -1., +1., 100, -1., +1.)),
+
+          de_x_mean(0.),
+          de_x_sigma(0.),
+          de_y_mean(0.),
+          de_y_sigma(0.) {
+      for (unsigned int nep = 0; nep <= n_exp_prot_max; ++nep) {
+        for (unsigned int nsi = 0; nsi < n_sigmas.size(); ++nsi) {
+          effPlots[nep][nsi] = EffPlots();
+        }
+      }
+    }
+
+    void write() const {
+      h_de_x->Write("h_de_x");
+      h_de_y->Write("h_de_y");
+      h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x");
+
+      char buf[100];
+
+      for (const auto &n_si : n_sigmas) {
+        sprintf(buf, "g_de_y_vs_de_x_n_si=%.1f", n_si);
+        TGraph *g = new TGraph();
+        g->SetName(buf);
+
+        g->SetPoint(0, de_x_mean - n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma);
+        g->SetPoint(1, de_x_mean + n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma);
+        g->SetPoint(2, de_x_mean + n_si * de_x_sigma, de_y_mean + n_si * de_y_sigma);
+        g->SetPoint(3, de_x_mean - n_si * de_x_sigma, de_y_mean + n_si * de_y_sigma);
+        g->SetPoint(4, de_x_mean - n_si * de_x_sigma, de_y_mean - n_si * de_y_sigma);
+
+        g->Write();
+      }
+
+      TDirectory *d_top = gDirectory;
+
+      for (const auto &nep_p : effPlots) {
+        if (nep_p.first == 0)
+          sprintf(buf, "exp prot all");
+        else
+          sprintf(buf, "exp prot %u", nep_p.first);
+
+        TDirectory *d_nep = d_top->mkdir(buf);
+
+        for (const auto &nsi_p : nep_p.second) {
+          sprintf(buf, "nsi = %.1f", n_sigmas[nsi_p.first]);
+
+          TDirectory *d_nsi = d_nep->mkdir(buf);
+
+          gDirectory = d_nsi;
+
+          nsi_p.second.write();
+        }
+      }
+
+      gDirectory = d_top;
+    }
+
+    void UpdateOptics(const LHCInterpolatedOpticalFunctionsSetCollection &ofc) {
+      const LHCInterpolatedOpticalFunctionsSet *of = NULL;
+
+      for (const auto &p : ofc) {
+        CTPPSDetId rpId(p.first);
+        unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
+
+        if (decRPId == rpId_N) {
+          of = &p.second;
+          break;
+        }
+      }
+
+      if (!of) {
+        edm::LogError("ArmData::UpdateOptics") << "Cannot find optical functions for RP " << rpId_N;
+        return;
+      }
+
+      std::vector<double> xiValues =
+          of->getXiValues();  // local copy made since the TSpline constructor needs non-const parameters
+      std::vector<double> xDValues = of->getFcnValues()[LHCOpticalFunctionsSet::exd];
+      s_x_to_xi_N = std::make_shared<TSpline3>("", xDValues.data(), xiValues.data(), xiValues.size());
+    }
+  };
+
+  std::map<unsigned int, ArmData> data_;
+
+  std::unique_ptr<TF1> ff_;
+};
+
+//----------------------------------------------------------------------------------------------------
+
+using namespace std;
+using namespace edm;
+
+//----------------------------------------------------------------------------------------------------
+
+CTPPSProtonReconstructionEfficiencyEstimatorData::CTPPSProtonReconstructionEfficiencyEstimatorData(
+    const edm::ParameterSet &iConfig)
+    : tokenTracks_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
+
+      tokenRecoProtonsMultiRP_(
+          consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
+
+      pixelDiscardBXShiftedTracks_(iConfig.getParameter<bool>("pixelDiscardBXShiftedTracks")),
+
+      localAngleXMin_(iConfig.getParameter<double>("localAngleXMin")),
+      localAngleXMax_(iConfig.getParameter<double>("localAngleXMax")),
+      localAngleYMin_(iConfig.getParameter<double>("localAngleYMin")),
+      localAngleYMax_(iConfig.getParameter<double>("localAngleYMax")),
+
+      opticsLabel_(iConfig.getParameter<std::string>("opticsLabel")),
+      n_prep_events_(iConfig.getParameter<unsigned int>("n_prep_events")),
+      n_exp_prot_max_(iConfig.getParameter<unsigned int>("n_exp_prot_max")),
+      n_sigmas_(iConfig.getParameter<std::vector<double>>("n_sigmas")),
+
+      outputFile_(iConfig.getParameter<string>("outputFile")),
+
+      verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity")),
+
+      ff_(new TF1("ff", "[0] * exp(- (x-[1])*(x-[1]) / 2 / [2]/[2]) + [4]")) {
+  data_[0].n_exp_prot_max = n_exp_prot_max_;
+  data_[0].n_sigmas = n_sigmas_;
+  data_[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
+  data_[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
+
+  data_[1].n_exp_prot_max = n_exp_prot_max_;
+  data_[1].n_sigmas = n_sigmas_;
+  data_[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
+  data_[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
+}
+
+//----------------------------------------------------------------------------------------------------
+
+void CTPPSProtonReconstructionEfficiencyEstimatorData::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
+  edm::ParameterSetDescription desc;
+
+  desc.add<edm::InputTag>("tagTracks", edm::InputTag())->setComment("input tag for local lite tracks");
+  desc.add<edm::InputTag>("tagRecoProtonsMultiRP", edm::InputTag())->setComment("input tag for multi-RP reco protons");
+
+  desc.add<bool>("pixelDiscardBXShiftedTracks", false)
+      ->setComment("whether to discard pixel tracks built from BX-shifted planes");
+
+  desc.add<double>("localAngleXMin", -0.03)->setComment("minimal accepted value of local horizontal angle (rad)");
+  desc.add<double>("localAngleXMax", +0.03)->setComment("maximal accepted value of local horizontal angle (rad)");
+  desc.add<double>("localAngleYMin", -0.04)->setComment("minimal accepted value of local vertical angle (rad)");
+  desc.add<double>("localAngleYMax", +0.04)->setComment("maximal accepted value of local vertical angle (rad)");
+
+  desc.add<std::string>("opticsLabel", "")->setComment("label of the optics records");
+
+  desc.add<unsigned int>("n_prep_events", 1000)
+      ->setComment("number of preparatory events (to determine de x and de y window)");
+
+  desc.add<unsigned int>("n_exp_prot_max", 5)->setComment("maximum number of expected protons per event and per arm");
+
+  desc.add<std::vector<double>>("n_sigmas", {3., 5., 7.})->setComment("list of n_sigma values");
+
+  desc.add<unsigned int>("rpId_45_N", 0)->setComment("decimal RP id for 45 near");
+  desc.add<unsigned int>("rpId_45_F", 0)->setComment("decimal RP id for 45 far");
+  desc.add<unsigned int>("rpId_56_N", 0)->setComment("decimal RP id for 56 near");
+  desc.add<unsigned int>("rpId_56_F", 0)->setComment("decimal RP id for 56 far");
+
+  desc.add<std::string>("outputFile", "output.root")->setComment("output file name");
+
+  desc.addUntracked<unsigned int>("verbosity", 0)->setComment("verbosity level");
+
+  descriptions.add("ctppsProtonReconstructionEfficiencyEstimatorData", desc);
+}
+
+//----------------------------------------------------------------------------------------------------
+
+void CTPPSProtonReconstructionEfficiencyEstimatorData::analyze(const edm::Event &iEvent,
+                                                               const edm::EventSetup &iSetup) {
+  std::ostringstream os;
+
+  // get conditions
+  edm::ESHandle<LHCInterpolatedOpticalFunctionsSetCollection> hOpticalFunctions;
+  iSetup.get<CTPPSInterpolatedOpticsRcd>().get(opticsLabel_, hOpticalFunctions);
+
+  // check optics change
+  if (opticsWatcher_.check(iSetup)) {
+    data_[0].UpdateOptics(*hOpticalFunctions);
+    data_[1].UpdateOptics(*hOpticalFunctions);
+  }
+
+  // get input
+  edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
+  iEvent.getByToken(tokenTracks_, hTracks);
+
+  Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
+  iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
+
+  // pre-selection of tracks
+  std::vector<unsigned int> sel_track_idc;
+  for (unsigned int idx = 0; idx < hTracks->size(); ++idx) {
+    const auto &tr = hTracks->at(idx);
+
+    if (tr.getTx() < localAngleXMin_ || tr.getTx() > localAngleXMax_ || tr.getTy() < localAngleYMin_ ||
+        tr.getTy() > localAngleYMax_)
+      continue;
+
+    if (pixelDiscardBXShiftedTracks_) {
+      if (tr.getPixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes ||
+          tr.getPixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes)
+        continue;
+    }
+
+    sel_track_idc.push_back(idx);
+  }
+
+  // debug print
+  if (verbosity_ > 1) {
+    os << "* tracks (pre-selected):" << std::endl;
+
+    for (const auto idx : sel_track_idc) {
+      const auto &tr = hTracks->at(idx);
+      CTPPSDetId rpId(tr.getRPId());
+      unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
+
+      os << "    [" << idx << "] RP=" << decRPId << ", x=" << tr.getX() << ", y=" << tr.getY() << std::endl;
+    }
+
+    os << "* protons:" << std::endl;
+    for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) {
+      const auto &pr = (*hRecoProtonsMultiRP)[i];
+      os << "    [" << i << "] track indices: ";
+      for (const auto &trr : pr.contributingLocalTracks())
+        os << trr.key() << ", ";
+      os << std::endl;
+    }
+  }
+
+  // make de_x and de_y plots
+  map<unsigned int, bool> hasTracksInArm;
+
+  for (const auto idx_i : sel_track_idc) {
+    const auto &tr_i = hTracks->at(idx_i);
+    CTPPSDetId rpId_i(tr_i.getRPId());
+    unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
+
+    for (const auto idx_j : sel_track_idc) {
+      const auto &tr_j = hTracks->at(idx_j);
+      CTPPSDetId rpId_j(tr_j.getRPId());
+      unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
+
+      // check whether desired RP combination
+      unsigned int arm = 123;
+      for (const auto &ap : data_) {
+        if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
+          arm = ap.first;
+      }
+
+      if (arm > 1)
+        continue;
+
+      // update flag
+      hasTracksInArm[arm] = true;
+
+      // update plots
+      auto &ad = data_[arm];
+      const double de_x = tr_j.getX() - tr_i.getX();
+      const double de_y = tr_j.getY() - tr_i.getY();
+
+      if (ad.n_events_with_tracks < n_prep_events_) {
+        ad.h_de_x->Fill(de_x);
+        ad.h_de_y->Fill(de_y);
+      }
+
+      ad.h2_de_y_vs_de_x->Fill(de_x, de_y);
+    }
+  }
+
+  // update event counter
+  for (auto &ap : data_) {
+    if (hasTracksInArm[ap.first])
+      ap.second.n_events_with_tracks++;
+  }
+
+  // if threshold reached do fits
+  for (auto &ap : data_) {
+    auto &ad = ap.second;
+
+    if (ad.n_events_with_tracks == n_prep_events_) {
+      if (ad.de_x_sigma > 0. && ad.de_y_sigma > 0.)
+        continue;
+
+      std::vector<std::pair<unsigned int, TH1D *>> m;
+      m.emplace_back(0, ad.h_de_x.get());
+      m.emplace_back(1, ad.h_de_y.get());
+
+      for (const auto &p : m) {
+        double max_pos = -1E100, max_val = -1.;
+        for (int bi = 1; bi < p.second->GetNbinsX(); ++bi) {
+          const double pos = p.second->GetBinCenter(bi);
+          const double val = p.second->GetBinContent(bi);
+
+          if (val > max_val) {
+            max_val = val;
+            max_pos = pos;
+          }
+        }
+
+        const double sig = 0.2;
+
+        ff_->SetParameters(max_val, max_pos, sig, 0.);
+        p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig);
+        p.second->Fit(ff_.get(), "Q", "", max_pos - 3. * sig, max_pos + 3. * sig);
+
+        if (p.first == 0) {
+          ad.de_x_mean = ff_->GetParameter(1);
+          ad.de_x_sigma = fabs(ff_->GetParameter(2));
+        }
+        if (p.first == 1) {
+          ad.de_y_mean = ff_->GetParameter(1);
+          ad.de_y_sigma = fabs(ff_->GetParameter(2));
+        }
+      }
+
+      if (verbosity_) {
+        os << "* fitting arm " << ap.first << std::endl
+           << "    de_x: mean = " << ad.de_x_mean << ", sigma = " << ad.de_x_sigma << std::endl
+           << "    de_y: mean = " << ad.de_y_mean << ", sigma = " << ad.de_y_sigma;
+      }
+    }
+  }
+
+  // data structures for efficiency analysis
+  struct ArmEventData {
+    std::map<unsigned int, std::set<unsigned int>> matched_track_idc;
+
+    std::set<unsigned int> reco_proton_idc;
+
+    std::map<unsigned int, std::set<unsigned int>> matched_track_with_prot_idc, matched_track_without_prot_idc;
+  };
+
+  std::map<unsigned int, ArmEventData> armEventData;
+
+  // determine the number of expected protons
+  for (const auto idx_i : sel_track_idc) {
+    const auto &tr_i = hTracks->at(idx_i);
+    CTPPSDetId rpId_i(tr_i.getRPId());
+    unsigned int decRPId_i = rpId_i.arm() * 100 + rpId_i.station() * 10 + rpId_i.rp();
+
+    for (const auto idx_j : sel_track_idc) {
+      const auto &tr_j = hTracks->at(idx_j);
+      CTPPSDetId rpId_j(tr_j.getRPId());
+      unsigned int decRPId_j = rpId_j.arm() * 100 + rpId_j.station() * 10 + rpId_j.rp();
+
+      // check whether desired RP combination
+      unsigned int arm = 123;
+      for (const auto &ap : data_) {
+        if (ap.second.rpId_N == decRPId_i && ap.second.rpId_F == decRPId_j)
+          arm = ap.first;
+      }
+
+      if (arm > 1)
+        continue;
+
+      // check near-far matching
+      auto &ad = data_[arm];
+      const double de_x = tr_j.getX() - tr_i.getX();
+      const double de_y = tr_j.getY() - tr_i.getY();
+
+      for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
+        const double n_si = ad.n_sigmas[nsi];
+        const bool match_x = fabs(de_x - ad.de_x_mean) < n_si * ad.de_x_sigma;
+        const bool match_y = fabs(de_y - ad.de_y_mean) < n_si * ad.de_y_sigma;
+        if (match_x && match_y)
+          armEventData[arm].matched_track_idc[nsi].insert(idx_i);
+      }
+    }
+  }
+
+  // determine the number of reconstructed protons
+  for (unsigned int i = 0; i < hRecoProtonsMultiRP->size(); ++i) {
+    const auto &proton = (*hRecoProtonsMultiRP)[i];
+
+    CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->getRPId());
+    unsigned int arm = rpId.arm();
+
+    if (proton.validFit())
+      armEventData[arm].reco_proton_idc.insert(i);
+  }
+
+  // compare matched tracks with reco protons
+  if (verbosity_ > 1)
+    os << "* cmp matched tracks vs. reco protons" << std::endl;
+
+  for (auto &ap : armEventData) {
+    auto &ad = data_[ap.first];
+
+    if (verbosity_ > 1)
+      os << "    arm " << ap.first << std::endl;
+
+    for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
+      if (verbosity_ > 1)
+        os << "        nsi = " << nsi << std::endl;
+
+      for (const auto &tri : ap.second.matched_track_idc[nsi]) {
+        const auto &track = hTracks->at(tri);
+
+        bool some_proton_matching = false;
+
+        if (verbosity_ > 1)
+          os << "            tri = " << tri << std::endl;
+
+        for (const auto &pri : ap.second.reco_proton_idc) {
+          const auto &proton = (*hRecoProtonsMultiRP)[pri];
+
+          bool proton_matching = false;
+
+          for (const auto &pr_tr : proton.contributingLocalTracks()) {
+            CTPPSDetId rpId(pr_tr->getRPId());
+            unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
+
+            if (decRPId != ad.rpId_N)
+              continue;
+
+            const double x = pr_tr->getX();
+            const double y = pr_tr->getY();
+            const double th = 1E-3;  // 1 um
+
+            const bool match = (fabs(x - track.getX()) < th) && (fabs(y - track.getY()) < th);
+
+            if (verbosity_ > 1)
+              os << "                pri = " << pri << ": x_tr = " << track.getX() << ", x_pr = " << x
+                 << ", match = " << match << std::endl;
+
+            if (match) {
+              proton_matching = true;
+              break;
+            }
+          }
+
+          if (proton_matching) {
+            some_proton_matching = true;
+            break;
+          }
+        }
+
+        if (verbosity_ > 1)
+          os << "            --> some_proton_matching " << some_proton_matching << std::endl;
+
+        if (some_proton_matching)
+          ap.second.matched_track_with_prot_idc[nsi].insert(tri);
+        else
+          ap.second.matched_track_without_prot_idc[nsi].insert(tri);
+      }
+    }
+  }
+
+  // debug print
+  if (verbosity_ > 1) {
+    for (auto &ap : armEventData) {
+      auto &ad = data_[ap.first];
+
+      if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
+        continue;
+
+      os << "* results for arm " << ap.first << std::endl;
+
+      os << "    reco_proton_idc: ";
+      for (const auto &idx : ap.second.reco_proton_idc)
+        os << idx << ", ";
+      os << std::endl;
+
+      for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
+        os << "    n_si = " << ad.n_sigmas[nsi] << std::endl;
+
+        os << "        matched_track_idc: ";
+        for (const auto &idx : ap.second.matched_track_idc[nsi])
+          os << idx << ", ";
+        os << std::endl;
+
+        os << "        matched_track_with_prot_idc: ";
+        for (const auto &idx : ap.second.matched_track_with_prot_idc[nsi])
+          os << idx << ", ";
+        os << std::endl;
+
+        os << "        matched_track_without_prot_idc: ";
+        for (const auto &idx : ap.second.matched_track_without_prot_idc[nsi])
+          os << idx << ", ";
+        os << std::endl;
+      }
+    }
+  }
+
+  // update efficiency plots
+  for (auto &ap : armEventData) {
+    auto &ad = data_[ap.first];
+
+    // stop if sigmas not yet determined
+    if (ad.de_x_sigma <= 0. && ad.de_y_sigma <= 0.)
+      continue;
+
+    // loop over n_sigma choices
+    for (unsigned int nsi = 0; nsi < ad.n_sigmas.size(); ++nsi) {
+      const unsigned int n_exp_prot = ap.second.matched_track_idc[nsi].size();
+      const unsigned int n_rec_prot = ap.second.reco_proton_idc.size();
+
+      // stop if N(expected protons) out of range
+      if (n_exp_prot < 1 || n_exp_prot > ad.n_exp_prot_max)
+        continue;
+
+      // update method 1 plots
+      const double eff = double(n_rec_prot) / n_exp_prot;
+
+      for (unsigned int tri : ap.second.matched_track_idc[nsi]) {
+        const double x_N = hTracks->at(tri).getX();
+        const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
+
+        ad.effPlots[0][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
+        ad.effPlots[0][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
+
+        ad.effPlots[n_exp_prot][nsi].p_eff1_vs_x_N->Fill(x_N, eff);
+        ad.effPlots[n_exp_prot][nsi].p_eff1_vs_xi_N->Fill(xi_N, eff);
+      }
+
+      // update method 2 plots
+      for (const auto &tri : ap.second.matched_track_with_prot_idc[nsi]) {
+        const double x_N = hTracks->at(tri).getX();
+        const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
+
+        ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
+        ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
+
+        ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 1.);
+        ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 1.);
+      }
+
+      for (const auto &tri : ap.second.matched_track_without_prot_idc[nsi]) {
+        const double x_N = hTracks->at(tri).getX();
+        const double xi_N = ad.s_x_to_xi_N->Eval(x_N * 1E-1);  // conversion mm to cm
+
+        ad.effPlots[0][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
+        ad.effPlots[0][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
+
+        ad.effPlots[n_exp_prot][nsi].p_eff2_vs_x_N->Fill(x_N, 0.);
+        ad.effPlots[n_exp_prot][nsi].p_eff2_vs_xi_N->Fill(xi_N, 0.);
+      }
+    }
+  }
+
+  if (verbosity_)
+    edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorData") << os.str();
+}
+
+//----------------------------------------------------------------------------------------------------
+
+void CTPPSProtonReconstructionEfficiencyEstimatorData::endJob() {
+  auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
+
+  for (const auto &ait : data_) {
+    char buf[100];
+    sprintf(buf, "arm %u", ait.first);
+    TDirectory *d_arm = f_out->mkdir(buf);
+    gDirectory = d_arm;
+
+    ait.second.write();
+  }
+}
+
+//----------------------------------------------------------------------------------------------------
+
+DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorData);
diff --git a/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimator.cc b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc
similarity index 86%
rename from Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimator.cc
rename to Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc
index c63ff4b730dc8..8db9af320344e 100644
--- a/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimator.cc
+++ b/Validation/CTPPS/plugins/CTPPSProtonReconstructionEfficiencyEstimatorMC.cc
@@ -36,9 +36,9 @@
 
 //----------------------------------------------------------------------------------------------------
 
-class CTPPSProtonReconstructionEfficiencyEstimator : public edm::one::EDAnalyzer<> {
+class CTPPSProtonReconstructionEfficiencyEstimatorMC : public edm::one::EDAnalyzer<> {
 public:
-  explicit CTPPSProtonReconstructionEfficiencyEstimator(const edm::ParameterSet &);
+  explicit CTPPSProtonReconstructionEfficiencyEstimatorMC(const edm::ParameterSet &);
 
 private:
   void analyze(const edm::Event &, const edm::EventSetup &) override;
@@ -62,6 +62,8 @@ class CTPPSProtonReconstructionEfficiencyEstimator : public edm::one::EDAnalyzer
 
   std::string outputFile_;
 
+  unsigned int verbosity_;
+
   struct PlotGroup {
     std::unique_ptr<TProfile> p_eff_vs_xi;
 
@@ -90,7 +92,7 @@ using namespace HepMC;
 
 //----------------------------------------------------------------------------------------------------
 
-CTPPSProtonReconstructionEfficiencyEstimator::CTPPSProtonReconstructionEfficiencyEstimator(
+CTPPSProtonReconstructionEfficiencyEstimatorMC::CTPPSProtonReconstructionEfficiencyEstimatorMC(
     const edm::ParameterSet &iConfig)
     : tokenHepMCAfterSmearing_(
           consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMCAfterSmearing"))),
@@ -108,7 +110,10 @@ CTPPSProtonReconstructionEfficiencyEstimator::CTPPSProtonReconstructionEfficienc
       rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
       rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
 
-      outputFile_(iConfig.getParameter<string>("outputFile")) {
+      outputFile_(iConfig.getParameter<string>("outputFile")),
+
+      verbosity_(iConfig.getUntrackedParameter<unsigned int>("verbosity", 0))
+{
   rpDecId_near_[0] = rpId_45_N_;
   rpDecId_far_[0] = rpId_45_F_;
 
@@ -124,17 +129,8 @@ CTPPSProtonReconstructionEfficiencyEstimator::CTPPSProtonReconstructionEfficienc
 
 //----------------------------------------------------------------------------------------------------
 
-void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
-  bool verbosity = false;
-
-  const auto eid = iEvent.id().event();
-  //if (eid == 46 || eid == 741 || eid == 1649 || eid == 4223 || eid == 4279)
-  //  verbosity = true;
-
-  if (verbosity) {
-    printf("--------------------------------------------------\n");
-    printf("event %llu\n", eid);
-  }
+void CTPPSProtonReconstructionEfficiencyEstimatorMC::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
+  std::ostringstream os;
 
   // get conditions
   edm::ESHandle<LHCInfo> hLHCInfo;
@@ -220,8 +216,8 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv
   std::map<unsigned int, bool> isStripRPNear, isStripRPFar;
 
   for (auto &pp : particleInfo) {
-    if (verbosity)
-      printf("* barcode=%i, arm=%u, xi=%.3f\n", pp.first, pp.second.arm, pp.second.xi);
+    if (verbosity_)
+      os << "* barcode=" << pp.first << ", arm=" << pp.second.arm << ", xi=" << pp.second.xi << std::endl;
 
     for (const auto &rpp : pp.second.recHitsPerRP) {
       CTPPSDetId rpId(rpp.first);
@@ -247,18 +243,14 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv
           pp.second.inAcceptanceFar = true;
       }
 
-      if (verbosity)
-        printf("    RP %u: %u hits\n", rpDecId, rpp.second);
+      if (verbosity_)
+        os << "    RP " << rpDecId << ": " << rpp.second << " hits" << std::endl;
     }
 
     pp.second.inAcceptance = pp.second.inAcceptanceNear && pp.second.inAcceptanceFar;
 
-    if (verbosity) {
-      printf("    inAcceptance: near=%u, far=%u, global=%u\n",
-             pp.second.inAcceptanceNear,
-             pp.second.inAcceptanceFar,
-             pp.second.inAcceptance);
-    }
+    if (verbosity_)
+      os << "    inAcceptance: near=" << pp.second.inAcceptanceNear << ", far=" << pp.second.inAcceptanceFar << ", global=" << pp.second.inAcceptance << std::endl;
   }
 
   // count particles in acceptance
@@ -308,16 +300,10 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv
     const auto &npa = nParticlesInAcceptance[arm];
     const auto &nrt = nReconstructedTracks[arm];
 
-    if (verbosity) {
-      printf("* arm %u: nRecoProtons=%u (tracks near=%u, far=%u), nAcc=%u (near=%u, far=%u)\n",
-             arm,
-             nReconstructedProtons[arm],
-             nReconstructedTracks[arm].near,
-             nReconstructedTracks[arm].far,
-             npa.global,
-             npa.near,
-             npa.far);
-    }
+    if (verbosity_)
+      os << "* arm " << arm << ": nRecoProtons="
+        << nReconstructedProtons[arm] << " (tracks near=" << nReconstructedTracks[arm].near << ", far=" << nReconstructedTracks[arm].far
+        << "), nAcc=" << npa.global << " (near=" << npa.near << ", far=" << npa.far << ")" << std::endl;
 
     // skip event if no track in global acceptance
     if (npa.global < 1)
@@ -334,8 +320,8 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv
 
     const double eff = double(nReconstructedProtons[arm]) / npa.global;
 
-    if (verbosity)
-      printf("    eff=%.3f\n", eff);
+    if (verbosity_)
+      os << "    eff=" << eff << std::endl;
 
     for (auto &pp : particleInfo) {
       if (pp.second.arm != arm || !pp.second.inAcceptance)
@@ -344,11 +330,14 @@ void CTPPSProtonReconstructionEfficiencyEstimator::analyze(const edm::Event &iEv
       p.p_eff_vs_xi->Fill(pp.second.xi, eff);
     }
   }
+
+  if (verbosity_)
+    edm::LogInfo("CTPPSProtonReconstructionEfficiencyEstimatorMC") << os.str();
 }
 
 //----------------------------------------------------------------------------------------------------
 
-void CTPPSProtonReconstructionEfficiencyEstimator::endJob() {
+void CTPPSProtonReconstructionEfficiencyEstimatorMC::endJob() {
   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
 
   for (const auto &ait : plots_) {
@@ -368,4 +357,4 @@ void CTPPSProtonReconstructionEfficiencyEstimator::endJob() {
 
 //----------------------------------------------------------------------------------------------------
 
-DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimator);
+DEFINE_FWK_MODULE(CTPPSProtonReconstructionEfficiencyEstimatorMC);
diff --git a/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc b/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc
index 4e13a9eff4716..32ac8d9406da3 100644
--- a/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc
+++ b/Validation/CTPPS/plugins/CTPPSProtonReconstructionPlotter.cc
@@ -46,6 +46,18 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
   unsigned int rpId_45_N_, rpId_45_F_;
   unsigned int rpId_56_N_, rpId_56_F_;
 
+  struct AssociationCuts {
+    double ti_tr_min;
+    double ti_tr_max;
+
+    void load(const edm::ParameterSet &ps) {
+      ti_tr_min = ps.getParameter<double>("ti_tr_min");
+      ti_tr_max = ps.getParameter<double>("ti_tr_max");
+    }
+  };
+
+  std::map<unsigned int, AssociationCuts> association_cuts_;
+
   std::string outputFile_;
 
   signed int maxNonEmptyEvents_;
@@ -82,17 +94,19 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
     std::unique_ptr<TProfile> p_th_y_vs_xi;
 
     std::map<unsigned int, TH1D *> m_h_xi_nTracks;
+    std::unique_ptr<TH1D> h_xi_n1f1;
 
     SingleRPPlots()
         : h_multiplicity(new TH1D("", ";reconstructed protons", 11, -0.5, 10.5)),
           h_xi(new TH1D("", ";#xi", 100, 0., 0.3)),
           h2_th_y_vs_xi(new TH2D("", ";#xi;#theta_{y}   (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)),
-          p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y}   (rad)", 100, 0., 0.3)) {
+          p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y}   (rad)", 100, 0., 0.3)),
+          h_xi_n1f1(new TH1D("", ";#xi", 100, 0., 0.3)) {
       for (unsigned int n = 2; n <= 10; ++n)
         m_h_xi_nTracks[n] = new TH1D(*h_xi);
     }
 
-    void fill(const reco::ForwardProton &p, unsigned int nTracks) {
+    void fill(const reco::ForwardProton &p, unsigned int nTracks, bool n1f1) {
       if (p.validFit()) {
         h_xi->Fill(p.xi());
 
@@ -103,6 +117,9 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
         auto it = m_h_xi_nTracks.find(nTracks);
         if (it != m_h_xi_nTracks.end())
           it->second->Fill(p.xi());
+
+        if (n1f1)
+          h_xi_n1f1->Fill(p.xi());
       }
     }
 
@@ -123,6 +140,8 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
       }
 
       gDirectory = d_top;
+
+      h_xi_n1f1->Write("h_xi_n1f1");
     }
   };
 
@@ -140,9 +159,11 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
     std::unique_ptr<TH2D> h2_timing_tracks_vs_prot_mult;
 
     std::map<unsigned int, TH1D *> m_h_xi_nTracks;
+    std::unique_ptr<TH1D> h_xi_n1f1;
 
     std::unique_ptr<TH1D> h_de_x_timing_vs_tracking, h_de_x_rel_timing_vs_tracking, h_de_x_match_timing_vs_tracking;
-    std::unique_ptr<TH1D> h_de_x_rel_timing_vs_tracking_ClCo;
+    std::unique_ptr<TH1D> h_de_x_timing_vs_tracking_ClCo, h_de_x_rel_timing_vs_tracking_ClCo,
+        h_de_x_match_timing_vs_tracking_ClCo;
 
     std::unique_ptr<TH2D> h2_y_vs_x_tt0_ClCo, h2_y_vs_x_tt1_ClCo, h2_y_vs_x_ttm_ClCo;
 
@@ -150,7 +171,7 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
         : h_multiplicity(new TH1D("", ";reconstructed protons per event", 11, -0.5, 10.5)),
           h_xi(new TH1D("", ";#xi", 100, 0., 0.3)),
           h_th_x(new TH1D("", ";#theta_{x}   (rad)", 250, -500E-6, +500E-6)),
-          h_th_y(new TH1D("", ";#theta_{y}   (rad)", 250, -500E-6, +500E-6)),
+          h_th_y(new TH1D("", ";#theta_{y}   (rad)", 500, -1000E-6, +1000E-6)),
           h_vtx_y(new TH1D("", ";vtx_{y}   (cm)", 100, -100E-3, +100E-3)),
           h_chi_sq(new TH1D("", ";#chi^{2}", 100, 0., 10.)),
           h_log_chi_sq(new TH1D("", ";log_{10} #chi^{2}", 100, -20., 5.)),
@@ -159,17 +180,24 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
           h_n_contrib_tracking_tracks(new TH1D("", ";n of contrib. tracking tracks per reco proton", 4, -0.5, +3.5)),
           h_n_contrib_timing_tracks(new TH1D("", ";n of contrib. timing tracks per reco proton", 4, -0.5, +3.5)),
           h2_th_x_vs_xi(new TH2D("", ";#xi;#theta_{x}   (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)),
-          h2_th_y_vs_xi(new TH2D("", ";#xi;#theta_{y}   (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)),
+          h2_th_y_vs_xi(new TH2D("", ";#xi;#theta_{y}   (rad)", 100, 0., 0.3, 200, -1000E-6, +1000E-6)),
           h2_vtx_y_vs_xi(new TH2D("", ";#xi;vtx_{y}   (cm)", 100, 0., 0.3, 100, -100E-3, +100E-3)),
           p_th_x_vs_xi(new TProfile("", ";#xi;#theta_{x}   (rad)", 100, 0., 0.3)),
           p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y}   (rad)", 100, 0., 0.3)),
           p_vtx_y_vs_xi(new TProfile("", ";#xi;vtx_{y}   (cm)", 100, 0., 0.3)),
           h2_timing_tracks_vs_prot_mult(
               new TH2D("", ";reco protons per event;timing tracks per event", 11, -0.5, 10.5, 11, -0.5, 10.5)),
+          h_xi_n1f1(new TH1D("", ";#xi", 100, 0., 0.3)),
+
           h_de_x_timing_vs_tracking(new TH1D("", ";#Delta x   (mm)", 200, -1., +1.)),
           h_de_x_rel_timing_vs_tracking(new TH1D("", ";#Delta x / #sigma(x)", 200, -20., +20.)),
           h_de_x_match_timing_vs_tracking(new TH1D("", ";match between tracking and timing tracks", 2, -0.5, +1.5)),
+
+          h_de_x_timing_vs_tracking_ClCo(new TH1D("", ";#Delta x   (mm)", 200, -1., +1.)),
           h_de_x_rel_timing_vs_tracking_ClCo(new TH1D("", ";#Delta x / #sigma(x)", 200, -20., +20.)),
+          h_de_x_match_timing_vs_tracking_ClCo(
+              new TH1D("", ";match between tracking and timing tracks", 2, -0.5, +1.5)),
+
           h2_y_vs_x_tt0_ClCo(new TH2D("", ";x   (mm);y   (mm)", 100, -5., 25., 100, -15., +15.)),
           h2_y_vs_x_tt1_ClCo(new TH2D("", ";x   (mm);y   (mm)", 100, -5., 25., 100, -15., +15.)),
           h2_y_vs_x_ttm_ClCo(new TH2D("", ";x   (mm);y   (mm)", 100, -5., 25., 100, -15., +15.)) {
@@ -191,7 +219,7 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
         m_h_xi_nTracks[n] = new TH1D(*h_xi);
     }
 
-    void fill(const reco::ForwardProton &p, unsigned int nTracks) {
+    void fill(const reco::ForwardProton &p, unsigned int nTracks, bool n1f1) {
       if (!p.validFit())
         return;
 
@@ -246,6 +274,9 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
       auto it = m_h_xi_nTracks.find(nTracks);
       if (it != m_h_xi_nTracks.end())
         it->second->Fill(p.xi());
+
+      if (n1f1)
+        h_xi_n1f1->Fill(p.xi());
     }
 
     void write() const {
@@ -306,10 +337,15 @@ class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
 
       gDirectory = d_top;
 
+      h_xi_n1f1->Write("h_xi_n1f1");
+
       h_de_x_timing_vs_tracking->Write("h_de_x_timing_vs_tracking");
       h_de_x_rel_timing_vs_tracking->Write("h_de_x_rel_timing_vs_tracking");
       h_de_x_match_timing_vs_tracking->Write("h_de_x_match_timing_vs_tracking");
+
+      h_de_x_timing_vs_tracking_ClCo->Write("h_de_x_timing_vs_tracking_ClCo");
       h_de_x_rel_timing_vs_tracking_ClCo->Write("h_de_x_rel_timing_vs_tracking_ClCo");
+      h_de_x_match_timing_vs_tracking_ClCo->Write("h_de_x_match_timing_vs_tracking_ClCo");
 
       h2_y_vs_x_tt0_ClCo->Write("h2_y_vs_x_tt0_ClCo");
       h2_y_vs_x_tt1_ClCo->Write("h2_y_vs_x_tt1_ClCo");
@@ -444,7 +480,12 @@ CTPPSProtonReconstructionPlotter::CTPPSProtonReconstructionPlotter(const edm::Pa
       p_y_L_diffNF_vs_y_L_N_(new TProfile("p_y_L_diffNF_vs_y_L_N", ";y_{LN};y_{LF} - y_{LN}", 100, -20., +20.)),
       p_y_R_diffNF_vs_y_R_N_(new TProfile("p_y_R_diffNF_vs_y_R_N", ";y_{RN};y_{RF} - y_{RN}", 100, -20., +20.)),
 
-      n_non_empty_events_(0) {}
+      n_non_empty_events_(0) {
+  for (const std::string &sector : {"45", "56"}) {
+    const unsigned int arm = (sector == "45") ? 0 : 1;
+    association_cuts_[arm].load(ps.getParameterSet("association_cuts_" + sector));
+  }
+}
 
 //----------------------------------------------------------------------------------------------------
 
@@ -461,10 +502,12 @@ void CTPPSProtonReconstructionPlotter::CalculateTimingTrackingDistance(const rec
   const double z_j = geometry.getRPTranslation(tr_j.getRPId()).z();
 
   // interpolation from tracking RPs
-  const double z_ti = -geometry.getRPTranslation(tr_ti.getRPId()).z();  // the minus sign fixes a bug in the diamond geometry
+  const double z_ti =
+      -geometry.getRPTranslation(tr_ti.getRPId()).z();  // the minus sign fixes a bug in the diamond geometry
   const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
   const double x_inter = f_i * tr_i.getX() + f_j * tr_j.getX();
-  const double x_inter_unc_sq = f_i * f_i * tr_i.getXUnc() * tr_i.getXUnc() + f_j * f_j * tr_j.getXUnc() * tr_j.getXUnc();
+  const double x_inter_unc_sq =
+      f_i * f_i * tr_i.getXUnc() * tr_i.getXUnc() + f_j * f_j * tr_j.getXUnc() * tr_j.getXUnc();
 
   // save distance
   de_x = tr_ti.getX() - x_inter;
@@ -500,19 +543,28 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
   const CTPPSLocalTrackLite *tr_R_N = nullptr;
   const CTPPSLocalTrackLite *tr_R_F = nullptr;
   std::map<unsigned int, unsigned int> armTrackCounter, armTimingTrackCounter;
+  std::map<unsigned int, unsigned int> armTrackCounter_N, armTrackCounter_F;
 
   for (const auto &tr : *hTracks) {
     CTPPSDetId rpId(tr.getRPId());
     unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
 
-    if (decRPId == rpId_45_N_)
+    if (decRPId == rpId_45_N_) {
       tr_L_N = &tr;
-    if (decRPId == rpId_45_F_)
+      armTrackCounter_N[0]++;
+    }
+    if (decRPId == rpId_45_F_) {
       tr_L_F = &tr;
-    if (decRPId == rpId_56_N_)
+      armTrackCounter_F[0]++;
+    }
+    if (decRPId == rpId_56_N_) {
       tr_R_N = &tr;
-    if (decRPId == rpId_56_F_)
+      armTrackCounter_N[1]++;
+    }
+    if (decRPId == rpId_56_F_) {
       tr_R_F = &tr;
+      armTrackCounter_F[1]++;
+    }
 
     armTrackCounter[rpId.arm()]++;
 
@@ -542,7 +594,10 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
   for (const auto &proton : *hRecoProtonsSingleRP) {
     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->getRPId());
     unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
-    singleRPPlots_[decRPId].fill(proton, armTrackCounter[rpId.arm()]);
+
+    const bool n1f1 = (armTrackCounter_N[rpId.arm()] == 1 && armTrackCounter_F[rpId.arm()] == 1);
+
+    singleRPPlots_[decRPId].fill(proton, armTrackCounter[rpId.arm()], n1f1);
 
     if (proton.validFit())
       singleRPMultiplicity[decRPId]++;
@@ -555,7 +610,10 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
   for (const auto &proton : *hRecoProtonsMultiRP) {
     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->getRPId());
     unsigned int armId = rpId.arm();
-    multiRPPlots_[armId].fill(proton, armTrackCounter[armId]);
+
+    const bool n1f1 = (armTrackCounter_N[armId] == 1 && armTrackCounter_F[armId] == 1);
+
+    multiRPPlots_[armId].fill(proton, armTrackCounter[armId], n1f1);
 
     if (proton.validFit())
       multiRPMultiplicity[armId]++;
@@ -594,15 +652,18 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
       double de_x = 0., de_x_unc = 0.;
       CalculateTimingTrackingDistance(proton, tr, *hGeometry, de_x, de_x_unc);
 
-      double rd = (de_x_unc > 0.) ? de_x / de_x_unc : -1E10;
+      const double rd = (de_x_unc > 0.) ? de_x / de_x_unc : -1E10;
+      const auto &ac = association_cuts_[armId];
+      const bool match = (ac.ti_tr_min <= fabs(rd) && fabs(rd) <= ac.ti_tr_max);
 
       pl.h_de_x_timing_vs_tracking->Fill(de_x);
       pl.h_de_x_rel_timing_vs_tracking->Fill(rd);
-      pl.h_de_x_match_timing_vs_tracking->Fill(fabs(de_x / de_x_unc) <= 1. ? 1. : 0.);
+      pl.h_de_x_match_timing_vs_tracking->Fill(match ? 1. : 0.);
 
-      if (clCo[armId]) {
-        if (armTimingTrackCounter[armId] == 1)
-          pl.h_de_x_rel_timing_vs_tracking_ClCo->Fill(rd);
+      if (clCo[armId] && armTimingTrackCounter[armId] == 1) {
+        pl.h_de_x_timing_vs_tracking_ClCo->Fill(de_x);
+        pl.h_de_x_rel_timing_vs_tracking_ClCo->Fill(rd);
+        pl.h_de_x_match_timing_vs_tracking_ClCo->Fill(match ? 1. : 0.);
       }
     }
   }
@@ -645,13 +706,6 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
 
     const reco::ForwardProton *p_s_N = nullptr, *p_s_F = nullptr;
 
-    /*
-    printf("multi-RP candidate: ");
-    for (const auto &tr_m : proton_m.contributingLocalTracks())
-      printf("%u, ", tr_m.key());
-    printf("\n");
-    */
-
     for (const auto &proton_s : *hRecoProtonsSingleRP) {
       // skip if source of single-RP reco not included in multi-RP reco
       const auto key_s = proton_s.contributingLocalTracks()[0].key();
@@ -663,8 +717,6 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
         }
       }
 
-      //printf("    key_s = %u --> compatible = %u\n", key_s, compatible);
-
       if (!compatible)
         continue;
 
@@ -673,8 +725,6 @@ void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const ed
       const unsigned int idx = rpId_s.arm() * 1000 + rpId_s.station() * 100 + rpId_s.rp() * 10 + rpId_s.arm();
       singleMultiCorrelationPlots_[idx].fill(proton_s, proton_m);
 
-      //printf("    arm_s = %u, arm_m = %u\n", rpId_s.arm(), arm);
-
       // select protons for arm-correlation plots
       const unsigned int rpDecId_s = rpId_s.arm() * 100 + rpId_s.station() * 10 + rpId_s.rp();
       if (rpDecId_s == rpId_45_N_ || rpDecId_s == rpId_56_N_)
diff --git a/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc b/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc
index 2b47adde51d29..0d65bcaf481d7 100644
--- a/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc
+++ b/Validation/CTPPS/plugins/CTPPSTrackDistributionPlotter.cc
@@ -23,6 +23,7 @@
 #include "TH2D.h"
 #include "TProfile.h"
 #include "TProfile2D.h"
+#include "TGraph.h"
 
 #include <map>
 
@@ -44,6 +45,9 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
 
   std::string outputFile_;
 
+  unsigned int events_total_;
+  std::map<unsigned int, unsigned int> events_per_arm_;
+
   struct RPPlots {
     bool initialized;
 
@@ -51,13 +55,14 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
     std::unique_ptr<TProfile> p_y_vs_x;
     std::unique_ptr<TH1D> h_x;
     std::unique_ptr<TH1D> h_y;
+    std::unique_ptr<TH1D> h_time;
 
     RPPlots() : initialized(false) {}
 
     void init(bool pixel, double pitch) {
       const double bin_size_x = (pixel) ? pitch * cos(18.4 / 180. * M_PI) : 100E-3;
 
-      h2_y_vs_x.reset(new TH2D("", "", 300, -10., +70., 300, -30., +30.));
+      h2_y_vs_x.reset(new TH2D("", "", 300, -10., +70., 600, -30., +30.));
       p_y_vs_x.reset(new TProfile("", "", 300, -10., +70.));
 
       int n_mi = ceil(10. / bin_size_x);
@@ -67,14 +72,17 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
 
       h_y.reset(new TH1D("", "", 300, -15., +15.));
 
+      h_time.reset(new TH1D("", ";time", 500, -50., +50.));
+
       initialized = true;
     }
 
-    void fill(double x, double y) {
+    void fill(double x, double y, double time) {
       h2_y_vs_x->Fill(x, y);
       p_y_vs_x->Fill(x, y);
       h_x->Fill(x);
       h_y->Fill(y);
+      h_time->Fill(time);
     }
 
     void write() const {
@@ -82,15 +90,26 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
       p_y_vs_x->Write("p_y_vs_x");
       h_x->Write("h_x");
       h_y->Write("h_y");
+      h_time->Write("h_time");
     }
   };
 
   std::map<unsigned int, RPPlots> rpPlots;
 
   struct ArmPlots {
+    unsigned int rpId_N, rpId_F;
+
     std::unique_ptr<TH1D> h_de_x, h_de_y;
     std::unique_ptr<TProfile> p_de_x_vs_x, p_de_y_vs_x;
     std::unique_ptr<TProfile2D> p2_de_x_vs_x_y, p2_de_y_vs_x_y;
+    std::unique_ptr<TH2D> h2_de_x_vs_x, h2_de_y_vs_x;
+    std::unique_ptr<TH2D> h2_de_y_vs_de_x;
+
+    struct Stat {
+      unsigned int sN = 0, s1 = 0;
+    };
+
+    std::map<unsigned int, std::map<unsigned int, Stat>> m_stat;
 
     ArmPlots()
         : h_de_x(new TH1D("", ";x^{F} - x^{N}", 100, -1., +1.)),
@@ -98,7 +117,10 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
           p_de_x_vs_x(new TProfile("", ";x^{N};x^{F} - x^{N}", 40, 0., 40.)),
           p_de_y_vs_x(new TProfile("", ";x^{N};y^{F} - y^{N}", 40, 0., 40.)),
           p2_de_x_vs_x_y(new TProfile2D("", ";x;y", 40, 0., 40., 40, -20., +20.)),
-          p2_de_y_vs_x_y(new TProfile2D("", ";x;y", 40, 0., 40., 40, -20., +20.)) {}
+          p2_de_y_vs_x_y(new TProfile2D("", ";x;y", 40, 0., 40., 40, -20., +20.)),
+          h2_de_x_vs_x(new TH2D("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., 100, -1., +1.)),
+          h2_de_y_vs_x(new TH2D("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., 100, -1., +1.)),
+          h2_de_y_vs_de_x(new TH2D("", ";x^{F} - x^{N};y^{F} - y^{N}", 100, -1., +1., 100, -1., +1.)) {}
 
     void fill(double x_N, double y_N, double x_F, double y_F) {
       h_de_x->Fill(x_F - x_N);
@@ -109,6 +131,11 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
 
       p2_de_x_vs_x_y->Fill(x_N, y_N, x_F - x_N);
       p2_de_y_vs_x_y->Fill(x_N, y_N, y_F - y_N);
+
+      h2_de_x_vs_x->Fill(x_N, x_F - x_N);
+      h2_de_y_vs_x->Fill(x_N, y_F - y_N);
+
+      h2_de_y_vs_de_x->Fill(x_F - x_N, y_F - y_N);
     }
 
     void write() const {
@@ -120,6 +147,27 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
 
       p2_de_x_vs_x_y->Write("p2_de_x_vs_x_y");
       p2_de_y_vs_x_y->Write("p2_de_y_vs_x_y");
+
+      h2_de_x_vs_x->Write("h2_de_x_vs_x");
+      h2_de_y_vs_x->Write("h2_de_y_vs_x");
+
+      h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x");
+
+      for (const auto& rp : m_stat) {
+        TGraph* g = new TGraph();
+
+        char buf[100];
+        sprintf(buf, "g_mean_track_mult_run_%u", rp.first);
+        g->SetName(buf);
+
+        for (const auto& lsp : rp.second) {
+          const int idx = g->GetN();
+          const double m = (lsp.second.s1 > 0) ? double(lsp.second.sN) / lsp.second.s1 : 0.;
+          g->SetPoint(idx, lsp.first, m);
+        }
+
+        g->Write();
+      }
     }
   };
 
@@ -131,7 +179,14 @@ class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
 CTPPSTrackDistributionPlotter::CTPPSTrackDistributionPlotter(const edm::ParameterSet& iConfig)
     : tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
       x_pitch_pixels_(iConfig.getUntrackedParameter<double>("x_pitch_pixels", 150E-3)),
-      outputFile_(iConfig.getParameter<std::string>("outputFile")) {}
+      outputFile_(iConfig.getParameter<std::string>("outputFile")),
+      events_total_(0) {
+  armPlots[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
+  armPlots[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
+
+  armPlots[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
+  armPlots[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
+}
 
 //----------------------------------------------------------------------------------------------------
 
@@ -141,6 +196,8 @@ void CTPPSTrackDistributionPlotter::analyze(const edm::Event& iEvent, const edm:
   iEvent.getByToken(tracksToken_, tracks);
 
   // process tracks
+  std::map<unsigned int, unsigned int> m_mult;
+
   for (const auto& trk : *tracks) {
     CTPPSDetId rpId(trk.getRPId());
     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
@@ -150,27 +207,53 @@ void CTPPSTrackDistributionPlotter::analyze(const edm::Event& iEvent, const edm:
     if (!pl.initialized)
       pl.init(rpPixel, x_pitch_pixels_);
 
-    pl.fill(trk.getX(), trk.getY());
+    pl.fill(trk.getX(), trk.getY(), trk.getTime());
+
+    m_mult[rpId.arm()]++;
+  }
+
+  for (unsigned int arm = 0; arm < 2; ++arm) {
+    auto& st = armPlots[arm].m_stat[iEvent.id().run()][iEvent.id().luminosityBlock()];
+    st.s1++;
+    st.sN += m_mult[arm];
   }
 
   for (const auto& t1 : *tracks) {
-    CTPPSDetId rpId1(t1.getRPId());
+    const CTPPSDetId rpId1(t1.getRPId());
 
     for (const auto& t2 : *tracks) {
-      CTPPSDetId rpId2(t2.getRPId());
+      const CTPPSDetId rpId2(t2.getRPId());
 
       if (rpId1.arm() != rpId2.arm())
         continue;
 
-      if (rpId1.station() == 0 && rpId2.station() == 2)
-        armPlots[rpId1.arm()].fill(t1.getX(), t1.getY(), t2.getX(), t2.getY());
+      auto& ap = armPlots[rpId1.arm()];
+
+      const unsigned int rpDecId1 = rpId1.arm() * 100 + rpId1.station() * 10 + rpId1.rp();
+      const unsigned int rpDecId2 = rpId2.arm() * 100 + rpId2.station() * 10 + rpId2.rp();
+
+      if (rpDecId1 == ap.rpId_N && rpDecId2 == ap.rpId_F)
+        ap.fill(t1.getX(), t1.getY(), t2.getX(), t2.getY());
     }
   }
+
+  // update counters
+  events_total_++;
+
+  if (m_mult[0] > 0)
+    events_per_arm_[0]++;
+  if (m_mult[1] > 0)
+    events_per_arm_[1]++;
 }
 
 //----------------------------------------------------------------------------------------------------
 
 void CTPPSTrackDistributionPlotter::endJob() {
+  edm::LogInfo("CTPPSTrackDistributionPlotter")
+      << "    events processed: " << events_total_ << " (" << std::scientific << double(events_total_) << ")\n"
+      << "    events with tracks in sector 45: " << events_per_arm_[0] << " (" << double(events_per_arm_[0]) << ")\n"
+      << "    events with tracks in sector 56: " << events_per_arm_[1] << " (" << double(events_per_arm_[1]) << ")";
+
   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
 
   for (const auto& it : rpPlots) {