TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc

0001 /****************************************************************************
0002  *
0003  * This is a part of PPS offline software.
0004  * Authors:
0005  *   Edoardo Bossini
0006  *   Piotr Maciej Cwiklicki
0007  *   Laurent Forthomme
0008  *
0009  ****************************************************************************/
0010
0011 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0012
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019
0020 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0021 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0022
0023 #include "CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h"
0024 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0025
0026 #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h"
0027 #include "CondFormats/PPSObjects/interface/PPSTimingCalibration.h"
0028 #include "TFitResult.h"
0029 //------------------------------------------------------------------------------
0030
0031 class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester {
0032 public:
0033   PPSTimingCalibrationPCLHarvester(const edm::ParameterSet&);
0034   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0035
0036   static void fillDescriptions(edm::ConfigurationDescriptions&);
0037
0038 private:
0039   void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override;
0040   edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> geomEsToken_;
0041   std::vector<CTPPSDiamondDetId> detids_;
0042   const std::string dqmDir_;
0043   const std::string formula_;
0044   const unsigned int min_entries_;
0045   static constexpr double resolution_ = 0.1;
0046   static constexpr double offset_ = 0.;
0047   TF1 interp_;
0048 };
0049
0050 //------------------------------------------------------------------------------
0051
0052 PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::ParameterSet& iConfig)
0053     : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0054       dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0055       formula_(iConfig.getParameter<std::string>("formula")),
0056       min_entries_(iConfig.getParameter<unsigned int>("minEntries")),
0057       interp_("interp", formula_.c_str(), 10.5, 25.) {
0058   // first ensure DB output service is available
0059   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0060   if (!poolDbService.isAvailable())
0061     throw cms::Exception("PPSTimingCalibrationPCLHarvester") << "PoolDBService required";
0062
0063   // constrain the min/max fit values
0064   interp_.SetParLimits(1, 9., 15.);
0065   interp_.SetParLimits(2, 0.2, 2.5);
0066 }
0067
0068 //------------------------------------------------------------------------------
0069
0070 void PPSTimingCalibrationPCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0071   const auto& geom = iSetup.getData(geomEsToken_);
0072   for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) {
0073     if (!CTPPSDiamondDetId::check(it->first))
0074       continue;
0075     const CTPPSDiamondDetId detid(it->first);
0076     detids_.emplace_back(detid);
0077   }
0078 }
0079
0080 //------------------------------------------------------------------------------
0081
0082 void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0083   // book the parameters containers
0084   PPSTimingCalibration::ParametersMap calib_params;
0085   PPSTimingCalibration::TimingMap calib_time;
0086
0087   iGetter.cd();
0088   iGetter.setCurrentFolder(dqmDir_);
0089
0090   // compute the fit parameters for all monitored channels
0091   TimingCalibrationHistograms hists;
0092   std::string ch_name;
0093   for (const auto& detid : detids_) {
0094     detid.channelName(ch_name);
0095     const auto chid = detid.rawId();
0096     const PPSTimingCalibration::Key key{
0097         (int)detid.arm(), (int)detid.station(), (int)detid.plane(), (int)detid.channel()};
0098
0099     calib_params[key] = {0, 0, 0, 0};
0100     calib_time[key] = std::make_pair(offset_, resolution_);
0101
0102     hists.leadingTime[chid] = iGetter.get(dqmDir_ + "/t_" + ch_name);
0103     if (hists.leadingTime[chid] == nullptr) {
0104       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0105           << "Failed to retrieve leading time monitor for channel (" << detid << ").";
0106       continue;
0107     }
0108     hists.toT[chid] = iGetter.get(dqmDir_ + "/tot_" + ch_name);
0109     if (hists.toT[chid] == nullptr) {
0110       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0111           << "Failed to retrieve time over threshold monitor for channel (" << detid << ").";
0112       continue;
0113     }
0114     hists.leadingTimeVsToT[chid] = iGetter.get(dqmDir_ + "/tvstot_" + ch_name);
0115     if (hists.leadingTimeVsToT[chid] == nullptr) {
0116       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0117           << "Failed to retrieve leading time vs. time over threshold monitor for channel (" << detid << ").";
0118       continue;
0119     }
0120     if (min_entries_ > 0 && hists.leadingTimeVsToT[chid]->getEntries() < min_entries_) {
0121       edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0122           << "Not enough entries for channel (" << detid << "): " << hists.leadingTimeVsToT[chid]->getEntries() << " < "
0123           << min_entries_ << ". Skipping calibration.";
0124       continue;
0125     }
0126     const double upper_tot_range = hists.toT[chid]->getMean() + 2.5;
0127     {  // scope for x-profile
0128
0129       std::string ch_name;
0130       detid.channelName(ch_name);
0131       auto profile = iBooker.bookProfile(ch_name + "_prof_x", ch_name + "_prof_x", 240, 0., 60., 450, -20., 25.);
0132
0133       std::unique_ptr<TProfile> prof(hists.leadingTimeVsToT[chid]->getTH2F()->ProfileX("_prof_x", 1, -1));
0134       *(profile->getTProfile()) = *((TProfile*)prof->Clone());
0135       profile->getTProfile()->SetTitle(ch_name.c_str());
0136       profile->getTProfile()->SetName(ch_name.c_str());
0137
0138       interp_.SetParameters(hists.leadingTime[chid]->getRMS(),
0139                             hists.toT[chid]->getMean(),
0140                             0.8,
0141                             hists.leadingTime[chid]->getMean() - hists.leadingTime[chid]->getRMS());
0142       const auto& res = profile->getTProfile()->Fit(&interp_, "B+", "", 10.4, upper_tot_range);
0143       if (!(bool)res) {
0144         calib_params[key] = {
0145             interp_.GetParameter(0), interp_.GetParameter(1), interp_.GetParameter(2), interp_.GetParameter(3)};
0146         calib_time[key] =
0147             std::make_pair(offset_, resolution_);  // hardcoded offset/resolution placeholder for the time being
0148         // can possibly do something with interp_.GetChiSquare() in the near future
0149
0150       } else
0151         edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0152             << "Fit did not converge for channel (" << detid << ").";
0153     }
0154   }
0155
0156   // fill the DB object record
0157   PPSTimingCalibration calib(formula_, calib_params, calib_time);
0158
0159   // write the object
0160   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0161   poolDbService->writeOneIOV(calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd_HPTDC");
0162 }
0163
0164 //------------------------------------------------------------------------------
0165
0166 void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0167   edm::ParameterSetDescription desc;
0168   desc.add<std::string>("dqmDir", "AlCaReco/PPSTimingCalibrationPCL")
0169       ->setComment("input path for the various DQM plots");
0170   desc.add<std::string>("formula", "[0]/(exp((x-[1])/[2])+1)+[3]")
0171       ->setComment("interpolation formula for the time walk component");
0172   desc.add<unsigned int>("minEntries", 100)->setComment("minimal number of hits to extract calibration");
0173   descriptions.addWithDefaultLabel(desc);
0174 }
0175
0176 DEFINE_FWK_MODULE(PPSTimingCalibrationPCLHarvester);