Merge pull request #29770 from CTPPS/pps_re_miniaod_11_1

PPS: updates for re-miniAOD
cms-sw · May 20, 2020 · ae84fe5 · ae84fe5
2 parents 94a56f1 + 00d5c66
commit ae84fe5
Show file tree

Hide file tree

Showing 10 changed files with 1,077 additions and 102 deletions.
diff --git a/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc b/RecoPPS/ProtonReconstruction/plugins/CTPPSProtonProducer.cc
@@ -44,6 +44,8 @@ class CTPPSProtonProducer : public edm::stream::EDProducer<> {
 
   edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
 
+  bool pixelDiscardBXShiftedTracks_;
+
   std::string lhcInfoLabel_;
   std::string opticsLabel_;
 
@@ -59,25 +61,32 @@ 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");
@@ -88,12 +97,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");
@@ -106,6 +122,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_;
 
@@ -117,6 +134,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)),
@@ -131,6 +151,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"),
@@ -157,6 +178,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");
 
@@ -187,6 +211,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");
 
@@ -253,6 +279,12 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe
             tr.ty() > localAngleYMax_)
           continue;
 
+        if (pixelDiscardBXShiftedTracks_) {
+          if (tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::allShiftedPlanes ||
+              tr.pixelTrackRecoInfo() == CTPPSpixelLocalTrackReconstructionInfo::mixedPlanes)
+            continue;
+        }
+
         const CTPPSDetId rpId(tr.rpId());
 
         if (verbosity_)
@@ -298,32 +330,36 @@ 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->rpTranslation(tr_i.rpId()).z();
+              const double z_j = hGeometry->rpTranslation(tr_j.rpId()).z();
+
               const auto &pr_i = singleRPResultsIndexed[i];
               const auto &pr_j = singleRPResultsIndexed[j];
 
               if (tr_i.rpId() == tr_j.rpId())
                 continue;
 
+              if (std::abs(z_i) >= std::abs(z_j))
+                continue;
+
               bool matching = true;
 
-              if (ac.x_cut_apply && std::abs(tr_i.x() - tr_j.x()) > ac.x_cut_value)
+              if (ac.x_cut_apply && std::abs(tr_i.x() - tr_j.x() - ac.x_cut_mean) > ac.x_cut_value)
                 matching = false;
-              if (ac.y_cut_apply && std::abs(tr_i.y() - tr_j.y()) > ac.y_cut_value)
+              else if (ac.y_cut_apply && std::abs(tr_i.y() - tr_j.y() - 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)
@@ -432,7 +468,7 @@ void CTPPSProtonProducer::produce(edm::Event &iEvent, const edm::EventSetup &iSe
               swt += w * tr.time();
             }
 
-            float time = 0., time_unc = 0.;
+            float time = default_time_, time_unc = 0.;
             if (sw > 0.) {
               time = swt / sw;
               time_unc = 1. / sqrt(sw);

diff --git a/RecoPPS/ProtonReconstruction/python/ctppsProtons_cff.py b/RecoPPS/ProtonReconstruction/python/ctppsProtons_cff.py
@@ -20,28 +20,73 @@ def applyDefaultSettings(ctppsProtons):
   ctppsProtons.association_cuts_45.xi_cut_apply = True
   ctppsProtons.association_cuts_45.xi_cut_value = 0.010
   ctppsProtons.association_cuts_45.th_y_cut_apply = False
+  ctppsProtons.association_cuts_45.ti_tr_min = -1.5
+  ctppsProtons.association_cuts_45.ti_tr_max = 2.0
 
   ctppsProtons.association_cuts_56.x_cut_apply = False
   ctppsProtons.association_cuts_56.y_cut_apply = False
   ctppsProtons.association_cuts_56.xi_cut_apply = True
   ctppsProtons.association_cuts_56.xi_cut_value = 0.015
   ctppsProtons.association_cuts_56.th_y_cut_apply = False
+  ctppsProtons.association_cuts_56.ti_tr_min = -1.5
+  ctppsProtons.association_cuts_56.ti_tr_max = 2.0
+
+  ctppsProtons.pixelDiscardBXShiftedTracks = True
+  ctppsProtons.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
+  ctppsProtons.association_cuts_45.x_cut_apply = False
+  ctppsProtons.association_cuts_45.y_cut_apply = False
+
+  ctppsProtons.association_cuts_45.xi_cut_apply = True
+  ctppsProtons.association_cuts_45.xi_cut_value = 5. * 0.00121
+  ctppsProtons.association_cuts_45.xi_cut_mean = +6.0695e-5
+
+  ctppsProtons.association_cuts_45.th_y_cut_apply = False
+
+  ctppsProtons.association_cuts_56.x_cut_apply = False
+
+  ctppsProtons.association_cuts_56.y_cut_apply = True
+  ctppsProtons.association_cuts_56.y_cut_value = 5. * 0.14777
+  ctppsProtons.association_cuts_56.y_cut_mean = -0.022612
+
+  ctppsProtons.association_cuts_56.xi_cut_apply = True
+  ctppsProtons.association_cuts_56.xi_cut_value = 5. * 0.0020627
+  ctppsProtons.association_cuts_56.xi_cut_mean = +8.012857e-5
+
+  ctppsProtons.association_cuts_56.th_y_cut_apply = False
 
 ctpps_2017.toModify(ctppsProtons, apply2017Settings)
 
 def apply2018Settings(ctppsProtons):
-  ctppsProtons.association_cuts_45.xi_cut_value = 0.013
-  ctppsProtons.association_cuts_45.th_y_cut_apply = True
-  ctppsProtons.association_cuts_45.th_y_cut_value = 30E-6
+  ctppsProtons.association_cuts_45.x_cut_apply = True
+  ctppsProtons.association_cuts_45.x_cut_value = 4. * 0.16008188
+  ctppsProtons.association_cuts_45.x_cut_mean = -0.065194856
+
+  ctppsProtons.association_cuts_45.y_cut_apply = True
+  ctppsProtons.association_cuts_45.y_cut_value = 4. * 0.1407986
+  ctppsProtons.association_cuts_45.y_cut_mean = +0.10973631
 
-  ctppsProtons.association_cuts_56.xi_cut_value = 0.013
-  ctppsProtons.association_cuts_56.th_y_cut_apply = True
-  ctppsProtons.association_cuts_56.th_y_cut_value = 20E-6
+  ctppsProtons.association_cuts_45.xi_cut_apply = True
+  ctppsProtons.association_cuts_45.xi_cut_value = 4. * 0.0012403586
+  ctppsProtons.association_cuts_45.xi_cut_mean = +3.113062e-5
+
+  ctppsProtons.association_cuts_45.th_y_cut_apply = False
+
+  ctppsProtons.association_cuts_56.x_cut_apply = True
+  ctppsProtons.association_cuts_56.x_cut_value = 5. * 0.18126434
+  ctppsProtons.association_cuts_56.x_cut_mean = +0.073016431
+
+  ctppsProtons.association_cuts_56.y_cut_apply = True
+  ctppsProtons.association_cuts_56.y_cut_value = 5. * 0.14990802
+  ctppsProtons.association_cuts_56.y_cut_mean = +0.064261029
+
+  ctppsProtons.association_cuts_56.xi_cut_apply = True
+  ctppsProtons.association_cuts_56.xi_cut_value = 5. * 0.002046409
+  ctppsProtons.association_cuts_56.xi_cut_mean = -1.1852528e-5
+
+  ctppsProtons.association_cuts_56.th_y_cut_apply = False
 
 ctpps_2018.toModify(ctppsProtons, apply2018Settings)
diff --git a/RecoPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc b/RecoPPS/ProtonReconstruction/src/ProtonReconstructionAlgorithm.cc
@@ -382,7 +382,8 @@ reco::ForwardProton ProtonReconstructionAlgorithm::reconstructFromMultiRP(const
 
   // print results
   if (verbosity_)
-    os << "    xi=" << xi << ", th_x=" << th_x << ", th_y=" << th_y << ", vtx_y=" << vtx_y << ", chiSq = " << chi2
+    os << "    fit valid=" << valid << std::endl
+       << "    xi=" << xi << ", th_x=" << th_x << ", th_y=" << th_y << ", vtx_y=" << vtx_y << ", chiSq = " << chi2
        << std::endl;
 
   // save reco candidate

diff --git 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
@@ -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();
 }
 
 //----------------------------------------------------------------------------------------------------