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     if (detid.station() == 1)  // for the time being, only compute for this station (run 2 diamond)
0077       detids_.emplace_back(detid);
0078   }
0079 }
0080
0081 //------------------------------------------------------------------------------
0082
0083 void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0084   // book the parameters containers
0085   PPSTimingCalibration::ParametersMap calib_params;
0086   PPSTimingCalibration::TimingMap calib_time;
0087
0088   iGetter.cd();
0089   iGetter.setCurrentFolder(dqmDir_);
0090
0091   // compute the fit parameters for all monitored channels
0092   TimingCalibrationHistograms hists;
0093   std::string ch_name;
0094   for (const auto& detid : detids_) {
0095     detid.channelName(ch_name);
0096     const auto chid = detid.rawId();
0097     const PPSTimingCalibration::Key key{
0098         (int)detid.arm(), (int)detid.station(), (int)detid.plane(), (int)detid.channel()};
0099
0100     calib_params[key] = {0, 0, 0, 0};
0101     calib_time[key] = std::make_pair(offset_, resolution_);
0102
0103     hists.leadingTime[chid] = iGetter.get(dqmDir_ + "/t_" + ch_name);
0104     if (hists.leadingTime[chid] == nullptr) {
0105       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0106           << "Failed to retrieve leading time monitor for channel (" << detid << ").";
0107       continue;
0108     }
0109     hists.toT[chid] = iGetter.get(dqmDir_ + "/tot_" + ch_name);
0110     if (hists.toT[chid] == nullptr) {
0111       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0112           << "Failed to retrieve time over threshold monitor for channel (" << detid << ").";
0113       continue;
0114     }
0115     hists.leadingTimeVsToT[chid] = iGetter.get(dqmDir_ + "/tvstot_" + ch_name);
0116     if (hists.leadingTimeVsToT[chid] == nullptr) {
0117       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0118           << "Failed to retrieve leading time vs. time over threshold monitor for channel (" << detid << ").";
0119       continue;
0120     }
0121     if (min_entries_ > 0 && hists.leadingTimeVsToT[chid]->getEntries() < min_entries_) {
0122       edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0123           << "Not enough entries for channel (" << detid << "): " << hists.leadingTimeVsToT[chid]->getEntries() << " < "
0124           << min_entries_ << ". Skipping calibration.";
0125       continue;
0126     }
0127     const double upper_tot_range = hists.toT[chid]->getMean() + 2.5;
0128     {  // scope for x-profile
0129
0130       std::string ch_name;
0131       detid.channelName(ch_name);
0132       auto profile = iBooker.bookProfile(ch_name + "_prof_x", ch_name + "_prof_x", 240, 0., 60., 450, -20., 25.);
0133
0134       std::unique_ptr<TProfile> prof(hists.leadingTimeVsToT[chid]->getTH2F()->ProfileX("_prof_x", 1, -1));
0135       *(profile->getTProfile()) = *((TProfile*)prof->Clone());
0136       profile->getTProfile()->SetTitle(ch_name.c_str());
0137       profile->getTProfile()->SetName(ch_name.c_str());
0138
0139       interp_.SetParameters(hists.leadingTime[chid]->getRMS(),
0140                             hists.toT[chid]->getMean(),
0141                             0.8,
0142                             hists.leadingTime[chid]->getMean() - hists.leadingTime[chid]->getRMS());
0143       const auto& res = profile->getTProfile()->Fit(&interp_, "B+", "", 10.4, upper_tot_range);
0144       if (!(bool)res) {
0145         calib_params[key] = {
0146             interp_.GetParameter(0), interp_.GetParameter(1), interp_.GetParameter(2), interp_.GetParameter(3)};
0147         calib_time[key] =
0148             std::make_pair(offset_, resolution_);  // hardcoded offset/resolution placeholder for the time being
0149         // can possibly do something with interp_.GetChiSquare() in the near future
0150
0151       } else
0152         edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0153             << "Fit did not converge for channel (" << detid << ").";
0154     }
0155   }
0156
0157   // fill the DB object record
0158   PPSTimingCalibration calib(formula_, calib_params, calib_time);
0159
0160   // write the object
0161   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0162   poolDbService->writeOneIOV(calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd_HPTDC");
0163 }
0164
0165 //------------------------------------------------------------------------------
0166
0167 void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0168   edm::ParameterSetDescription desc;
0169   desc.add<std::string>("dqmDir", "AlCaReco/PPSTimingCalibrationPCL")
0170       ->setComment("input path for the various DQM plots");
0171   desc.add<std::string>("formula", "[0]/(exp((x-[1])/[2])+1)+[3]")
0172       ->setComment("interpolation formula for the time walk component");
0173   desc.add<unsigned int>("minEntries", 100)->setComment("minimal number of hits to extract calibration");
0174   descriptions.addWithDefaultLabel(desc);
0175 }
0176
0177 DEFINE_FWK_MODULE(PPSTimingCalibrationPCLHarvester);