diff --git a/CalibPPS/TimingCalibration/BuildFile.xml b/CalibPPS/TimingCalibration/BuildFile.xml new file mode 100644 index 0000000000000..feb1ce7e23f4a --- /dev/null +++ b/CalibPPS/TimingCalibration/BuildFile.xml @@ -0,0 +1,6 @@ + + + + + + diff --git a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h new file mode 100644 index 0000000000000..c230bec63b193 --- /dev/null +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -0,0 +1,34 @@ +#ifndef CalibPPS_TimingCalibration_DoublePeakCorrection_h +#define CalibPPS_TimingCalibration_DoublePeakCorrection_h + +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" + +#include "TFile.h" +#include "TH2F.h" + +#include + +class DoublePeakCorrection { +public: + void extractLsAndTimeOffset(const std::string&, const unsigned int, const std::vector&); + void fillLsAndTimeOffset(const TH2F*, const PlaneKey&); + bool isCorrectionNeeded(const PlaneKey&) const; + double getCorrectedLeadingTime(const double, const unsigned int, const PlaneKey&) const; + double getEncodedLsAndTimeOffset(const PlaneKey&) const; + + static double GetCorrectedLeadingTime(const double, const unsigned int, const double); + +private: + const TH2F* getTVsLs(TFile&, const std::string&, const CTPPSDiamondDetId&); + std::tuple findLsAndTimePeaks(const TH2F*, const PlaneKey&) const; + double findTimeOffset(const TH2F*, const unsigned int, const double, const double) const; + double findGaussianMean(const std::unique_ptr&, const double) const; + + std::unordered_map, PlaneKeyHash> lsAndTimeOffsets_; + + constexpr static unsigned int LsEncodingMultiple_{100}; +}; + +#endif diff --git a/CalibPPS/TimingCalibration/interface/PlaneMap.h b/CalibPPS/TimingCalibration/interface/PlaneMap.h new file mode 100644 index 0000000000000..5b0c28e3719b0 --- /dev/null +++ b/CalibPPS/TimingCalibration/interface/PlaneMap.h @@ -0,0 +1,17 @@ +#ifndef CalibPPS_TimingCalibration_PlaneMap_h +#define CalibPPS_TimingCalibration_PlaneMap_h + +#include +#include +#include + +using PlaneKey = std::tuple; + +struct PlaneKeyHash { + std::size_t operator()(const PlaneKey& planeKey) const noexcept { + return std::hash()(std::get<0>(planeKey)) ^ std::hash()(std::get<1>(planeKey)) ^ + std::hash()(std::get<2>(planeKey)); + } +}; + +#endif diff --git a/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h new file mode 100644 index 0000000000000..9cfc57bc18413 --- /dev/null +++ b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h @@ -0,0 +1,21 @@ +#ifndef CalibPPS_TimingCalibration_TimingCalibrationData_h +#define CalibPPS_TimingCalibration_TimingCalibrationData_h + +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" + +#include "DQMServices/Core/interface/MonitorElement.h" + +struct TimingCalibrationData { + using PlaneMonitorMap = std::unordered_map; + using ChannelMonitorMap = std::unordered_map; + + ChannelMonitorMap leadingTime; + ChannelMonitorMap toT; + ChannelMonitorMap leadingTimeVsToT; + PlaneMonitorMap leadingTimeVsLs; + + DoublePeakCorrection doublePeakCorrection; +}; + +#endif diff --git a/CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h b/CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h deleted file mode 100644 index 2d847dcdfe6f4..0000000000000 --- a/CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h +++ /dev/null @@ -1,17 +0,0 @@ -#ifndef CalibPPS_TimingCalibration_TimingCalibrationStruct_h -#define CalibPPS_TimingCalibration_TimingCalibrationStruct_h - -#include "DQMServices/Core/interface/DQMStore.h" -#include - -struct TimingCalibrationHistograms { -public: - TimingCalibrationHistograms() = default; - - using MonitorMap = std::unordered_map; - - MonitorMap leadingTime, toT; - MonitorMap leadingTimeVsToT; -}; - -#endif diff --git a/CalibPPS/TimingCalibration/plugins/BuildFile.xml b/CalibPPS/TimingCalibration/plugins/BuildFile.xml index 83f19f4dc0ff3..c59fd1bfa2bb0 100644 --- a/CalibPPS/TimingCalibration/plugins/BuildFile.xml +++ b/CalibPPS/TimingCalibration/plugins/BuildFile.xml @@ -1,10 +1,11 @@ - - + + + + - diff --git a/CalibPPS/TimingCalibration/plugins/PPSDiamondSampicTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSDiamondSampicTimingCalibrationPCLHarvester.cc index 291c5bcf5dbed..eff84a774d304 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSDiamondSampicTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSDiamondSampicTimingCalibrationPCLHarvester.cc @@ -33,7 +33,6 @@ #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" diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index c75e2403a7568..5c6f320db455e 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -5,69 +5,76 @@ * Edoardo Bossini * Piotr Maciej Cwiklicki * Laurent Forthomme + * Tomasz Ostafin * ****************************************************************************/ +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" + +#include "CondCore/DBOutputService/interface/PoolDBOutputService.h" + +#include "CondFormats/PPSObjects/interface/PPSTimingCalibration.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" + #include "DQMServices/Core/interface/DQMEDHarvester.h" +#include "FWCore/Framework/interface/EventSetup.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 "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.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" +#include "Math/MinimizerOptions.h" #include "TFitResult.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 geomEsToken_; - std::vector detids_; + + bool fetchWorkerHistograms( + DQMStore::IGetter&, const CTPPSDiamondDetId&, const PlaneKey&, const uint32_t, const std::string&); + TProfile* createTVsTotProfile(DQMStore::IBooker&, dqm::reco::MonitorElement*, const std::string&) const; + std::pair findFitRange(const TH1F*, const int, const double, const double) const; + + TimingCalibrationData timingCalibartionData_; const std::string dqmDir_; const std::string formula_; - const unsigned int min_entries_; - const double threshold_fraction_of_max_; - static constexpr double upper_limit_max_search_ = 20; - static constexpr double upper_limit_range_search_ = 20; - static constexpr double lower_limit_range_search_ = 8; - static constexpr double resolution_ = 0.1; - static constexpr double offset_ = 0.; - TF1 interp_; + const std::string tVsLsFilename_; + std::vector detIds_; + const edm::ESGetToken geomEsToken_; + const unsigned int minEntries_; + + static constexpr double FixedFitBoundIndication_{-1.0}; }; //------------------------------------------------------------------------------ PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::ParameterSet& iConfig) - : geomEsToken_(esConsumes()), - dqmDir_(iConfig.getParameter("dqmDir")), - formula_(iConfig.getParameter("formula")), - min_entries_(iConfig.getParameter("minEntries")), - threshold_fraction_of_max_(iConfig.getParameter("thresholdFractionOfMax")), - interp_("interp", formula_.c_str(), 10.5, 25.) { + : dqmDir_{iConfig.getParameter("dqmDir")}, + formula_{iConfig.getParameter("formula")}, + tVsLsFilename_{iConfig.getParameter("tVsLsFilename")}, + geomEsToken_{esConsumes()}, + minEntries_{iConfig.getParameter("minEntries")} { // first ensure DB output service is available edm::Service 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); + if (!poolDbService.isAvailable()) { + throw cms::Exception{"PPSTimingCalibrationPCLHarvester"} << "PoolDBService required"; + } } //------------------------------------------------------------------------------ @@ -75,122 +82,217 @@ PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::Pa 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); - detids_.emplace_back(detid); + if (CTPPSDiamondDetId::check(it->first)) { + const CTPPSDiamondDetId detId{it->first}; + detIds_.push_back(detId); + } } + timingCalibartionData_.doublePeakCorrection.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds_); } //------------------------------------------------------------------------------ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) { // book the parameters containers - PPSTimingCalibration::ParametersMap calib_params; - PPSTimingCalibration::TimingMap calib_time; + PPSTimingCalibration::ParametersMap calibParams; + PPSTimingCalibration::TimingMap calibTime; + + TF1 interp{"interp", formula_.c_str()}; + interp.SetParLimits(0, 0.5, 5.0); + interp.SetParLimits(1, 4.0, 15.0); + interp.SetParLimits(2, 0.1, 4.0); + interp.SetParLimits(3, 0.1, 15.0); + + // set a higher max function calls limit for the ROOT fit algorithm + ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(5'000); iGetter.cd(); iGetter.setCurrentFolder(dqmDir_); + constexpr double defaultOffset{0.0}; + constexpr double defaultResolution{0.1}; + constexpr std::array thresholds{ + {FixedFitBoundIndication_, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07}}; + // 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()}; - - calib_params[key] = {0, 0, 0, 0}; - calib_time[key] = std::make_pair(offset_, resolution_); - - hists.leadingTime[chid] = iGetter.get(dqmDir_ + "/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(dqmDir_ + "/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(dqmDir_ + "/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; - } + std::string channelName; + for (const auto& detId : detIds_) { + const uint32_t channelId{detId.rawId()}; + detId.channelName(channelName); + const PlaneKey planeKey{detId.arm(), detId.station(), detId.plane()}; + const PPSTimingCalibration::Key channelKey{static_cast(detId.arm()), + static_cast(detId.station()), + static_cast(detId.plane()), + static_cast(detId.channel())}; + calibParams[channelKey] = {0.0, 0.0, 0.0, 0.0}; + calibTime[channelKey] = {defaultOffset, defaultResolution}; + if (fetchWorkerHistograms(iGetter, detId, planeKey, channelId, channelName)) { + timingCalibartionData_.doublePeakCorrection.fillLsAndTimeOffset( + timingCalibartionData_.leadingTimeVsLs[planeKey]->getTH2F(), planeKey); + if (tVsLsFilename_.empty() && timingCalibartionData_.doublePeakCorrection.isCorrectionNeeded(planeKey)) { + calibTime[channelKey] = {0.0, 0.0}; + } else { + TProfile* const tVsTotProfile{ + createTVsTotProfile(iBooker, timingCalibartionData_.leadingTimeVsToT.at(channelId), channelName)}; - //find max - int max_bin_pos = 1; - for (int i = 0; i < hists.toT[chid]->getNbinsX(); i++) { - double bin_value = hists.toT[chid]->getBinContent(i); - int bin_x_pos = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(i); - if (bin_x_pos > upper_limit_max_search_) - break; - if (bin_value > hists.toT[chid]->getBinContent(max_bin_pos)) - max_bin_pos = i; - } - //find ranges - int upper_limit_pos = max_bin_pos; - int lower_limit_pos = max_bin_pos; - double threshold = threshold_fraction_of_max_ * hists.toT[chid]->getBinContent(max_bin_pos); - while (hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos) < upper_limit_range_search_) { - upper_limit_pos++; - if (hists.toT[chid]->getBinContent(upper_limit_pos) < threshold) - break; + const TH1F* tot{timingCalibartionData_.toT.at(channelId)->getTH1F()}; + const int maxTotBin{tot->GetMaximumBin()}; + + const double defaultUpperLowerAsymptotesDiff{timingCalibartionData_.leadingTime.at(channelId)->getRMS()}; + const double defaultCenterOfDistribution{tot->GetMean()}; + const double defaultLowerAsymptote{timingCalibartionData_.leadingTime.at(channelId)->getMean() - defaultUpperLowerAsymptotesDiff}; + + double bestChiSqDivNdf{std::numeric_limits::max()}; + double bestLowerTotRange{0.0}; + double bestUpperTotRange{0.0}; + for (const double upperThresholdFractionOfMax : thresholds) { + for (const double lowerThresholdFractionOfMax : thresholds) { + constexpr double defaultFitSlope{0.8}; + interp.SetParameters( + defaultUpperLowerAsymptotesDiff, defaultCenterOfDistribution, defaultFitSlope, defaultLowerAsymptote); + const auto [lowerTotRange, upperTotRange] = + findFitRange(tot, maxTotBin, lowerThresholdFractionOfMax, upperThresholdFractionOfMax); + + const TFitResultPtr& tVsTotProfileFitResult{ + tVsTotProfile->Fit(&interp, "BNS", "", lowerTotRange, upperTotRange)}; + if (tVsTotProfileFitResult->IsValid()) { + const double chiSqDivNdf{tVsTotProfileFitResult->Chi2() / tVsTotProfileFitResult->Ndf()}; + if (chiSqDivNdf < bestChiSqDivNdf) { + bestChiSqDivNdf = chiSqDivNdf; + bestUpperTotRange = upperTotRange; + bestLowerTotRange = lowerTotRange; + } + } + } + } + + if (bestUpperTotRange != 0.0) { + tVsTotProfile->Fit(&interp, "B", "", bestLowerTotRange, bestUpperTotRange); + calibParams[channelKey] = { + interp.GetParameter(0), interp.GetParameter(1), interp.GetParameter(2), interp.GetParameter(3)}; + calibTime[channelKey] = {timingCalibartionData_.doublePeakCorrection.getEncodedLsAndTimeOffset(planeKey), + defaultResolution}; + } else { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:dqmEndJob"} << "Fit did not converge for channel (" << detId + << ")."; + } + } } - while (hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos) > lower_limit_range_search_) { - lower_limit_pos--; - if (hists.toT[chid]->getBinContent(lower_limit_pos) < threshold) - break; + + // fill the DB object record + PPSTimingCalibration calib{formula_, calibParams, calibTime}; + + // write the object + edm::Service poolDbService; + poolDbService->writeOneIOV(calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd_HPTDC"); + } +} + +//------------------------------------------------------------------------------ + +bool PPSTimingCalibrationPCLHarvester::fetchWorkerHistograms(DQMStore::IGetter& iGetter, + const CTPPSDiamondDetId& detId, + const PlaneKey& planeKey, + const uint32_t channelId, + const std::string& channelName) { + timingCalibartionData_.leadingTime[channelId] = iGetter.get(dqmDir_ + "/t_" + channelName); + if (!timingCalibartionData_.leadingTime.at(channelId)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve leading time monitor for channel (" << detId << "). Skipping calibration."; + return false; + } + + timingCalibartionData_.toT[channelId] = iGetter.get(dqmDir_ + "/tot_" + channelName); + if (!timingCalibartionData_.toT.at(channelId)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve time over threshold monitor for channel (" << detId << "). Skipping calibration."; + return false; + } + + timingCalibartionData_.leadingTimeVsToT[channelId] = iGetter.get(dqmDir_ + "/tvstot_" + channelName); + if (!timingCalibartionData_.leadingTimeVsToT.at(channelId)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve leading time vs. time over threshold monitor for channel (" << detId + << "). Skipping calibration."; + return false; + } + + const auto tVsTotEntries = + static_cast(timingCalibartionData_.leadingTimeVsToT.at(channelId)->getEntries()); + if (tVsTotEntries < minEntries_) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Not enough entries for channel (" << detId << "): " << tVsTotEntries << " < " << minEntries_ + << ". Skipping calibration."; + return false; + } + + if (!timingCalibartionData_.leadingTimeVsLs.contains(planeKey)) { + std::string planeName; + detId.planeName(planeName); + timingCalibartionData_.leadingTimeVsLs[planeKey] = iGetter.get(dqmDir_ + "/tvsls_" + planeName); + if (!timingCalibartionData_.leadingTimeVsLs.at(planeKey)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve leading time vs. time over threshold monitor for plane (" << planeName + << "). Skipping calibration."; + return false; } - double upper_tot_range = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos); - double lower_tot_range = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos); - - { // scope for x-profile - - std::string ch_name; - detid.channelName(ch_name); - auto profile = iBooker.bookProfile(ch_name + "_prof_x", ch_name + "_prof_x", 240, 0., 60., 450, -20., 25.); - - std::unique_ptr prof(hists.leadingTimeVsToT[chid]->getTH2F()->ProfileX("_prof_x", 1, -1)); - *(profile->getTProfile()) = *((TProfile*)prof->Clone()); - profile->getTProfile()->SetTitle(ch_name.c_str()); - profile->getTProfile()->SetName(ch_name.c_str()); - - interp_.SetParameters(hists.leadingTime[chid]->getRMS(), - hists.toT[chid]->getMean(), - 0.8, - hists.leadingTime[chid]->getMean() - hists.leadingTime[chid]->getRMS()); - const auto& res = profile->getTProfile()->Fit(&interp_, "B+", "", lower_tot_range, 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(offset_, resolution_); // hardcoded offset/resolution 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 << ")."; + } + + return true; +} + +//------------------------------------------------------------------------------ + +TProfile* PPSTimingCalibrationPCLHarvester::createTVsTotProfile(DQMStore::IBooker& iBooker, + dqm::reco::MonitorElement* tVsTot, + const std::string& channelName) const { + const MonitorElement* const monitorProfile{ + iBooker.bookProfile(channelName, channelName, 240, 0.0, 60.0, 450, -20.0, 25.0)}; + TProfile* const monitorTProfile{monitorProfile->getTProfile()}; + const std::unique_ptr tVsTotProfile{tVsTot->getTH2F()->ProfileX()}; + *(monitorTProfile) = *(static_cast(tVsTotProfile->Clone())); + + const char* const profileName = channelName.c_str(); + monitorTProfile->SetTitle(profileName); + monitorTProfile->SetName(profileName); + monitorTProfile->SetYTitle("Average t (ns)"); + + return monitorTProfile; +} + +std::pair PPSTimingCalibrationPCLHarvester::findFitRange( + const TH1F* tot, + const int maxTotBin, + const double lowerThresholdFractionOfMax, + const double upperThresholdFractionOfMax) const { + constexpr double totLowerLimitRangeSearch{8.0}; + const TAxis* const totXAxis{tot->GetXaxis()}; + const double maxTotCount{tot->GetBinContent(maxTotBin)}; + + double lowerTotRange{8}; + if (lowerThresholdFractionOfMax != FixedFitBoundIndication_) { + int lowerLimitPos{maxTotBin}; + const double lowerThreshold{lowerThresholdFractionOfMax * maxTotCount}; + while (tot->GetBinContent(lowerLimitPos) >= lowerThreshold && + totXAxis->GetBinCenter(lowerLimitPos) > totLowerLimitRangeSearch) { + --lowerLimitPos; } + lowerTotRange = totXAxis->GetBinCenter(lowerLimitPos); } - // fill the DB object record - PPSTimingCalibration calib(formula_, calib_params, calib_time); + constexpr double totUpperLimitRangeSearch{20.0}; + double upperTotRange{15}; + if (upperThresholdFractionOfMax != FixedFitBoundIndication_) { + int upperLimitPos{maxTotBin}; + const double upperThreshold{upperThresholdFractionOfMax * maxTotCount}; + while (tot->GetBinContent(upperLimitPos) >= upperThreshold && + totXAxis->GetBinCenter(upperLimitPos) < totUpperLimitRangeSearch) { + ++upperLimitPos; + } + upperTotRange = totXAxis->GetBinCenter(upperLimitPos); + } - // write the object - edm::Service poolDbService; - poolDbService->writeOneIOV(calib, poolDbService->currentTime(), "PPSTimingCalibrationRcd_HPTDC"); + return {lowerTotRange, upperTotRange}; } //------------------------------------------------------------------------------ @@ -201,9 +303,9 @@ void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescri ->setComment("input path for the various DQM plots"); desc.add("formula", "[0]/(exp((x-[1])/[2])+1)+[3]") ->setComment("interpolation formula for the time walk component"); + desc.add("tVsLsFilename", "") + ->setComment("ROOT filename with t vs LS histogram for double peak correction"); desc.add("minEntries", 100)->setComment("minimal number of hits to extract calibration"); - desc.add("thresholdFractionOfMax", 0.05) - ->setComment("threshold for TOT fit, defined as percent of max count in 1D TOT"); descriptions.addWithDefaultLabel(desc); } diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc index e0196ac412aac..0e51e6a162b48 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc @@ -5,63 +5,86 @@ * Edoardo Bossini * Piotr Maciej Cwiklicki * Laurent Forthomme + * Tomasz Ostafin * ****************************************************************************/ +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" + +#include "DataFormats/Common/interface/DetSetVector.h" +#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" +#include "DataFormats/CTPPSReco/interface/CTPPSDiamondRecHit.h" + #include "DQMServices/Core/interface/DQMGlobalEDAnalyzer.h" #include "DQMServices/Core/interface/DQMStore.h" +#include "FWCore/Framework/interface/Event.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" +#include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h" //------------------------------------------------------------------------------ -class PPSTimingCalibrationPCLWorker : public DQMGlobalEDAnalyzer { +class PPSTimingCalibrationPCLWorker : public DQMGlobalEDAnalyzer { public: explicit PPSTimingCalibrationPCLWorker(const edm::ParameterSet&); - void dqmAnalyze(const edm::Event&, const edm::EventSetup&, const TimingCalibrationHistograms&) const override; + void dqmBeginRun(const edm::Run&, const edm::EventSetup&, TimingCalibrationData&) const override; + void dqmAnalyze(const edm::Event&, const edm::EventSetup&, const TimingCalibrationData&) const override; static void fillDescriptions(edm::ConfigurationDescriptions&); private: + using DiamondRecHitVector = edm::DetSetVector; + void bookHistograms(DQMStore::IBooker&, const edm::Run&, const edm::EventSetup&, - TimingCalibrationHistograms&) const override; + TimingCalibrationData&) const override; - template - bool searchForProduct(edm::Event const& iEvent, - const std::vector>& tokens, + void searchForProduct(const edm::Event& iEvent, + const std::vector>& tokens, const std::vector& tags, - edm::Handle& handle) const; - - const std::vector RecHitTags_; - std::vector>> diamondRecHitTokens_; - const edm::ESGetToken geomEsToken_; + edm::Handle& handle) const; const std::string dqmDir_; + const std::string tVsLsFilename_; + std::vector> diamondRecHitTokens_; + const std::vector recHitTags_; + const edm::ESGetToken geomEsToken_; }; //------------------------------------------------------------------------------ PPSTimingCalibrationPCLWorker::PPSTimingCalibrationPCLWorker(const edm::ParameterSet& iConfig) - : RecHitTags_(iConfig.getParameter>("diamondRecHitTags")), - geomEsToken_(esConsumes()), - dqmDir_(iConfig.getParameter("dqmDir")) { - for (auto& tag : RecHitTags_) - diamondRecHitTokens_.push_back(consumes>(tag)); + : dqmDir_{iConfig.getParameter("dqmDir")}, + tVsLsFilename_{iConfig.getParameter("tVsLsFilename")}, + recHitTags_{iConfig.getParameter>("diamondRecHitTags")}, + geomEsToken_{esConsumes()} { + for (const auto& tag : recHitTags_) { + diamondRecHitTokens_.push_back(consumes(tag)); + } +} + +//------------------------------------------------------------------------------ + +void PPSTimingCalibrationPCLWorker::dqmBeginRun(const edm::Run& iRun, + const edm::EventSetup& iSetup, + TimingCalibrationData& timingCalibrationData) const { + std::vector detIds; + const auto& geom = iSetup.getData(geomEsToken_); + for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) { + if (CTPPSDiamondDetId::check(it->first)) { + const CTPPSDiamondDetId detId{it->first}; + detIds.push_back(detId); + } + } + timingCalibrationData.doublePeakCorrection.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds); } //------------------------------------------------------------------------------ @@ -69,22 +92,32 @@ PPSTimingCalibrationPCLWorker::PPSTimingCalibrationPCLWorker(const edm::Paramete void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, const edm::Run& iRun, const edm::EventSetup& iSetup, - TimingCalibrationHistograms& iHists) const { + TimingCalibrationData& timingCalibrationData) const { iBooker.cd(); iBooker.setCurrentFolder(dqmDir_); - std::string ch_name; + std::string planeName; + std::string channelName; 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); - - 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.); + if (CTPPSDiamondDetId::check(it->first)) { + const CTPPSDiamondDetId detId{it->first}; + const uint32_t channelId{detId.rawId()}; + detId.channelName(channelName); + const PlaneKey planeKey{detId.arm(), detId.station(), detId.plane()}; + + timingCalibrationData.leadingTime[channelId] = + iBooker.book1D("t_" + channelName, channelName + ";t (ns);Entries", 1200, -60.0, 60.0); + timingCalibrationData.toT[channelId] = + iBooker.book1D("tot_" + channelName, channelName + ";ToT (ns);Entries", 160, -20.0, 20.0); + timingCalibrationData.leadingTimeVsToT[channelId] = + iBooker.book2D("tvstot_" + channelName, channelName + ";ToT (ns);t (ns)", 240, 0.0, 60.0, 450, -20.0, 25.0); + if (!timingCalibrationData.leadingTimeVsLs.contains(planeKey)) { + detId.planeName(planeName); + timingCalibrationData.leadingTimeVsLs[planeKey] = + iBooker.book2D("tvsls_" + planeName, planeName + ";LS;t (ns)", 3000, 1.0, 3001.0, 500, 0.0, 25.0); + } + } } } @@ -92,64 +125,77 @@ void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, void PPSTimingCalibrationPCLWorker::dqmAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, - const TimingCalibrationHistograms& iHists) const { - edm::Handle> dsv_rechits; + const TimingCalibrationData& timingCalibrationData) const { + edm::Handle dsvRechits; // then extract the rechits information for later processing - searchForProduct(iEvent, diamondRecHitTokens_, RecHitTags_, dsv_rechits); + searchForProduct(iEvent, diamondRecHitTokens_, recHitTags_, dsvRechits); // 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."; + if (dsvRechits->empty()) { + edm::LogWarning{"PPSTimingCalibrationPCLWorker:dqmAnalyze"} + << "No rechits retrieved from the event content for event " << iEvent.id() << '.'; 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."; + + const unsigned int ls{iEvent.luminosityBlock()}; + std::string planeName; + for (const auto& dsRechits : *dsvRechits) { + const CTPPSDiamondDetId detId{dsRechits.detId()}; + const uint32_t channelId{detId.rawId()}; + if (!timingCalibrationData.leadingTimeVsToT.contains(channelId)) { + edm::LogWarning{"PPSTimingCalibrationPCLWorker:dqmAnalyze"} << "Pad with Detector ID =" << detId + << " is not set to be monitored."; continue; } - for (const auto& rechit : ds_rechits) { + + detId.planeName(planeName); + const PlaneKey planeKey{detId.arm(), detId.station(), detId.plane()}; + for (const auto& rechit : dsRechits) { // 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()); + if (rechit.time() != 0.0 && rechit.toT() >= 0.0) { + const double correctedLeadingTime{ + timingCalibrationData.doublePeakCorrection.getCorrectedLeadingTime(rechit.time(), ls, planeKey)}; + timingCalibrationData.leadingTime.at(channelId)->Fill(correctedLeadingTime); + timingCalibrationData.toT.at(channelId)->Fill(rechit.toT()); + timingCalibrationData.leadingTimeVsToT.at(channelId)->Fill(rechit.toT(), correctedLeadingTime); + timingCalibrationData.leadingTimeVsLs.at(planeKey)->Fill(ls, correctedLeadingTime); + } } } } //------------------------------------------------------------------------------ -void PPSTimingCalibrationPCLWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { - edm::ParameterSetDescription desc; - desc.add>("diamondRecHitTags", {edm::InputTag("ctppsDiamondRecHits")}) - ->setComment("input tag for the PPS diamond detectors rechits"); - desc.add("dqmDir", "AlCaReco/PPSTimingCalibrationPCL") - ->setComment("output path for the various DQM plots"); - - descriptions.addWithDefaultLabel(desc); -} - -template -bool PPSTimingCalibrationPCLWorker::searchForProduct(edm::Event const& iEvent, - const std::vector>& tokens, +void PPSTimingCalibrationPCLWorker::searchForProduct(const edm::Event& iEvent, + const std::vector>& tokens, const std::vector& tags, - edm::Handle& handle) const { - bool foundProduct = false; - for (unsigned int i = 0; i < tokens.size(); i++) - if (auto h = iEvent.getHandle(tokens[i])) { + edm::Handle& handle) const { + bool foundProduct{false}; + for (unsigned int i{0}; i < tokens.size(); ++i) + if (const auto h = iEvent.getHandle(tokens[i])) { handle = h; foundProduct = true; - edm::LogInfo("searchForProduct") << "Found a product with " << tags[i]; + edm::LogInfo{"searchForProduct"} << "Found a product with " << tags[i] << '.'; break; } - if (!foundProduct) - throw edm::Exception(edm::errors::ProductNotFound) << "Could not find a product with any of the selected labels."; + if (!foundProduct) { + throw edm::Exception{edm::errors::ProductNotFound} << "Could not find a product with any of the selected labels."; + } +} - return foundProduct; +//------------------------------------------------------------------------------ + +void PPSTimingCalibrationPCLWorker::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("dqmDir", "AlCaReco/PPSTimingCalibrationPCL") + ->setComment("output path for the various DQM plots"); + desc.add("tVsLsFilename", "") + ->setComment("ROOT filename with t vs LS histogram for double peak correction"); + desc.add>("diamondRecHitTags", {edm::InputTag{"ctppsDiamondRecHits"}}) + ->setComment("input tag for the PPS diamond detectors rechits"); + + descriptions.addWithDefaultLabel(desc); } DEFINE_FWK_MODULE(PPSTimingCalibrationPCLWorker); diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc new file mode 100644 index 0000000000000..c4da6018d0cb8 --- /dev/null +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -0,0 +1,151 @@ +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" + +#include "FWCore/Utilities/interface/EDMException.h" + +#include "TF1.h" +#include "TFitResult.h" + +void DoublePeakCorrection::extractLsAndTimeOffset(const std::string& tVsLsFilename, + const unsigned int run, + const std::vector& detIds) { + if (!tVsLsFilename.empty()) { + TFile tVsLsFile{tVsLsFilename.c_str()}; + if (tVsLsFile.IsOpen()) { + const std::string runNumber{std::to_string(run)}; + for (const auto& detId : detIds) { + const PlaneKey planeKey{detId.arm(), detId.station(), detId.plane()}; + fillLsAndTimeOffset(getTVsLs(tVsLsFile, runNumber, detId), planeKey); + } + } else { + throw edm::Exception{edm::errors::FileOpenError} << "Can't open the file with t vs LS: " << tVsLsFilename << '.'; + } + } +} + +const TH2F* DoublePeakCorrection::getTVsLs(TFile& tVsLsFile, + const std::string& runNumber, + const CTPPSDiamondDetId& detId) { + std::string planeName; + detId.planeName(planeName); + std::string tVsLsHistPath{"DQMData/Run " + runNumber + "/AlCaReco/Run summary/PPSTimingCalibrationPCL/tvsls_" + + planeName}; + const auto* tVsLs = tVsLsFile.Get(tVsLsHistPath.c_str()); + if (tVsLs) { + return tVsLs; + } + throw edm::Exception{edm::errors::FileReadError} << "Can't open the t vs LS histogram: " << tVsLsHistPath << '.'; +} + +void DoublePeakCorrection::fillLsAndTimeOffset(const TH2F* tVsLs, const PlaneKey& planeKey) { + if (!lsAndTimeOffsets_.contains(planeKey)) { + const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); + if (doublePeakLs != 1) { + lsAndTimeOffsets_[planeKey] = { + doublePeakLs, findTimeOffset(tVsLs, doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; + } + } +} + +std::tuple DoublePeakCorrection::findLsAndTimePeaks(const TH2F* tVsLs, + const PlaneKey& planeKey) const { + auto numOfLsBins = static_cast(tVsLs->GetNbinsX()); + auto numOfTBins = static_cast(tVsLs->GetNbinsY()); + double firstPeakTWithMaxCount{0.0}; + double secondPeakTWithMaxCount{0.0}; + unsigned int tDiffCount{0}; + unsigned int doublePeakLs{1}; + for (unsigned int lsBin{1}; lsBin <= numOfLsBins; ++lsBin) { + double tMaxCount{0}; + for (unsigned int tBin{1}; tBin <= numOfTBins; ++tBin) { + const double tCount{tVsLs->GetBinContent(lsBin, tBin)}; + const double tBinCenter{tVsLs->GetYaxis()->GetBinCenter(tBin)}; + constexpr double minTCount{10.0}; + if (tCount >= minTCount && tCount > tMaxCount) { + tMaxCount = tCount; + secondPeakTWithMaxCount = tBinCenter; + } + } + + constexpr double tMinDiff{1.0}; + if (firstPeakTWithMaxCount != 0.0 && std::abs(secondPeakTWithMaxCount - firstPeakTWithMaxCount) > tMinDiff) { + constexpr unsigned int minTDiffCount{5}; + ++tDiffCount; + if (tDiffCount == 1) { + doublePeakLs = lsBin; + } else if (tDiffCount == minTDiffCount) { + return {doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; + } + } else { + tDiffCount = 0; + firstPeakTWithMaxCount = secondPeakTWithMaxCount; + } + } + + return {1, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; +} + +double DoublePeakCorrection::findTimeOffset(const TH2F* tVsLs, + const unsigned int doublePeakLs, + const double firstPeakEstimatedMean, + const double secondPeakEstimatedMean) const { + const std::unique_ptr firstPeak{tVsLs->ProjectionY("_py", 1, doublePeakLs - 1)}; + const std::unique_ptr secondPeak{tVsLs->ProjectionY("_py", doublePeakLs, -1)}; + return findGaussianMean(secondPeak, secondPeakEstimatedMean) - findGaussianMean(firstPeak, firstPeakEstimatedMean); +} + +double DoublePeakCorrection::findGaussianMean(const std::unique_ptr& peak, const double estimatedMean) const { + constexpr unsigned int meanParamIndex{1}; + constexpr double fitSigma{2.5}; + const double fitLeftBound{estimatedMean - fitSigma}; + const double fitRightBound{estimatedMean + fitSigma}; + TF1 fitFunction{"peak", "gaus"}; + fitFunction.SetParLimits(meanParamIndex, fitLeftBound, fitRightBound); + fitFunction.SetParameter(meanParamIndex, estimatedMean); + const TFitResultPtr& peakFit{peak->Fit(&fitFunction, "NS", "", fitLeftBound, fitRightBound)}; + if (peakFit->IsValid()) { + return peakFit->Parameter(meanParamIndex); + } + throw edm::Exception{edm::errors::FatalRootError} << "Double peak Gaussian fit not valid."; +} + +bool DoublePeakCorrection::isCorrectionNeeded(const PlaneKey& planeKey) const { + return lsAndTimeOffsets_.contains(planeKey); +} + +double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, + const unsigned int ls, + const PlaneKey& planeKey) const { + if (auto it = lsAndTimeOffsets_.find(planeKey); it != std::end(lsAndTimeOffsets_)) { + const auto [doublePeakLs, doublePeakTimeOffset] = it->second; + if (ls >= doublePeakLs) { + return leadingTime - doublePeakTimeOffset; + } + } + return leadingTime; +} + +double DoublePeakCorrection::getEncodedLsAndTimeOffset(const PlaneKey& planeKey) const { + if (auto it = lsAndTimeOffsets_.find(planeKey); it != std::end(lsAndTimeOffsets_)) { + const auto [doublePeakLs, doublePeakTimeOffset] = it->second; + if (doublePeakTimeOffset > 0.0) { + return doublePeakLs * LsEncodingMultiple_ + doublePeakTimeOffset; + } + return -(doublePeakLs * LsEncodingMultiple_ - doublePeakTimeOffset); + } + return 0.0; +} + +double DoublePeakCorrection::GetCorrectedLeadingTime(const double leadingTime, + const unsigned int ls, + const double encodedLsAndTimeOffset) { + if (leadingTime != 0.0 && encodedLsAndTimeOffset != 0.0) { + const auto doublePeakLs = static_cast(std::abs(encodedLsAndTimeOffset) / LsEncodingMultiple_); + if (ls >= doublePeakLs) { + const double doublePeakTimeOffset = encodedLsAndTimeOffset >= 0 + ? encodedLsAndTimeOffset - doublePeakLs * LsEncodingMultiple_ + : encodedLsAndTimeOffset + doublePeakLs * LsEncodingMultiple_; + return leadingTime - doublePeakTimeOffset; + } + } + return leadingTime; +} diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py index abf1088fe5fce..bd8da12600723 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py @@ -1,54 +1,81 @@ -run = 357902 -input_file=['/store/express/Run2022D/StreamALCAPPSExpress/ALCAPROMPT/PromptCalibProdPPSTimingCalib-Express-v2/000/357/902/00000/123861b4-b632-4835-91b6-1e586d34509e.root'] import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing + +from Configuration.AlCa.GlobalTag import GlobalTag from Configuration.StandardSequences.Eras import eras + + process = cms.Process("harvester", eras.Run3) +options = VarParsing.VarParsing() + +options.register("globalTag", + "", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Global Tag" +) +options.register("inputFiles", + "", + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string +) +options.register("tVsLsFilename", + "", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "ROOT filename with t vs LS histogram for double peak correction" +) + +options.parseArguments() + process.load("FWCore.MessageService.MessageLogger_cfi") +process.MessageLogger.cerr.FwkReport.reportEvery = 1_000_000 -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(-1)) -process.load('Configuration.StandardSequences.Services_cff') -process.load('Configuration.EventContent.EventContent_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') -from Configuration.AlCa.GlobalTag import GlobalTag -from Configuration.AlCa.autoCond import autoCond -process.GlobalTag = GlobalTag(process.GlobalTag, autoCond['run3_data_prompt'], '') +process.load("Configuration.StandardSequences.Services_cff") +process.load("Configuration.EventContent.EventContent_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.GlobalTag = GlobalTag(process.GlobalTag, options.globalTag, "") # Source (histograms) +fileList = [ + f"file:{f}" + if not (f.startswith("/store/") or f.startswith("file:") or f.startswith("root:")) + else f for f in options.inputFiles +] + process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - input_file - ), + fileNames = cms.untracked.vstring(fileList) ) # output service for database -process.load('CondCore.CondDB.CondDB_cfi') -process.CondDB.connect = 'sqlite_file:ppsDiamondTiming_calibration'+str(run)+'.sqlite' # SQLite output +process.load("CondCore.CondDB.CondDB_cfi") +process.CondDB.connect = "sqlite_file:ppsDiamondTiming_calibration.sqlite" # SQLite output -process.PoolDBOutputService = cms.Service('PoolDBOutputService', +process.PoolDBOutputService = cms.Service("PoolDBOutputService", process.CondDB, - timetype = cms.untracked.string('runnumber'), + timetype = cms.untracked.string("runnumber"), toPut = cms.VPSet( cms.PSet( - record = cms.string('PPSTimingCalibrationRcd_HPTDC'), - tag = cms.string('DiamondTimingCalibration'), + record = cms.string("PPSTimingCalibrationRcd_HPTDC"), + tag = cms.string("DiamondTimingCalibration") ) ) ) process.load("CalibPPS.TimingCalibration.ppsTimingCalibrationPCLHarvester_cfi") -#process.PPSDiamondSampicTimingCalibrationPCLHarvester.jsonCalibFile="initial.cal.json" +process.ppsTimingCalibrationPCLHarvester.tVsLsFilename = options.tVsLsFilename # load DQM framework process.load("DQMServices.Core.DQMStore_cfi") process.load("DQMServices.Components.DQMEnvironment_cfi") process.dqmEnv.subSystemFolder = "CalibPPS" -process.dqmSaver.convention = 'Offline' -process.dqmSaver.workflow = "/CalibPPS/TimingCalibration/CMSSW_12_6_0_pre2" +process.dqmSaver.convention = "Offline" +process.dqmSaver.workflow = "/CalibPPS/TimingCalibration/CMSSW_14_1_4" process.dqmSaver.saveByRun = -1 process.dqmSaver.saveAtJobEnd = True -process.dqmSaver.forceRunNumber = run process.DQMStore = cms.Service("DQMStore") @@ -61,11 +88,8 @@ process.EDMtoMEConverter.lumiInputTag = "MEtoEDMConvertPPSTimingCalib:MEtoEDMConverterLumi" process.EDMtoMEConverter.runInputTag = "MEtoEDMConvertPPSTimingCalib:MEtoEDMConverterRun" -#import FWCore.PythonUtilities.LumiList as LumiList -#process.source.lumisToProcess = LumiList.LumiList(filename = 'allrunsSB-PPS-forCalib.json').getVLuminosityBlockRange() - process.p = cms.Path( - process.EDMtoMEConverter* + process.EDMtoMEConverter * process.ppsTimingCalibrationPCLHarvester ) @@ -78,6 +102,3 @@ process.p, process.end_path ) - - - diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py index 140576769fe18..a454f20396144 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py @@ -1,65 +1,82 @@ import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing +import FWCore.PythonUtilities.LumiList as LumiList + +from Configuration.AlCa.GlobalTag import GlobalTag + process = cms.Process("worker") +options = VarParsing.VarParsing() +options.register("globalTag", + "", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Global Tag" +) +options.register("inputFiles", + "", + VarParsing.VarParsing.multiplicity.list, + VarParsing.VarParsing.varType.string +) +options.register("tVsLsFilename", + "", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "ROOT filename with t vs LS histogram for double peak correction" +) +options.register("jsonFileName", + "", + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "Certification JSON filename" +) + +options.parseArguments() + process.load("FWCore.MessageService.MessageLogger_cfi") -process.load('RecoPPS.Local.totemTimingLocalReconstruction_cff') +process.MessageLogger.cerr.FwkReport.reportEvery = 1_000_000 + +process.load("RecoPPS.Local.totemTimingLocalReconstruction_cff") process.source = cms.Source("EmptySource") -process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(10000) ) +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(-1)) process.a1 = cms.EDAnalyzer("StreamThingAnalyzer", - product_to_get = cms.string('m1') + product_to_get = cms.string("m1") ) -process.load('Configuration.StandardSequences.Services_cff') -process.load('Configuration.EventContent.EventContent_cff') -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') -from Configuration.AlCa.GlobalTag import GlobalTag -from Configuration.AlCa.autoCond import autoCond -process.GlobalTag = GlobalTag(process.GlobalTag, autoCond['run3_data_prompt'], '') -process.load("EventFilter.CTPPSRawToDigi.ctppsRawToDigi_cff") +process.load("Configuration.StandardSequences.Services_cff") +process.load("Configuration.EventContent.EventContent_cff") +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +process.GlobalTag = GlobalTag(process.GlobalTag, options.globalTag, "") process.load("RecoPPS.Configuration.recoCTPPS_cff") -#process.load('CondCore.CondDB.CondDB_cfi') -#process.CondDB.connect = 'sqlite_file:ppsDiamondTiming_calibration.sqlite' # SQLite input -#process.PoolDBESSource = cms.ESSource('PoolDBESSource', -# process.CondDB, -# DumpStats = cms.untracked.bool(True), -# toGet = cms.VPSet( -# cms.PSet( -# record = cms.string('PPSTimingCalibrationRcd'), -# tag = cms.string('DiamondTimingCalibration') -# ) -# ) -#) - process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - '/store/data/Run2022B/AlCaPPS/RAW/v1/000/355/207/00000/c23440f4-49c0-44aa-b8f6-f40598fb4705.root', - - -), + fileNames = cms.untracked.vstring(options.inputFiles) ) +if options.jsonFileName != "": + process.source.lumisToProcess = LumiList.LumiList(filename=options.jsonFileName).getVLuminosityBlockRange() + print(f"Using JSON file: {options.jsonFileName}") + process.load("CalibPPS.TimingCalibration.ppsTimingCalibrationPCLWorker_cfi") +process.ppsTimingCalibrationPCLWorker.tVsLsFilename = options.tVsLsFilename + process.DQMStore = cms.Service("DQMStore") -process.dqmOutput = cms.OutputModule("DQMRootOutputModule", - fileName = cms.untracked.string("worker_output.root") +process.dqmOutput = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string("file:worker_output.root"), + outputCommands = cms.untracked.vstring( + "drop *", + "keep *_MEtoEDMConvertPPSTimingCalib_*_*" + ) ) process.load("CalibPPS.TimingCalibration.ALCARECOPromptCalibProdPPSTimingCalib_cff") -process.ctppsPixelDigis.inputLabel = "hltPPSCalibrationRaw" -process.ctppsDiamondRawToDigi.rawDataTag = "hltPPSCalibrationRaw" -process.totemRPRawToDigi.rawDataTag = "hltPPSCalibrationRaw" -process.totemTimingRawToDigi.rawDataTag = "hltPPSCalibrationRaw" - process.path = cms.Path( - #process.a1* - process.ctppsRawToDigi * - process.recoCTPPS * - process.ppsTimingCalibrationPCLWorker + process.ppsTimingCalibrationPCLWorker * + process.MEtoEDMConvertPPSTimingCalib ) process.end_path = cms.EndPath( diff --git a/RecoPPS/Local/BuildFile.xml b/RecoPPS/Local/BuildFile.xml index 8ed2b4e1bb8eb..545c11cc96c10 100644 --- a/RecoPPS/Local/BuildFile.xml +++ b/RecoPPS/Local/BuildFile.xml @@ -11,6 +11,7 @@ + diff --git a/RecoPPS/Local/interface/CTPPSDiamondRecHitProducerAlgorithm.h b/RecoPPS/Local/interface/CTPPSDiamondRecHitProducerAlgorithm.h index 7505d9659f908..f77e7ceffeee6 100644 --- a/RecoPPS/Local/interface/CTPPSDiamondRecHitProducerAlgorithm.h +++ b/RecoPPS/Local/interface/CTPPSDiamondRecHitProducerAlgorithm.h @@ -25,7 +25,9 @@ class CTPPSDiamondRecHitProducerAlgorithm using TimingRecHitProducerAlgorithm::TimingRecHitProducerAlgorithm; void build(const CTPPSGeometry&, const edm::DetSetVector&, - edm::DetSetVector&) override; + edm::DetSetVector&, + const unsigned int, + const bool) override; private: static constexpr unsigned short MAX_CHANNEL = 20; diff --git a/RecoPPS/Local/interface/TimingRecHitProducerAlgorithm.h b/RecoPPS/Local/interface/TimingRecHitProducerAlgorithm.h index ec96d45c02496..a46059560d1fb 100644 --- a/RecoPPS/Local/interface/TimingRecHitProducerAlgorithm.h +++ b/RecoPPS/Local/interface/TimingRecHitProducerAlgorithm.h @@ -31,7 +31,7 @@ class TimingRecHitProducerAlgorithm { calibLUT_ = &calibLUT; calib_fct_ = std::make_unique(calib_->formula()); } - virtual void build(const G&, const D&, R&) = 0; + virtual void build(const G&, const D&, R&, const unsigned int, const bool) = 0; protected: /// Conversion constant between time slice and absolute time (in ns) diff --git a/RecoPPS/Local/interface/TotemT2RecHitProducerAlgorithm.h b/RecoPPS/Local/interface/TotemT2RecHitProducerAlgorithm.h index bf21da364e5a8..92fa317006338 100644 --- a/RecoPPS/Local/interface/TotemT2RecHitProducerAlgorithm.h +++ b/RecoPPS/Local/interface/TotemT2RecHitProducerAlgorithm.h @@ -24,7 +24,9 @@ class TotemT2RecHitProducerAlgorithm : public TimingRecHitProducerAlgorithm&, - edmNew::DetSetVector&) override; + edmNew::DetSetVector&, + const unsigned int, + const bool) override; }; #endif diff --git a/RecoPPS/Local/interface/TotemTimingRecHitProducerAlgorithm.h b/RecoPPS/Local/interface/TotemTimingRecHitProducerAlgorithm.h index cbb4214044ffc..b5ff37b39cecd 100644 --- a/RecoPPS/Local/interface/TotemTimingRecHitProducerAlgorithm.h +++ b/RecoPPS/Local/interface/TotemTimingRecHitProducerAlgorithm.h @@ -31,7 +31,9 @@ class TotemTimingRecHitProducerAlgorithm : public TimingRecHitProducerAlgorithm< void setCalibration(const PPSTimingCalibration&); void build(const CTPPSGeometry&, const edm::DetSetVector&, - edm::DetSetVector&) override; + edm::DetSetVector&, + const unsigned int, + const bool) override; private: struct RegressionResults { diff --git a/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc b/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc index 9f62a6383c6a6..d2a27c275ee3b 100644 --- a/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc +++ b/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc @@ -40,6 +40,7 @@ class CTPPSDiamondRecHitProducer : public edm::stream::EDProducer<> { static void fillDescriptions(edm::ConfigurationDescriptions&); private: + void beginRun(edm::Run const&, edm::EventSetup const&) override; void produce(edm::Event&, const edm::EventSetup&) override; edm::EDGetTokenT > digiToken_; @@ -51,15 +52,17 @@ class CTPPSDiamondRecHitProducer : public edm::stream::EDProducer<> { edm::ESWatcher calibWatcher_; edm::ESWatcher lutWatcher_; - bool applyCalib_; CTPPSDiamondRecHitProducerAlgorithm algo_; + bool applyCalib_; + bool isRun2_; }; CTPPSDiamondRecHitProducer::CTPPSDiamondRecHitProducer(const edm::ParameterSet& iConfig) - : digiToken_(consumes >(iConfig.getParameter("digiTag"))), - geometryToken_(esConsumes()), - applyCalib_(iConfig.getParameter("applyCalibration")), - algo_(iConfig) { + : digiToken_{consumes >(iConfig.getParameter("digiTag"))}, + geometryToken_{esConsumes()}, + algo_{iConfig}, + applyCalib_{iConfig.getParameter("applyCalibration")}, + isRun2_{false} { if (applyCalib_) { timingCalibrationToken_ = esConsumes( edm::ESInputTag(iConfig.getParameter("timingCalibrationTag"))); @@ -68,6 +71,11 @@ CTPPSDiamondRecHitProducer::CTPPSDiamondRecHitProducer(const edm::ParameterSet& produces >(); } +void CTPPSDiamondRecHitProducer::beginRun(const edm::Run& iRun, const edm::EventSetup&) { + const unsigned int run{iRun.run()}; + isRun2_ = run <= 299'999 && run >= 200'000; +} + void CTPPSDiamondRecHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { auto pOut = std::make_unique >(); @@ -75,11 +83,11 @@ void CTPPSDiamondRecHitProducer::produce(edm::Event& iEvent, const edm::EventSet const auto& digis = iEvent.get(digiToken_); if (!digis.empty()) { - if (applyCalib_ && (calibWatcher_.check(iSetup) or lutWatcher_.check(iSetup))) + if (applyCalib_ && (calibWatcher_.check(iSetup) || lutWatcher_.check(iSetup))) algo_.setCalibration(iSetup.getData(timingCalibrationToken_), iSetup.getData(timingCalibrationLUTToken_)); // produce the rechits collection - algo_.build(iSetup.getData(geometryToken_), digis, *pOut); + algo_.build(iSetup.getData(geometryToken_), digis, *pOut, iEvent.luminosityBlock(), isRun2_); } iEvent.put(std::move(pOut)); diff --git a/RecoPPS/Local/plugins/TotemT2RecHitProducer.cc b/RecoPPS/Local/plugins/TotemT2RecHitProducer.cc index 2167c0c246df0..f52798b8e95c1 100644 --- a/RecoPPS/Local/plugins/TotemT2RecHitProducer.cc +++ b/RecoPPS/Local/plugins/TotemT2RecHitProducer.cc @@ -61,7 +61,7 @@ void TotemT2RecHitProducer::produce(edm::Event& iEvent, const edm::EventSetup& i if (!digis.empty()) { // produce the rechits collection - algo_.build(iSetup.getData(geometryToken_), digis, *pOut); + algo_.build(iSetup.getData(geometryToken_), digis, *pOut, 0.0, false); } iEvent.put(std::move(pOut)); diff --git a/RecoPPS/Local/plugins/TotemTimingRecHitProducer.cc b/RecoPPS/Local/plugins/TotemTimingRecHitProducer.cc index aa8327fb0a1ea..40db77be3b981 100644 --- a/RecoPPS/Local/plugins/TotemTimingRecHitProducer.cc +++ b/RecoPPS/Local/plugins/TotemTimingRecHitProducer.cc @@ -78,7 +78,7 @@ void TotemTimingRecHitProducer::produce(edm::Event& iEvent, const edm::EventSetu algo_.setCalibration(iSetup.getData(timingCalibrationToken_)); // produce the rechits collection - algo_.build(iSetup.getData(geometryToken_), digis, *pOut); + algo_.build(iSetup.getData(geometryToken_), digis, *pOut, 0.0, false); } iEvent.put(std::move(pOut)); diff --git a/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc b/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc index 6421f6f3a91c2..f4627acb110b4 100644 --- a/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc +++ b/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc @@ -6,16 +6,21 @@ * ****************************************************************************/ +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" + +#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" + #include "FWCore/Utilities/interface/isFinite.h" #include "RecoPPS/Local/interface/CTPPSDiamondRecHitProducerAlgorithm.h" -#include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" //---------------------------------------------------------------------------------------------------- void CTPPSDiamondRecHitProducerAlgorithm::build(const CTPPSGeometry& geom, const edm::DetSetVector& input, - edm::DetSetVector& output) { + edm::DetSetVector& output, + const unsigned int ls, + const bool is_run_2) { for (const auto& vec : input) { const CTPPSDiamondDetId detid(vec.detId()); @@ -69,11 +74,15 @@ void CTPPSDiamondRecHitProducerAlgorithm::build(const CTPPSGeometry& geom, } } - const int time_slice = - (t_lead != 0) ? (t_lead - ch_t_offset / ts_to_ns_) / 1024 : CTPPSDiamondRecHit::TIMESLICE_WITHOUT_LEADING; + // In Run2 time_offset was used for the time window. From Run 3 it's used for correcting the leading edge double peak. + const int time_slice = is_run_2 ? (t_lead != 0) ? (t_lead - ch_t_offset / ts_to_ns_) / 1024 + : CTPPSDiamondRecHit::TIMESLICE_WITHOUT_LEADING + : 0; // calibrated time of arrival const double t0 = (t_lead % 1024) * ts_to_ns_ + lut[t_lead % 1024] * ts_to_ns_ - ch_t_twc; + // DP-corrected time of arrival if needed (no DPs in Run 2) + const double toa = is_run_2 ? t0 : DoublePeakCorrection::GetCorrectedLeadingTime(t0, ls, ch_t_offset); rec_hits.emplace_back( // spatial information x_pos, @@ -83,7 +92,7 @@ void CTPPSDiamondRecHitProducerAlgorithm::build(const CTPPSGeometry& geom, z_pos, z_width, // timing information - t0, + toa, tot, ch_t_precis, time_slice, diff --git a/RecoPPS/Local/src/TotemT2RecHitProducerAlgorithm.cc b/RecoPPS/Local/src/TotemT2RecHitProducerAlgorithm.cc index 18cfe6633dff5..19ea3326b676a 100644 --- a/RecoPPS/Local/src/TotemT2RecHitProducerAlgorithm.cc +++ b/RecoPPS/Local/src/TotemT2RecHitProducerAlgorithm.cc @@ -13,7 +13,9 @@ void TotemT2RecHitProducerAlgorithm::build(const TotemGeometry& geom, const edmNew::DetSetVector& input, - edmNew::DetSetVector& output) { + edmNew::DetSetVector& output, + const unsigned int, + const bool) { for (const auto& vec : input) { const TotemT2DetId detid(vec.detId()); diff --git a/RecoPPS/Local/src/TotemTimingRecHitProducerAlgorithm.cc b/RecoPPS/Local/src/TotemTimingRecHitProducerAlgorithm.cc index 165ebf15f5022..2048f66b6c56d 100644 --- a/RecoPPS/Local/src/TotemTimingRecHitProducerAlgorithm.cc +++ b/RecoPPS/Local/src/TotemTimingRecHitProducerAlgorithm.cc @@ -39,7 +39,9 @@ void TotemTimingRecHitProducerAlgorithm::setCalibration(const PPSTimingCalibrati void TotemTimingRecHitProducerAlgorithm::build(const CTPPSGeometry& geom, const edm::DetSetVector& input, - edm::DetSetVector& output) { + edm::DetSetVector& output, + const unsigned int, + const bool) { for (const auto& vec : input) { const CTPPSDetId detid(vec.detId());