File indexing completed on 2024-04-06 11:58:37
0001
0002
0003
0004
0005
0006
0007
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 const double threshold_fraction_of_max_;
0046 static constexpr double upper_limit_max_search_ = 20;
0047 static constexpr double upper_limit_range_search_ = 20;
0048 static constexpr double lower_limit_range_search_ = 8;
0049 static constexpr double resolution_ = 0.1;
0050 static constexpr double offset_ = 0.;
0051 TF1 interp_;
0052 };
0053
0054
0055
0056 PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::ParameterSet& iConfig)
0057 : geomEsToken_(esConsumes<edm::Transition::BeginRun>()),
0058 dqmDir_(iConfig.getParameter<std::string>("dqmDir")),
0059 formula_(iConfig.getParameter<std::string>("formula")),
0060 min_entries_(iConfig.getParameter<unsigned int>("minEntries")),
0061 threshold_fraction_of_max_(iConfig.getParameter<double>("thresholdFractionOfMax")),
0062 interp_("interp", formula_.c_str(), 10.5, 25.) {
0063
0064 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0065 if (!poolDbService.isAvailable())
0066 throw cms::Exception("PPSTimingCalibrationPCLHarvester") << "PoolDBService required";
0067
0068
0069 interp_.SetParLimits(1, 9., 15.);
0070 interp_.SetParLimits(2, 0.2, 2.5);
0071 }
0072
0073
0074
0075 void PPSTimingCalibrationPCLHarvester::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0076 const auto& geom = iSetup.getData(geomEsToken_);
0077 for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) {
0078 if (!CTPPSDiamondDetId::check(it->first))
0079 continue;
0080 const CTPPSDiamondDetId detid(it->first);
0081 detids_.emplace_back(detid);
0082 }
0083 }
0084
0085
0086
0087 void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {
0088
0089 PPSTimingCalibration::ParametersMap calib_params;
0090 PPSTimingCalibration::TimingMap calib_time;
0091
0092 iGetter.cd();
0093 iGetter.setCurrentFolder(dqmDir_);
0094
0095
0096 TimingCalibrationHistograms hists;
0097 std::string ch_name;
0098 for (const auto& detid : detids_) {
0099 detid.channelName(ch_name);
0100 const auto chid = detid.rawId();
0101 const PPSTimingCalibration::Key key{
0102 (int)detid.arm(), (int)detid.station(), (int)detid.plane(), (int)detid.channel()};
0103
0104 calib_params[key] = {0, 0, 0, 0};
0105 calib_time[key] = std::make_pair(offset_, resolution_);
0106
0107 hists.leadingTime[chid] = iGetter.get(dqmDir_ + "/t_" + ch_name);
0108 if (hists.leadingTime[chid] == nullptr) {
0109 edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0110 << "Failed to retrieve leading time monitor for channel (" << detid << ").";
0111 continue;
0112 }
0113 hists.toT[chid] = iGetter.get(dqmDir_ + "/tot_" + ch_name);
0114 if (hists.toT[chid] == nullptr) {
0115 edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0116 << "Failed to retrieve time over threshold monitor for channel (" << detid << ").";
0117 continue;
0118 }
0119 hists.leadingTimeVsToT[chid] = iGetter.get(dqmDir_ + "/tvstot_" + ch_name);
0120 if (hists.leadingTimeVsToT[chid] == nullptr) {
0121 edm::LogInfo("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0122 << "Failed to retrieve leading time vs. time over threshold monitor for channel (" << detid << ").";
0123 continue;
0124 }
0125 if (min_entries_ > 0 && hists.leadingTimeVsToT[chid]->getEntries() < min_entries_) {
0126 edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0127 << "Not enough entries for channel (" << detid << "): " << hists.leadingTimeVsToT[chid]->getEntries() << " < "
0128 << min_entries_ << ". Skipping calibration.";
0129 continue;
0130 }
0131
0132
0133 int max_bin_pos = 1;
0134 for (int i = 0; i < hists.toT[chid]->getNbinsX(); i++) {
0135 double bin_value = hists.toT[chid]->getBinContent(i);
0136 int bin_x_pos = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(i);
0137 if (bin_x_pos > upper_limit_max_search_)
0138 break;
0139 if (bin_value > hists.toT[chid]->getBinContent(max_bin_pos))
0140 max_bin_pos = i;
0141 }
0142
0143 int upper_limit_pos = max_bin_pos;
0144 int lower_limit_pos = max_bin_pos;
0145 double threshold = threshold_fraction_of_max_ * hists.toT[chid]->getBinContent(max_bin_pos);
0146 while (hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos) < upper_limit_range_search_) {
0147 upper_limit_pos++;
0148 if (hists.toT[chid]->getBinContent(upper_limit_pos) < threshold)
0149 break;
0150 }
0151 while (hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos) > lower_limit_range_search_) {
0152 lower_limit_pos--;
0153 if (hists.toT[chid]->getBinContent(lower_limit_pos) < threshold)
0154 break;
0155 }
0156 double upper_tot_range = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos);
0157 double lower_tot_range = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos);
0158
0159 {
0160
0161 std::string ch_name;
0162 detid.channelName(ch_name);
0163 auto profile = iBooker.bookProfile(ch_name + "_prof_x", ch_name + "_prof_x", 240, 0., 60., 450, -20., 25.);
0164
0165 std::unique_ptr<TProfile> prof(hists.leadingTimeVsToT[chid]->getTH2F()->ProfileX("_prof_x", 1, -1));
0166 *(profile->getTProfile()) = *((TProfile*)prof->Clone());
0167 profile->getTProfile()->SetTitle(ch_name.c_str());
0168 profile->getTProfile()->SetName(ch_name.c_str());
0169
0170 interp_.SetParameters(hists.leadingTime[chid]->getRMS(),
0171 hists.toT[chid]->getMean(),
0172 0.8,
0173 hists.leadingTime[chid]->getMean() - hists.leadingTime[chid]->getRMS());
0174 const auto& res = profile->getTProfile()->Fit(&interp_, "B+", "", lower_tot_range, upper_tot_range);
0175 if (!(bool)res) {
0176 calib_params[key] = {
0177 interp_.GetParameter(0), interp_.GetParameter(1), interp_.GetParameter(2), interp_.GetParameter(3)};
0178 calib_time[key] =
0179 std::make_pair(offset_, resolution_);
0180
0181
0182 } else
0183 edm::LogWarning("PPSTimingCalibrationPCLHarvester:dqmEndJob")
0184 << "Fit did not converge for channel (" << detid << ").";
0185 }
0186 }
0187
0188
0189 PPSTimingCalibration calib(formula_, calib_params, calib_time);
0190
0191
0192 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0193 poolDbService->writeOneIOV(calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd_HPTDC");
0194 }
0195
0196
0197
0198 void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0199 edm::ParameterSetDescription desc;
0200 desc.add<std::string>("dqmDir", "AlCaReco/PPSTimingCalibrationPCL")
0201 ->setComment("input path for the various DQM plots");
0202 desc.add<std::string>("formula", "[0]/(exp((x-[1])/[2])+1)+[3]")
0203 ->setComment("interpolation formula for the time walk component");
0204 desc.add<unsigned int>("minEntries", 100)->setComment("minimal number of hits to extract calibration");
0205 desc.add<double>("thresholdFractionOfMax", 0.05)
0206 ->setComment("threshold for TOT fit, defined as percent of max count in 1D TOT");
0207 descriptions.addWithDefaultLabel(desc);
0208 }
0209
0210 DEFINE_FWK_MODULE(PPSTimingCalibrationPCLHarvester);
