[PPS] Prompt calibration loop component for timing detector

CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h

@@ -0,0 +1,17 @@
+#ifndef CalibPPS_TimingCalibration_TimingCalibrationStruct_h
+#define CalibPPS_TimingCalibration_TimingCalibrationStruct_h
+
+#include "DQMServices/Core/interface/DQMStore.h"
+#include <unordered_map>
+
+struct TimingCalibrationHistograms {
+public:
+  TimingCalibrationHistograms() = default;
+
+  using MonitorMap = std::unordered_map<uint32_t, dqm::reco::MonitorElement*>;
+
+  MonitorMap leadingTime, toT;
+  MonitorMap leadingTimeVsToT;
+};
+
+#endif

CalibPPS/TimingCalibration/plugins/BuildFile.xml

@@ -0,0 +1,9 @@
+<library name="CalibPPSTimingCalibrationPlugins" file="*.cc">
+  <flags EDM_PLUGIN="1"/>
+  <use name="DQMServices/Core"/>
+  <use name="FWCore/ServiceRegistry"/>
+  <use name="CondCore/DBOutputService"/>
+  <use name="CondFormats/PPSObjects"/>
+  <use name="Geometry/Records"/>
+  <use name="Geometry/VeryForwardGeometryBuilder"/>
+</library>

CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc

@@ -0,0 +1,160 @@
+/****************************************************************************
+ *
+ * This is a part of PPS offline software.
+ * Authors:
+ *   Edoardo Bossini
+ *   Piotr Maciej Cwiklicki
+ *   Laurent Forthomme
+ *
+ ****************************************************************************/
+
+#include "DQMServices/Core/interface/DQMEDHarvester.h"
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+
+#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
+#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
+
+#include "CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h"
+#include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
+
+#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h"
+#include "CondFormats/PPSObjects/interface/PPSTimingCalibration.h"
+
+//------------------------------------------------------------------------------
+
+class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester {
+public:
+  PPSTimingCalibrationPCLHarvester(const edm::ParameterSet&);
+  void beginRun(const edm::Run&, const edm::EventSetup&) override;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+private:
+  void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
+  edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> geomEsToken_;
+  std::vector<CTPPSDiamondDetId> detids_;
+  const std::string dqmDir_;
+  const std::string formula_;
+  const unsigned int min_entries_;
+  TF1 interp_;
+};
+
+//------------------------------------------------------------------------------
+
+PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::ParameterSet& iConfig)
+    : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
+      dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
+      formula_(iConfig.getParameter<std::string>("formula")),
+      min_entries_(iConfig.getParameter<unsigned int>("minEntries")),
+      interp_("interp", formula_.c_str(), 10.5, 25.) {
+  // first ensure DB output service is available
+  edm::Service<cond::service::PoolDBOutputService> poolDbService;
+  if (!poolDbService.isAvailable())
+    throw cms::Exception("PPSTimingCalibrationPCLHarvester") << "PoolDBService required";
+
+  // constrain the min/max fit values
+  interp_.SetParLimits(1, 9., 15.);
+  interp_.SetParLimits(2, 0.2, 2.5);
+}
+
+//------------------------------------------------------------------------------
+
+void PPSTimingCalibrationPCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
+  const auto& geom = iSetup.getData(geomEsToken_);
+  for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) {
+    if (!CTPPSDiamondDetId::check(it->first))
+      continue;
+    const CTPPSDiamondDetId detid(it->first);
+    if (detid.station() == 1)  // for the time being, only compute for this station (run 2 diamond)
+      detids_.emplace_back(detid);
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
+  // book the parameters containers
+  PPSTimingCalibration::ParametersMap calib_params;
+  PPSTimingCalibration::TimingMap calib_time;
+
+  iGetter.cd();
+  iGetter.setCurrentFolder(dqmDir_);
+
+  // compute the fit parameters for all monitored channels
+  TimingCalibrationHistograms hists;
+  std::string ch_name;
+  for (const auto& detid : detids_) {
+    detid.channelName(ch_name);
+    const auto chid = detid.rawId();
+    const PPSTimingCalibration::Key key{
+        (int)detid.arm(), (int)detid.station(), (int)detid.plane(), (int)detid.channel()};
+    hists.leadingTime[chid] = iGetter.get("t_" + ch_name);
+    if (hists.leadingTime[chid] == nullptr) {
+      edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
+          << "Failed to retrieve leading time monitor for channel (" << detid << ").";
+      continue;
+    }
+    hists.toT[chid] = iGetter.get("tot_" + ch_name);
+    if (hists.toT[chid] == nullptr) {
+      edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
+          << "Failed to retrieve time over threshold monitor for channel (" << detid << ").";
+      continue;
+    }
+    hists.leadingTimeVsToT[chid] = iGetter.get("tvstot_" + ch_name);
+    if (hists.leadingTimeVsToT[chid] == nullptr) {
+      edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
+          << "Failed to retrieve leading time vs. time over threshold monitor for channel (" << detid << ").";
+      continue;
+    }
+    if (min_entries_ > 0 && hists.leadingTimeVsToT[chid]->getEntries() < min_entries_) {
+      edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
+          << "Not enough entries for channel (" << detid << "): " << hists.leadingTimeVsToT[chid]->getEntries() << " < "
+          << min_entries_ << ". Skipping calibration.";
+      continue;
+    }
+    const double upper_tot_range = hists.toT[chid]->getMean() + 2.5;
+    {  // scope for x-profile
+      std::unique_ptr<TProfile> prof(hists.leadingTimeVsToT[chid]->getTH2D()->ProfileX("_prof_x", 1, -1));
+      interp_.SetParameters(hists.leadingTime[chid]->getRMS(),
+                            hists.toT[chid]->getMean(),
+                            0.8,
+                            hists.leadingTime[chid]->getMean() - hists.leadingTime[chid]->getRMS());
+      const auto& res = prof->Fit(&interp_, "B+", "", 10.4, upper_tot_range);
+      if ((bool)res) {
+        calib_params[key] = {
+            interp_.GetParameter(0), interp_.GetParameter(1), interp_.GetParameter(2), interp_.GetParameter(3)};
+        calib_time[key] = std::make_pair(0.1, 0.);  // hardcoded resolution/offset placeholder for the time being
+        // can possibly do something with interp_.GetChiSquare() in the near future
+      } else
+        edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
+            << "Fit did not converge for channel (" << detid << ").";
+    }
+  }
+
+  // fill the DB object record
+  PPSTimingCalibration calib(formula_, calib_params, calib_time);
+
+  // write the object
+  edm::Service<cond::service::PoolDBOutputService> poolDbService;
+  poolDbService->writeOne(&calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd");
+}
+
+//------------------------------------------------------------------------------
+
+void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<std::string>("dqmDir", "AlCaReco/PPSTimingCalibrationPCL")
+      ->setComment("input path for the various DQM plots");
+  desc.add<std::string>("formula", "[0]/(exp((x-[1])/[2])+1)+[3]")
+      ->setComment("interpolation formula for the time walk component");
+  desc.add<unsigned int>("minEntries", 100)->setComment("minimal number of hits to extract calibration");
+  descriptions.addWithDefaultLabel(desc);
+}
+
+DEFINE_FWK_MODULE(PPSTimingCalibrationPCLHarvester);

CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc

@@ -0,0 +1,126 @@
+/****************************************************************************
+ *
+ * This is a part of PPS offline software.
+ * Authors:
+ *   Edoardo Bossini
+ *   Piotr Maciej Cwiklicki
+ *   Laurent Forthomme
+ *
+ ****************************************************************************/
+
+#include "DQMServices/Core/interface/DQMGlobalEDAnalyzer.h"
+#include "DQMServices/Core/interface/DQMStore.h"
+
+#include "FWCore/Framework/interface/Frameworkfwd.h"
+#include "FWCore/Framework/interface/MakerMacros.h"
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+
+#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
+#include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
+
+#include "DataFormats/Common/interface/DetSetVector.h"
+#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h"
+#include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h"
+
+#include "CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h"
+
+//------------------------------------------------------------------------------
+
+class PPSTimingCalibrationPCLWorker : public DQMGlobalEDAnalyzer<TimingCalibrationHistograms> {
+public:
+  explicit PPSTimingCalibrationPCLWorker(const edm::ParameterSet&);
+
+  void dqmAnalyze(const edm::Event&, const edm::EventSetup&, const TimingCalibrationHistograms&) const override;
+
+  static void fillDescriptions(edm::ConfigurationDescriptions&);
+
+private:
+  void bookHistograms(DQMStore::IBooker&,
+                      const edm::Run&,
+                      const edm::EventSetup&,
+                      TimingCalibrationHistograms&) const override;
+
+  edm::EDGetTokenT<edm::DetSetVector<CTPPSDiamondRecHit>> diamondRecHitToken_;
+  edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> geomEsToken_;
+
+  const std::string dqmDir_;
+};
+
+//------------------------------------------------------------------------------
+
+PPSTimingCalibrationPCLWorker::PPSTimingCalibrationPCLWorker(const edm::ParameterSet& iConfig)
+    : diamondRecHitToken_(
+          consumes<edm::DetSetVector<CTPPSDiamondRecHit>>(iConfig.getParameter<edm::InputTag>("diamondRecHitTag"))),
+      geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
+      dqmDir_(iConfig.getParameter<std::string>("dqmDir")) {}
+
+//------------------------------------------------------------------------------
+
+void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker,
+                                                   const edm::Run& iRun,
+                                                   const edm::EventSetup& iSetup,
+                                                   TimingCalibrationHistograms& iHists) const {
+  iBooker.cd();
+  iBooker.setCurrentFolder(dqmDir_);
+  std::string ch_name;
+
+  const auto& geom = iSetup.getData(geomEsToken_);
+  for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) {
+    if (!CTPPSDiamondDetId::check(it->first))
+      continue;
+    const CTPPSDiamondDetId detid(it->first);
+    if (detid.station() != 1)
+      continue;
+    detid.channelName(ch_name);
+    iHists.leadingTime[detid.rawId()] = iBooker.book1D("t_" + ch_name, ch_name + ";t (ns);Entries", 1200, -60., 60.);
+    iHists.toT[detid.rawId()] = iBooker.book1D("tot_" + ch_name, ch_name + ";ToT (ns);Entries", 100, -20., 20.);
+    iHists.leadingTimeVsToT[detid.rawId()] =
+        iBooker.book2D("tvstot_" + ch_name, ch_name + ";ToT (ns);t (ns)", 240, 0., 60., 450, -20., 25.);
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void PPSTimingCalibrationPCLWorker::dqmAnalyze(const edm::Event& iEvent,
+                                               const edm::EventSetup& iSetup,
+                                               const TimingCalibrationHistograms& iHists) const {
+  // then extract the rechits information for later processing
+  edm::Handle<edm::DetSetVector<CTPPSDiamondRecHit>> dsv_rechits;
+  iEvent.getByToken(diamondRecHitToken_, dsv_rechits);
+  // ensure timing detectors rechits are found in the event content
+  if (dsv_rechits->empty()) {
+    edm::LogWarning("PPSTimingCalibrationPCLWorker:dqmAnalyze") << "No rechits retrieved from the event content.";
+    return;
+  }
+  for (const auto& ds_rechits : *dsv_rechits) {
+    const CTPPSDiamondDetId detid(ds_rechits.detId());
+    if (iHists.leadingTimeVsToT.count(detid.rawId()) == 0) {
+      edm::LogWarning("PPSTimingCalibrationPCLWorker:dqmAnalyze")
+          << "Pad with detId=" << detid << " is not set to be monitored.";
+      continue;
+    }
+    for (const auto& rechit : ds_rechits) {
+      // skip invalid rechits
+      if (rechit.time() == 0. || rechit.toT() < 0.)
+        continue;
+      iHists.leadingTime.at(detid.rawId())->Fill(rechit.time());
+      iHists.toT.at(detid.rawId())->Fill(rechit.toT());
+      iHists.leadingTimeVsToT.at(detid.rawId())->Fill(rechit.toT(), rechit.time());
+    }
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void PPSTimingCalibrationPCLWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
+  edm::ParameterSetDescription desc;
+  desc.add<edm::InputTag>("diamondRecHitTag", edm::InputTag("ctppsDiamondUncalibRecHits"))
+      ->setComment("input tag for the PPS diamond detectors rechits");
+  desc.add<std::string>("dqmDir", "AlCaReco/PPSTimingCalibrationPCL")
+      ->setComment("output path for the various DQM plots");
+
+  descriptions.addWithDefaultLabel(desc);
+}
+
+DEFINE_FWK_MODULE(PPSTimingCalibrationPCLWorker);

CalibPPS/TimingCalibration/python/ALCARECOPPSCalTrackBasedSel_Output_cff.py

@@ -0,0 +1,15 @@
+import FWCore.ParameterSet.Config as cms
+
+OutALCARECOPPSCalTrackBasedSel_noDrop = cms.PSet(
+    SelectEvents = cms.untracked.PSet(
+        SelectEvents = cms.vstring('pathALCARECOPPSCalTrackBasedSel')
+    ),
+    outputCommands = cms.untracked.vstring(
+        'keep *_ALCARECOPPSCalTrackBasedSel_*_*',
+        'keep *_ctppsDiamondRawToDigi_*_*'
+        # will be updated to add the spatial alignment required collections
+    )
+)
+
+OutALCARECOPPSCalTrackBasedSel = OutALCARECOPPSCalTrackBasedSel_noDrop.clone()
+OutALCARECOPPSCalTrackBasedSel.outputCommands.insert(0, 'drop *')

CalibPPS/TimingCalibration/python/ALCARECOPPSCalTrackBasedSel_cff.py

@@ -0,0 +1,48 @@
+import FWCore.ParameterSet.Config as cms
+
+# define the HLT base path
+from HLTrigger.HLTfilters.hltHighLevel_cfi import hltHighLevel as _hlt
+ALCARECOPPSCalTrackBasedSelHLT = _hlt.clone(
+    andOr = True,
+    HLTPaths = ['HLT_ZeroBias_v*'],
+    #eventSetupPathKey = 'SiStripCalZeroBias', # in case we have a proper base key
+    throw = False
+)
+
+# perform basic PPS reconstruction
+from EventFilter.CTPPSRawToDigi.ctppsRawToDigi_cff import *
+from RecoPPS.Configuration.recoCTPPS_cff import *
+
+# select events passing the filter on pixel tracks
+from HLTrigger.special.hltPPSPerPotTrackFilter_cfi import hltPPSPerPotTrackFilter as _filter
+hltPPSPerPotTrackFilter = _filter.clone(
+    pixelFilter = cms.VPSet(
+        cms.PSet( # sector 45, near pot
+            detid = cms.uint32(2022703104),
+            minTracks = cms.int32(1),
+            maxTracks = cms.int32(6),
+        ),
+        cms.PSet( # sector 45, far pot
+            detid = cms.uint32(2023227392),
+            minTracks = cms.int32(1),
+            maxTracks = cms.int32(6),
+        ),
+        cms.PSet( # sector 56, near pot
+            detid = cms.uint32(2039480320),
+            minTracks = cms.int32(1),
+            maxTracks = cms.int32(6),
+        ),
+        cms.PSet( # sector 56, far pot
+            detid = cms.uint32(2040004608),
+            minTracks = cms.int32(1),
+            maxTracks = cms.int32(6),
+        ),
+    )
+)
+
+seqALCARECOPPSCalTrackBasedSel = cms.Sequence(
+    ctppsRawToDigi *
+    recoCTPPS *
+    ALCARECOPPSCalTrackBasedSelHLT *
+    hltPPSPerPotTrackFilter
+)

CalibPPS/TimingCalibration/python/ALCARECOPPSTimingCalib_Output_cff.py

@@ -0,0 +1,13 @@
+import FWCore.ParameterSet.Config as cms
+
+OutALCARECOPPSTimingCalib_noDrop = cms.PSet(
+    SelectEvents = cms.untracked.PSet(
+        SelectEvents = cms.vstring('pathALCARECOPPSTimingCalib')
+    ),
+    outputCommands = cms.untracked.vstring(
+        'keep *_MEtoEDMConvertPPSTimingCalib_*_*',
+    )
+)
+
+OutALCARECOPPSTimingCalib = OutALCARECOPPSTimingCalib_noDrop.clone()
+OutALCARECOPPSTimingCalib.outputCommands.insert(0, 'drop *')