Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-16 23:26:25

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 
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   TF1 interp_;
0046 };
0047 
0048 //------------------------------------------------------------------------------
0049 
0050 PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::ParameterSet& iConfig)
0051     : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0052       dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0053       formula_(iConfig.getParameter<std::string>("formula")),
0054       min_entries_(iConfig.getParameter<unsigned int>("minEntries")),
0055       interp_("interp", formula_.c_str(), 10.5, 25.) {
0056   // first ensure DB output service is available
0057   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0058   if (!poolDbService.isAvailable())
0059     throw cms::Exception("PPSTimingCalibrationPCLHarvester") << "PoolDBService required";
0060 
0061   // constrain the min/max fit values
0062   interp_.SetParLimits(1, 9., 15.);
0063   interp_.SetParLimits(2, 0.2, 2.5);
0064 }
0065 
0066 //------------------------------------------------------------------------------
0067 
0068 void PPSTimingCalibrationPCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0069   const auto& geom = iSetup.getData(geomEsToken_);
0070   for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) {
0071     if (!CTPPSDiamondDetId::check(it->first))
0072       continue;
0073     const CTPPSDiamondDetId detid(it->first);
0074     if (detid.station() == 1)  // for the time being, only compute for this station (run 2 diamond)
0075       detids_.emplace_back(detid);
0076   }
0077 }
0078 
0079 //------------------------------------------------------------------------------
0080 
0081 void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0082   // book the parameters containers
0083   PPSTimingCalibration::ParametersMap calib_params;
0084   PPSTimingCalibration::TimingMap calib_time;
0085 
0086   iGetter.cd();
0087   iGetter.setCurrentFolder(dqmDir_);
0088 
0089   // compute the fit parameters for all monitored channels
0090   TimingCalibrationHistograms hists;
0091   std::string ch_name;
0092   for (const auto& detid : detids_) {
0093     detid.channelName(ch_name);
0094     const auto chid = detid.rawId();
0095     const PPSTimingCalibration::Key key{
0096         (int)detid.arm(), (int)detid.station(), (int)detid.plane(), (int)detid.channel()};
0097     hists.leadingTime[chid] = iGetter.get(dqmDir_ + "/t_" + ch_name);
0098     if (hists.leadingTime[chid] == nullptr) {
0099       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0100           << "Failed to retrieve leading time monitor for channel (" << detid << ").";
0101       continue;
0102     }
0103     hists.toT[chid] = iGetter.get(dqmDir_ + "/tot_" + ch_name);
0104     if (hists.toT[chid] == nullptr) {
0105       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0106           << "Failed to retrieve time over threshold monitor for channel (" << detid << ").";
0107       continue;
0108     }
0109     hists.leadingTimeVsToT[chid] = iGetter.get(dqmDir_ + "/tvstot_" + ch_name);
0110     if (hists.leadingTimeVsToT[chid] == nullptr) {
0111       edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0112           << "Failed to retrieve leading time vs. time over threshold monitor for channel (" << detid << ").";
0113       continue;
0114     }
0115     if (min_entries_ > 0 && hists.leadingTimeVsToT[chid]->getEntries() < min_entries_) {
0116       edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0117           << "Not enough entries for channel (" << detid << "): " << hists.leadingTimeVsToT[chid]->getEntries() << " < "
0118           << min_entries_ << ". Skipping calibration.";
0119       continue;
0120     }
0121     const double upper_tot_range = hists.toT[chid]->getMean() + 2.5;
0122     {  // scope for x-profile
0123       std::unique_ptr<TProfile> prof(hists.leadingTimeVsToT[chid]->getTH2F()->ProfileX("_prof_x", 1, -1));
0124       interp_.SetParameters(hists.leadingTime[chid]->getRMS(),
0125                             hists.toT[chid]->getMean(),
0126                             0.8,
0127                             hists.leadingTime[chid]->getMean() - hists.leadingTime[chid]->getRMS());
0128       const auto& res = prof->Fit(&interp_, "B+", "", 10.4, upper_tot_range);
0129       if ((bool)res) {
0130         calib_params[key] = {
0131             interp_.GetParameter(0), interp_.GetParameter(1), interp_.GetParameter(2), interp_.GetParameter(3)};
0132         calib_time[key] = std::make_pair(0.1, 0.);  // hardcoded resolution/offset placeholder for the time being
0133         // can possibly do something with interp_.GetChiSquare() in the near future
0134       } else
0135         edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0136             << "Fit did not converge for channel (" << detid << ").";
0137     }
0138   }
0139 
0140   // fill the DB object record
0141   PPSTimingCalibration calib(formula_, calib_params, calib_time);
0142 
0143   // write the object
0144   edm::Service<cond::service::PoolDBOutputService> poolDbService;
0145   poolDbService->writeOneIOV(calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd");
0146 }
0147 
0148 //------------------------------------------------------------------------------
0149 
0150 void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0151   edm::ParameterSetDescription desc;
0152   desc.add<std::string>("dqmDir", "AlCaReco/PPSTimingCalibrationPCL")
0153       ->setComment("input path for the various DQM plots");
0154   desc.add<std::string>("formula", "[0]/(exp((x-[1])/[2])+1)+[3]")
0155       ->setComment("interpolation formula for the time walk component");
0156   desc.add<unsigned int>("minEntries", 100)->setComment("minimal number of hits to extract calibration");
0157   descriptions.addWithDefaultLabel(desc);
0158 }
0159 
0160 DEFINE_FWK_MODULE(PPSTimingCalibrationPCLHarvester);