From 33c74548f5ee8fe0e7cabb73337544c5b5fee78d Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 13 Dec 2023 15:17:03 +0100 Subject: [PATCH 01/20] PPS timing calibration automation, worker configuration --- .../test/DiamondCalibrationWorker_cfg.py | 51 ++++++++++++------- 1 file changed, 33 insertions(+), 18 deletions(-) diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py index 140576769fe18..16f41b0c9853a 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py @@ -1,7 +1,22 @@ import FWCore.ParameterSet.Config as cms +import FWCore.ParameterSet.VarParsing as VarParsing 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.parseArguments() + + process.load("FWCore.MessageService.MessageLogger_cfi") process.load('RecoPPS.Local.totemTimingLocalReconstruction_cff') process.source = cms.Source("EmptySource") @@ -16,8 +31,9 @@ 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.GlobalTag = GlobalTag(process.GlobalTag, autoCond['run3_data_prompt'], '') +process.GlobalTag = GlobalTag(process.GlobalTag, options.globalTag,'') +#process.load("EventFilter.CTPPSRawToDigi.ctppsRawToDigi_cff") process.load("RecoPPS.Configuration.recoCTPPS_cff") #process.load('CondCore.CondDB.CondDB_cfi') @@ -34,32 +50,31 @@ #) 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)) process.load("CalibPPS.TimingCalibration.ppsTimingCalibrationPCLWorker_cfi") 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.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.ctppsRawToDigi * + #process.recoCTPPS * + process.ppsTimingCalibrationPCLWorker * + process.MEtoEDMConvertPPSTimingCalib ) process.end_path = cms.EndPath( From ab25dbd70d2113de95edcf15621df9e66918841c Mon Sep 17 00:00:00 2001 From: Christopher Date: Thu, 14 Dec 2023 13:16:53 +0100 Subject: [PATCH 02/20] Add automation harvester configuration --- .../test/DiamondCalibrationHarvester_cfg.py | 28 ++++++++++++------- .../test/DiamondCalibrationWorker_cfg.py | 2 +- 2 files changed, 19 insertions(+), 11 deletions(-) diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py index abf1088fe5fce..fc20758c58e39 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py @@ -1,9 +1,21 @@ -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 from Configuration.StandardSequences.Eras import eras +import FWCore.ParameterSet.VarParsing as VarParsing 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.parseArguments() + process.load("FWCore.MessageService.MessageLogger_cfi") process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) @@ -13,18 +25,15 @@ 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.GlobalTag = GlobalTag(process.GlobalTag, options.globalTag,'') # Source (histograms) process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring( - input_file - ), -) + fileNames = cms.untracked.vstring (options.inputFiles)) # output service for database process.load('CondCore.CondDB.CondDB_cfi') -process.CondDB.connect = 'sqlite_file:ppsDiamondTiming_calibration'+str(run)+'.sqlite' # SQLite output +process.CondDB.connect = 'sqlite_file:ppsDiamondTiming_calibration.sqlite' # SQLite output process.PoolDBOutputService = cms.Service('PoolDBOutputService', process.CondDB, @@ -45,10 +54,9 @@ 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.workflow = "/CalibPPS/TimingCalibration/CMSSW_13_3_0" process.dqmSaver.saveByRun = -1 process.dqmSaver.saveAtJobEnd = True -process.dqmSaver.forceRunNumber = run process.DQMStore = cms.Service("DQMStore") diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py index 16f41b0c9853a..80930401c4e64 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py @@ -20,7 +20,7 @@ process.load("FWCore.MessageService.MessageLogger_cfi") 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') From e47ba05aae32e9ebdf4af6e45197058dceb743c4 Mon Sep 17 00:00:00 2001 From: Christopher Date: Fri, 15 Dec 2023 15:14:33 +0100 Subject: [PATCH 03/20] Harvester argument automatically adds file: for local files --- .../TimingCalibration/test/DiamondCalibrationHarvester_cfg.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py index fc20758c58e39..8ad3f70a7fb4a 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py @@ -28,8 +28,10 @@ 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 (options.inputFiles)) + fileNames = cms.untracked.vstring(fileList)) # output service for database process.load('CondCore.CondDB.CondDB_cfi') From 6960c524a750cc58660e6c85ec77d94132ce28e8 Mon Sep 17 00:00:00 2001 From: Christopher Date: Sat, 20 Jan 2024 16:48:15 +0100 Subject: [PATCH 04/20] Update calibration fit --- .../PPSTimingCalibrationPCLHarvester.cc | 36 +++++++++---------- 1 file changed, 17 insertions(+), 19 deletions(-) diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index c75e2403a7568..800e71c717a55 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -42,10 +42,7 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester { 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; + const double threshold_percent_of_max_; static constexpr double resolution_ = 0.1; static constexpr double offset_ = 0.; TF1 interp_; @@ -58,8 +55,9 @@ PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::Pa 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.) { + threshold_percent_of_max_(iConfig.getParameter("thresholdPercentOfMax")), + interp_("interp", formula_.c_str(), 10.5, 25.){ + // first ensure DB output service is available edm::Service poolDbService; if (!poolDbService.isAvailable()) @@ -130,32 +128,33 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM } //find max - int max_bin_pos = 1; - for (int i = 0; i < hists.toT[chid]->getNbinsX(); i++) { + int max_bin_pos=1; + for (int i = 0; igetNbinsX(); 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_) + int bin_x_pos = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(bin_value); + if(bin_x_pos > 20) break; - if (bin_value > hists.toT[chid]->getBinContent(max_bin_pos)) - max_bin_pos = i; + if(bin_value > hists.toT[chid]->getBinContent(i)) + 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_) { + double threshold = threshold_percent_of_max_ * hists.toT[chid]->getBinContent(max_bin_pos); + while(upper_limit_pos < hists.toT[chid]->getNbinsX()){ upper_limit_pos++; - if (hists.toT[chid]->getBinContent(upper_limit_pos) < threshold) + if(hists.toT[chid]->getBinContent(upper_limit_pos)getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos) > lower_limit_range_search_) { + while(lower_limit_pos > 0){ lower_limit_pos--; - if (hists.toT[chid]->getBinContent(lower_limit_pos) < threshold) + if(hists.toT[chid]->getBinContent(upper_limit_pos)getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos); double lower_tot_range = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos); + //const double upper_tot_range = hists.toT[chid]->getMean() + 2.5; { // scope for x-profile std::string ch_name; @@ -202,8 +201,7 @@ void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescri desc.add("formula", "[0]/(exp((x-[1])/[2])+1)+[3]") ->setComment("interpolation formula for the time walk component"); 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"); + desc.add("thresholdPercentOfMax", 0.1)->setComment("threshold for TOT fit, defined as percent of max count in 1D TOT"); descriptions.addWithDefaultLabel(desc); } From 84d76f849dd95896536fce023ed4fb1c51eb7e74 Mon Sep 17 00:00:00 2001 From: Christopher Date: Sat, 20 Jan 2024 17:33:37 +0100 Subject: [PATCH 05/20] Update the limits of the calibration fit --- .../plugins/PPSTimingCalibrationPCLHarvester.cc | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index 800e71c717a55..002aa42db2302 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -131,24 +131,24 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM int max_bin_pos=1; for (int i = 0; igetNbinsX(); i++) { double bin_value = hists.toT[chid]->getBinContent(i); - int bin_x_pos = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(bin_value); + int bin_x_pos = hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(i); if(bin_x_pos > 20) break; - if(bin_value > hists.toT[chid]->getBinContent(i)) + 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_percent_of_max_ * hists.toT[chid]->getBinContent(max_bin_pos); - while(upper_limit_pos < hists.toT[chid]->getNbinsX()){ + while(hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos) < 18){ upper_limit_pos++; if(hists.toT[chid]->getBinContent(upper_limit_pos) 0){ + while(hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos) > 8){ lower_limit_pos--; - if(hists.toT[chid]->getBinContent(upper_limit_pos)getBinContent(lower_limit_pos)getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos); From 1ff0ea773e4df1ac51d189ad6d0018de8c443836 Mon Sep 17 00:00:00 2001 From: Christopher Date: Wed, 13 Mar 2024 20:08:20 +0100 Subject: [PATCH 06/20] CMSSW PR changes and threshold parameter parsing --- .../PPSTimingCalibrationPCLHarvester.cc | 34 ++++++++++--------- .../plugins/PPSTimingCalibrationPCLWorker.cc | 2 +- .../test/DiamondCalibrationHarvester_cfg.py | 7 ++++ 3 files changed, 26 insertions(+), 17 deletions(-) diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index 002aa42db2302..c75e2403a7568 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -42,7 +42,10 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester { const std::string dqmDir_; const std::string formula_; const unsigned int min_entries_; - const double threshold_percent_of_max_; + 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_; @@ -55,9 +58,8 @@ PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::Pa dqmDir_(iConfig.getParameter("dqmDir")), formula_(iConfig.getParameter("formula")), min_entries_(iConfig.getParameter("minEntries")), - threshold_percent_of_max_(iConfig.getParameter("thresholdPercentOfMax")), - interp_("interp", formula_.c_str(), 10.5, 25.){ - + threshold_fraction_of_max_(iConfig.getParameter("thresholdFractionOfMax")), + interp_("interp", formula_.c_str(), 10.5, 25.) { // first ensure DB output service is available edm::Service poolDbService; if (!poolDbService.isAvailable()) @@ -128,33 +130,32 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM } //find max - int max_bin_pos=1; - for (int i = 0; igetNbinsX(); i++) { + 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 > 20) + if (bin_x_pos > upper_limit_max_search_) break; - if(bin_value > hists.toT[chid]->getBinContent(max_bin_pos)) - max_bin_pos=i; + 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_percent_of_max_ * hists.toT[chid]->getBinContent(max_bin_pos); - while(hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(upper_limit_pos) < 18){ + 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)getBinContent(upper_limit_pos) < threshold) break; } - while(hists.toT[chid]->getTH1()->GetXaxis()->GetBinCenter(lower_limit_pos) > 8){ + 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)getBinContent(lower_limit_pos) < threshold) break; } 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); - //const double upper_tot_range = hists.toT[chid]->getMean() + 2.5; { // scope for x-profile std::string ch_name; @@ -201,7 +202,8 @@ void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescri desc.add("formula", "[0]/(exp((x-[1])/[2])+1)+[3]") ->setComment("interpolation formula for the time walk component"); desc.add("minEntries", 100)->setComment("minimal number of hits to extract calibration"); - desc.add("thresholdPercentOfMax", 0.1)->setComment("threshold for TOT fit, defined as percent of max count in 1D TOT"); + 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..a858617c505bf 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc @@ -82,7 +82,7 @@ void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, 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.toT[detid.rawId()] = iBooker.book1D("tot_" + ch_name, ch_name + ";ToT (ns);Entries", 160, -20., 20.); iHists.leadingTimeVsToT[detid.rawId()] = iBooker.book2D("tvstot_" + ch_name, ch_name + ";ToT (ns);t (ns)", 240, 0., 60., 450, -20., 25.); } diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py index 8ad3f70a7fb4a..6d7b13baafc92 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py @@ -14,6 +14,12 @@ '', VarParsing.VarParsing.multiplicity.list, VarParsing.VarParsing.varType.string) + +options.register('thresholdFractionOfMax', + '', + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.float) + options.parseArguments() process.load("FWCore.MessageService.MessageLogger_cfi") @@ -50,6 +56,7 @@ process.load("CalibPPS.TimingCalibration.ppsTimingCalibrationPCLHarvester_cfi") #process.PPSDiamondSampicTimingCalibrationPCLHarvester.jsonCalibFile="initial.cal.json" +process.ppsTimingCalibrationPCLHarvester.thresholdFractionOfMax = options.thresholdFractionOfMax # load DQM framework process.load("DQMServices.Core.DQMStore_cfi") From f2c8b0f0c2a6783c74e9f02473d9cba64ab373fb Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin <92521351+tostafin@users.noreply.github.com> Date: Mon, 15 Jul 2024 16:13:25 +0200 Subject: [PATCH 07/20] Add LumiList filter from JSON for worker (#88) Co-authored-by: Tomasz Ostafin --- .../test/DiamondCalibrationWorker_cfg.py | 18 +++++++++++++++--- 1 file changed, 15 insertions(+), 3 deletions(-) diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py index 80930401c4e64..f81dce351c90c 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py @@ -5,15 +5,22 @@ 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('jsonFileName', + '', + VarParsing.VarParsing.multiplicity.singleton, + VarParsing.VarParsing.varType.string, + "JSON file list name") + options.parseArguments() @@ -52,6 +59,11 @@ process.source = cms.Source("PoolSource", fileNames = cms.untracked.vstring (options.inputFiles)) +if options.jsonFileName != '': + import FWCore.PythonUtilities.LumiList as LumiList + print(f"Using JSON file: {options.jsonFileName}") + process.source.lumisToProcess = LumiList.LumiList(filename=options.jsonFileName).getVLuminosityBlockRange() + process.load("CalibPPS.TimingCalibration.ppsTimingCalibrationPCLWorker_cfi") process.DQMStore = cms.Service("DQMStore") @@ -70,7 +82,7 @@ #process.totemTimingRawToDigi.rawDataTag = "hltPPSCalibrationRaw" process.path = cms.Path( - #process.a1* + #process.a1* #process.ctppsRawToDigi * #process.recoCTPPS * process.ppsTimingCalibrationPCLWorker * From f5fd4818dea9ae29c1dbc65ed4156dbc7ebc0ea9 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Sat, 9 Nov 2024 01:04:54 +0100 Subject: [PATCH 08/20] Perform t vs TOT fitting in iterations The thresholdFractionOfMax parameter is removed and replaced by iterative thresholds plus fixed bounds to find the best fit based on chi^2 and ndf. --- .../PPSTimingCalibrationPCLHarvester.cc | 339 +++++++++++------- .../test/DiamondCalibrationHarvester_cfg.py | 86 +++-- 2 files changed, 252 insertions(+), 173 deletions(-) diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index c75e2403a7568..df751a5550a62 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -5,69 +5,74 @@ * Edoardo Bossini * Piotr Maciej Cwiklicki * Laurent Forthomme + * Tomasz Ostafin * ****************************************************************************/ +#include "CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.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&, + TimingCalibrationHistograms&, + const CTPPSDiamondDetId&, + const uint32_t, + const std::string&) const; + TProfile* createTVsTotProfile(DQMStore::IBooker&, dqm::reco::MonitorElement*, const std::string&) const; + std::pair findFitRange(dqm::reco::MonitorElement*, const double, const double) const; + 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_; + 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")}, + 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 +80,202 @@ 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); + } } } //------------------------------------------------------------------------------ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) { + 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); + // book the parameters containers - PPSTimingCalibration::ParametersMap calib_params; - PPSTimingCalibration::TimingMap calib_time; + PPSTimingCalibration::ParametersMap calibParams; + PPSTimingCalibration::TimingMap calibTime; iGetter.cd(); iGetter.setCurrentFolder(dqmDir_); + constexpr double defaultFitSlope{0.8}; + 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); + if (fetchWorkerHistograms(iGetter, hists, detId, channelId, channelName)) { + const PPSTimingCalibration::Key armKey{static_cast(detId.arm()), + static_cast(detId.station()), + static_cast(detId.plane()), + static_cast(detId.channel())}; - //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; + TProfile* const tVsTotProfile{createTVsTotProfile(iBooker, hists.leadingTimeVsToT.at(channelId), channelName)}; + + const double defaultUpperLowerAsymptotesDiff = hists.leadingTime[channelId]->getRMS(); + const double defaultCenterOfDistribution = hists.toT[channelId]->getMean(); + const double defaultLowerAsymptote = hists.leadingTime[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) { + interp.SetParameters( + defaultUpperLowerAsymptotesDiff, defaultCenterOfDistribution, defaultFitSlope, defaultLowerAsymptote); + const auto [lowerTotRange, upperTotRange] = + findFitRange(hists.toT.at(channelId), 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; + } + } + } + } + + calibParams[armKey] = {0.0, 0.0, 0.0, 0.0}; + calibTime[armKey] = {defaultOffset, defaultResolution}; + if (bestUpperTotRange != 0.0) { + tVsTotProfile->Fit(&interp, "B", "", bestLowerTotRange, bestUpperTotRange); + calibParams[armKey] = { + interp.GetParameter(0), interp.GetParameter(1), interp.GetParameter(2), interp.GetParameter(3)}; + calibTime[armKey] = {defaultOffset, + defaultResolution}; // hardcoded offset/resolution placeholder for the time being + } 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, + TimingCalibrationHistograms& hists, + const CTPPSDiamondDetId& detId, + const uint32_t channelId, + const std::string& channelName) const { + hists.leadingTime[channelId] = iGetter.get(dqmDir_ + "/t_" + channelName); + if (!hists.leadingTime.at(channelId)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve leading time monitor for channel (" << detId << "). Skipping calibration."; + return false; + } + + hists.toT[channelId] = iGetter.get(dqmDir_ + "/tot_" + channelName); + if (!hists.toT.at(channelId)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve time over threshold monitor for channel (" << detId << "). Skipping calibration."; + return false; + } + + hists.leadingTimeVsToT[channelId] = iGetter.get(dqmDir_ + "/tvstot_" + channelName); + if (!hists.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(hists.leadingTimeVsToT.at(channelId)->getEntries()); + if (tVsTotEntries < minEntries_) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Not enough entries for channel (" << detId << "): " << tVsTotEntries << " < " << minEntries_ + << ". Skipping calibration."; + return false; + } + + 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( + dqm::reco::MonitorElement* tot, + const double lowerThresholdFractionOfMax, + const double upperThresholdFractionOfMax) const { + constexpr double totUpperLimitMaxSearch{20.0}; + int maxTotBin{1}; + const int numOfToTBins{tot->getNbinsX()}; + const TAxis* const totXAxis{tot->getTH1()->GetXaxis()}; + for (int i{2}; i <= numOfToTBins || totXAxis->GetBinCenter(i) <= totUpperLimitMaxSearch; ++i) { + if (tot->getBinContent(i) > tot->getBinContent(maxTotBin)) { + maxTotBin = i; } - 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 << ")."; + } + + constexpr double totLowerLimitRangeSearch{8.0}; + double lowerTotRange{8}; + if (lowerThresholdFractionOfMax != FixedFitBoundIndication_) { + int lowerLimitPos{maxTotBin}; + const double lowerThreshold{lowerThresholdFractionOfMax * tot->getBinContent(maxTotBin)}; + 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 * tot->getBinContent(maxTotBin)}; + 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}; } //------------------------------------------------------------------------------ @@ -202,8 +287,6 @@ void PPSTimingCalibrationPCLHarvester::fillDescriptions(edm::ConfigurationDescri desc.add("formula", "[0]/(exp((x-[1])/[2])+1)+[3]") ->setComment("interpolation formula for the time walk component"); 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/test/DiamondCalibrationHarvester_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py index 6d7b13baafc92..63f34e5983a68 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py @@ -1,69 +1,71 @@ import FWCore.ParameterSet.Config as cms -from Configuration.StandardSequences.Eras import eras 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('thresholdFractionOfMax', - '', - VarParsing.VarParsing.multiplicity.singleton, - VarParsing.VarParsing.varType.float) - +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.parseArguments() process.load("FWCore.MessageService.MessageLogger_cfi") -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, options.globalTag,'') +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] - +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(fileList)) + fileNames = cms.untracked.vstring(fileList) +) # output service for database -process.load('CondCore.CondDB.CondDB_cfi') -process.CondDB.connect = 'sqlite_file:ppsDiamondTiming_calibration.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.thresholdFractionOfMax = options.thresholdFractionOfMax # 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_13_3_0" +process.dqmSaver.convention = "Offline" +process.dqmSaver.workflow = "/CalibPPS/TimingCalibration/CMSSW_14_1_4" process.dqmSaver.saveByRun = -1 process.dqmSaver.saveAtJobEnd = True @@ -78,11 +80,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 ) @@ -95,6 +94,3 @@ process.p, process.end_path ) - - - From b551352ccf0a0d82ad78117ad5f526070ac89768 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Sat, 9 Nov 2024 23:38:15 +0100 Subject: [PATCH 09/20] Add correction for leading edge double peak Adding implementation for correcting double peaks which sometimes occur in the leading edge histograms. The algorithm finds the LS at which the double peak happens, computes the required time offset and passes this information to the worker to correct t and t vs TOT histograms, as well as to the harvester to save it into the SQLite file. --- CalibPPS/TimingCalibration/BuildFile.xml | 6 + .../interface/DoublePeakCorrection.h | 33 +++ .../TimingCalibration/interface/PlaneMap.h | 17 ++ .../interface/TimingCalibrationData.h | 13 ++ .../interface/TimingCalibrationHistograms.h | 19 ++ .../interface/TimingCalibrationStruct.h | 17 -- .../TimingCalibration/plugins/BuildFile.xml | 7 +- ...mondSampicTimingCalibrationPCLHarvester.cc | 1 - .../PPSTimingCalibrationPCLHarvester.cc | 142 ++++++++----- .../plugins/PPSTimingCalibrationPCLWorker.cc | 192 +++++++++++------- .../src/DoublePeakCorrection.cc | 127 ++++++++++++ .../test/DiamondCalibrationHarvester_cfg.py | 8 + .../test/DiamondCalibrationWorker_cfg.py | 102 +++++----- 13 files changed, 479 insertions(+), 205 deletions(-) create mode 100644 CalibPPS/TimingCalibration/BuildFile.xml create mode 100644 CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h create mode 100644 CalibPPS/TimingCalibration/interface/PlaneMap.h create mode 100644 CalibPPS/TimingCalibration/interface/TimingCalibrationData.h create mode 100644 CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h delete mode 100644 CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h create mode 100644 CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc 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..c75f675062276 --- /dev/null +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -0,0 +1,33 @@ +#ifndef CalibPPS_TimingCalibration_DoublePeakCorrection_h +#define CalibPPS_TimingCalibration_DoublePeakCorrection_h + +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "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&); + bool isCorrectionNeeded(const TH2F* tVsLs, const PlaneKey& planeKey) const; + double getCorrectedLeadingTime(const double, const unsigned int, const PlaneKey&) const; + double getEncodedLsAndTimeOffset(const PlaneKey& planeKey) const; + +private: + const TH2F* getTVsLs(TFile&, const std::string&, const CTPPSDiamondDetId&); + void fillLsAndTimeOffset(const TH2F*, const PlaneKey&); + std::tuple findLsAndTimePeaks(const TH2F* tVsLs, const PlaneKey& planeKey) const; + double findTimeOffset(const TH2F*, const double, const double) const; + double findGaussianMean(const std::unique_ptr& tProjection, const double estimatedMean) const; + + std::unordered_map, PlaneKeyHash> lsAndTimeOffsets_; + + static constexpr double TMaxDiff_ = 1.5; +}; + +#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..9d0626de1d04a --- /dev/null +++ b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h @@ -0,0 +1,13 @@ +#ifndef CalibPPS_TimingCalibration_TimingCalibrationData_h +#define CalibPPS_TimingCalibration_TimingCalibrationData_h + +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "DoublePeakCorrection.h" +// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h" +#include "TimingCalibrationHistograms.h" + +struct TimingCalibrationData : TimingCalibrationHistograms { + DoublePeakCorrection doublePeakCorrection; +}; + +#endif diff --git a/CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h b/CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h new file mode 100644 index 0000000000000..564a3d0716da0 --- /dev/null +++ b/CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h @@ -0,0 +1,19 @@ +#ifndef CalibPPS_TimingCalibration_TimingCalibrationHistograms_h +#define CalibPPS_TimingCalibration_TimingCalibrationHistograms_h + +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "PlaneMap.h" + +#include "DQMServices/Core/interface/MonitorElement.h" + +struct TimingCalibrationHistograms { + using PlaneMonitorMap = std::unordered_map; + using ChannelMonitorMap = std::unordered_map; + + ChannelMonitorMap leadingTime; + ChannelMonitorMap toT; + ChannelMonitorMap leadingTimeVsToT; + PlaneMonitorMap leadingTimeVsLs; +}; + +#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 df751a5550a62..2e7c9e0ab5294 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -9,7 +9,12 @@ * ****************************************************************************/ -#include "CalibPPS/TimingCalibration/interface/TimingCalibrationStruct.h" +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "../interface/DoublePeakCorrection.h" +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "../interface/PlaneMap.h" +// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h" +#include "../interface/TimingCalibrationHistograms.h" #include "CondCore/DBOutputService/interface/PoolDBOutputService.h" @@ -47,13 +52,16 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester { bool fetchWorkerHistograms(DQMStore::IGetter&, TimingCalibrationHistograms&, const CTPPSDiamondDetId&, + const PlaneKey&, const uint32_t, const std::string&) const; TProfile* createTVsTotProfile(DQMStore::IBooker&, dqm::reco::MonitorElement*, const std::string&) const; std::pair findFitRange(dqm::reco::MonitorElement*, const double, const double) const; + DoublePeakCorrection doublePeakCorrection_; const std::string dqmDir_; const std::string formula_; + const std::string tVsLsFilename_; std::vector detIds_; const edm::ESGetToken geomEsToken_; const unsigned int minEntries_; @@ -66,6 +74,7 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester { PPSTimingCalibrationPCLHarvester::PPSTimingCalibrationPCLHarvester(const edm::ParameterSet& iConfig) : 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 @@ -85,11 +94,16 @@ void PPSTimingCalibrationPCLHarvester::beginRun(const edm::Run& iRun, const edm: detIds_.push_back(detId); } } + doublePeakCorrection_.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds_); } //------------------------------------------------------------------------------ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) { + // book the parameters containers + 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); @@ -99,10 +113,6 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM // set a higher max function calls limit for the ROOT fit algorithm ROOT::Math::MinimizerOptions::SetDefaultMaxFunctionCalls(5'000); - // book the parameters containers - PPSTimingCalibration::ParametersMap calibParams; - PPSTimingCalibration::TimingMap calibTime; - iGetter.cd(); iGetter.setCurrentFolder(dqmDir_); @@ -113,57 +123,62 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM {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; + TimingCalibrationHistograms histograms; std::string channelName; for (const auto& detId : detIds_) { const uint32_t channelId{detId.rawId()}; detId.channelName(channelName); - if (fetchWorkerHistograms(iGetter, hists, detId, channelId, channelName)) { - const PPSTimingCalibration::Key armKey{static_cast(detId.arm()), - static_cast(detId.station()), - static_cast(detId.plane()), - static_cast(detId.channel())}; - - TProfile* const tVsTotProfile{createTVsTotProfile(iBooker, hists.leadingTimeVsToT.at(channelId), channelName)}; - - const double defaultUpperLowerAsymptotesDiff = hists.leadingTime[channelId]->getRMS(); - const double defaultCenterOfDistribution = hists.toT[channelId]->getMean(); - const double defaultLowerAsymptote = hists.leadingTime[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) { - interp.SetParameters( - defaultUpperLowerAsymptotesDiff, defaultCenterOfDistribution, defaultFitSlope, defaultLowerAsymptote); - const auto [lowerTotRange, upperTotRange] = - findFitRange(hists.toT.at(channelId), 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; + 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, histograms, detId, planeKey, channelId, channelName)) { + if (tVsLsFilename_.empty() && doublePeakCorrection_.isCorrectionNeeded(histograms.leadingTimeVsLs[planeKey]->getTH2F(), planeKey)) { + calibTime[channelKey] = {0.0, 0.0}; + } else { + TProfile* const tVsTotProfile{ + createTVsTotProfile(iBooker, histograms.leadingTimeVsToT.at(channelId), channelName)}; + + const double defaultUpperLowerAsymptotesDiff = histograms.leadingTime.at(channelId)->getRMS(); + const double defaultCenterOfDistribution = histograms.toT.at(channelId)->getMean(); + const double defaultLowerAsymptote = + histograms.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) { + interp.SetParameters( + defaultUpperLowerAsymptotesDiff, defaultCenterOfDistribution, defaultFitSlope, defaultLowerAsymptote); + const auto [lowerTotRange, upperTotRange] = + findFitRange(histograms.toT.at(channelId), 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; + } } } } - } - calibParams[armKey] = {0.0, 0.0, 0.0, 0.0}; - calibTime[armKey] = {defaultOffset, defaultResolution}; - if (bestUpperTotRange != 0.0) { - tVsTotProfile->Fit(&interp, "B", "", bestLowerTotRange, bestUpperTotRange); - calibParams[armKey] = { - interp.GetParameter(0), interp.GetParameter(1), interp.GetParameter(2), interp.GetParameter(3)}; - calibTime[armKey] = {defaultOffset, - defaultResolution}; // hardcoded offset/resolution placeholder for the time being - } else { - edm::LogWarning{"PPSTimingCalibrationPCLHarvester:dqmEndJob"} << "Fit did not converge for channel (" << detId - << ")."; + 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] = {doublePeakCorrection_.getEncodedLsAndTimeOffset(planeKey), defaultResolution}; + } else { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:dqmEndJob"} << "Fit did not converge for channel (" << detId + << ")."; + } } } @@ -179,33 +194,34 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM //------------------------------------------------------------------------------ bool PPSTimingCalibrationPCLHarvester::fetchWorkerHistograms(DQMStore::IGetter& iGetter, - TimingCalibrationHistograms& hists, + TimingCalibrationHistograms& histograms, const CTPPSDiamondDetId& detId, + const PlaneKey& planeKey, const uint32_t channelId, const std::string& channelName) const { - hists.leadingTime[channelId] = iGetter.get(dqmDir_ + "/t_" + channelName); - if (!hists.leadingTime.at(channelId)) { + histograms.leadingTime[channelId] = iGetter.get(dqmDir_ + "/t_" + channelName); + if (!histograms.leadingTime.at(channelId)) { edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} << "Failed to retrieve leading time monitor for channel (" << detId << "). Skipping calibration."; return false; } - hists.toT[channelId] = iGetter.get(dqmDir_ + "/tot_" + channelName); - if (!hists.toT.at(channelId)) { + histograms.toT[channelId] = iGetter.get(dqmDir_ + "/tot_" + channelName); + if (!histograms.toT.at(channelId)) { edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} << "Failed to retrieve time over threshold monitor for channel (" << detId << "). Skipping calibration."; return false; } - hists.leadingTimeVsToT[channelId] = iGetter.get(dqmDir_ + "/tvstot_" + channelName); - if (!hists.leadingTimeVsToT.at(channelId)) { + histograms.leadingTimeVsToT[channelId] = iGetter.get(dqmDir_ + "/tvstot_" + channelName); + if (!histograms.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(hists.leadingTimeVsToT.at(channelId)->getEntries()); + const auto tVsTotEntries = static_cast(histograms.leadingTimeVsToT.at(channelId)->getEntries()); if (tVsTotEntries < minEntries_) { edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} << "Not enough entries for channel (" << detId << "): " << tVsTotEntries << " < " << minEntries_ @@ -213,6 +229,18 @@ bool PPSTimingCalibrationPCLHarvester::fetchWorkerHistograms(DQMStore::IGetter& return false; } + if (!histograms.leadingTimeVsLs.contains(planeKey)) { + std::string planeName; + detId.planeName(planeName); + histograms.leadingTimeVsLs[planeKey] = iGetter.get(dqmDir_ + "/tvsls_" + planeName); + if (!histograms.leadingTimeVsLs.at(planeKey)) { + edm::LogWarning{"PPSTimingCalibrationPCLHarvester:fetchWorkerHistograms"} + << "Failed to retrieve leading time vs. time over threshold monitor for plane (" << planeName + << "). Skipping calibration."; + return false; + } + } + return true; } @@ -286,6 +314,8 @@ 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"); descriptions.addWithDefaultLabel(desc); } diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc index a858617c505bf..34c8bed4f8204 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc @@ -5,63 +5,89 @@ * Edoardo Bossini * Piotr Maciej Cwiklicki * Laurent Forthomme + * Tomasz Ostafin * ****************************************************************************/ +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "../interface/DoublePeakCorrection.h" +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "../interface/PlaneMap.h" +// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" +#include "../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& tcData) 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); + } + } + tcData.doublePeakCorrection.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds); } //------------------------------------------------------------------------------ @@ -69,22 +95,31 @@ PPSTimingCalibrationPCLWorker::PPSTimingCalibrationPCLWorker(const edm::Paramete void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, const edm::Run& iRun, const edm::EventSetup& iSetup, - TimingCalibrationHistograms& iHists) const { + TimingCalibrationData& tcData) 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", 160, -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()}; + + tcData.leadingTime[channelId] = + iBooker.book1D("t_" + channelName, channelName + ";t (ns);Entries", 1200, -60.0, 60.0); + tcData.toT[channelId] = iBooker.book1D("tot_" + channelName, channelName + ";ToT (ns);Entries", 160, -20.0, 20.0); + tcData.leadingTimeVsToT[channelId] = + iBooker.book2D("tvstot_" + channelName, channelName + ";ToT (ns);t (ns)", 240, 0.0, 60.0, 450, -20.0, 25.0); + if (!tcData.leadingTimeVsLs.contains(planeKey)) { + detId.planeName(planeName); + tcData.leadingTimeVsLs[planeKey] = + iBooker.book2D("tvsls_" + planeName, planeName + ";LS;t (ns)", 3000, 0.0, 3000.0, 500, 0.0, 20.0); + } + } } } @@ -92,64 +127,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& tcData) 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 (!tcData.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{ + tcData.doublePeakCorrection.getCorrectedLeadingTime(rechit.time(), ls, planeKey)}; + tcData.leadingTime.at(channelId)->Fill(correctedLeadingTime); + tcData.toT.at(channelId)->Fill(rechit.toT()); + tcData.leadingTimeVsToT.at(channelId)->Fill(rechit.toT(), correctedLeadingTime); + tcData.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..ee6709f308baa --- /dev/null +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -0,0 +1,127 @@ +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "../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()}; + if (!lsAndTimeOffsets_.contains(planeKey)) { + 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) { + const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); + + if (doublePeakLs != 1) { + lsAndTimeOffsets_[planeKey] = {doublePeakLs, + findTimeOffset(tVsLs, 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}; + double prevTWithMaxCount{0.0}; + 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; + } + } + + if (prevTWithMaxCount != 0.0 && secondPeakTWithMaxCount - prevTWithMaxCount > TMaxDiff_) { + return {lsBin, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; + } + + firstPeakTWithMaxCount = prevTWithMaxCount; + prevTWithMaxCount = secondPeakTWithMaxCount; + } + + return {1, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; +} + +double DoublePeakCorrection::findTimeOffset(const TH2F* tVsLs, + const double firstPeakEstimatedMean, + const double secondPeakEstimatedMean) const { + const std::unique_ptr tProjection{tVsLs->ProjectionY()}; + return findGaussianMean(tProjection, secondPeakEstimatedMean) - findGaussianMean(tProjection, firstPeakEstimatedMean); +} + +double DoublePeakCorrection::findGaussianMean(const std::unique_ptr& tProjection, + const double estimatedMean) const { + constexpr unsigned int meanParamIndex{1}; + const double fitLeftBound{estimatedMean - TMaxDiff_}; + const double fitRightBound{estimatedMean + TMaxDiff_}; + TF1 fitFunction{"peak", "gaus"}; + fitFunction.SetParLimits(meanParamIndex, fitLeftBound, fitRightBound); + fitFunction.SetParameter(meanParamIndex, estimatedMean); + const TFitResultPtr& peakFit{tProjection->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 TH2F* tVsLs, const PlaneKey& planeKey) const { + const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); + return doublePeakLs != 1; +} + +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_)) { + constexpr double encodingMultiple = 100'000.0; + const auto [doublePeakLs, doublePeakTimeOffset] = it->second; + return doublePeakLs * encodingMultiple + doublePeakTimeOffset; + } + return 0.0; +} diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py index 63f34e5983a68..bd8da12600723 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationHarvester_cfg.py @@ -20,10 +20,17 @@ 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)) @@ -59,6 +66,7 @@ ) process.load("CalibPPS.TimingCalibration.ppsTimingCalibrationPCLHarvester_cfi") +process.ppsTimingCalibrationPCLHarvester.tVsLsFilename = options.tVsLsFilename # load DQM framework process.load("DQMServices.Core.DQMStore_cfi") diff --git a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py index f81dce351c90c..a454f20396144 100644 --- a/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py +++ b/CalibPPS/TimingCalibration/test/DiamondCalibrationWorker_cfg.py @@ -1,90 +1,80 @@ import FWCore.ParameterSet.Config as cms import FWCore.ParameterSet.VarParsing as VarParsing +import FWCore.PythonUtilities.LumiList as LumiList -process = cms.Process("worker") +from Configuration.AlCa.GlobalTag import GlobalTag -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) +process = cms.Process("worker") -options.register('jsonFileName', - '', - VarParsing.VarParsing.multiplicity.singleton, - VarParsing.VarParsing.varType.string, - "JSON file list name") +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(-1) ) +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.GlobalTag = GlobalTag(process.GlobalTag, options.globalTag,'') -#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 (options.inputFiles)) + fileNames = cms.untracked.vstring(options.inputFiles) +) -if options.jsonFileName != '': - import FWCore.PythonUtilities.LumiList as LumiList - print(f"Using JSON file: {options.jsonFileName}") +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("PoolOutputModule", fileName = cms.untracked.string("file:worker_output.root"), outputCommands = cms.untracked.vstring( - 'drop *', - 'keep *_MEtoEDMConvertPPSTimingCalib_*_*', - )) + "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.MEtoEDMConvertPPSTimingCalib ) From 57b687596d2a65693fbb6f3e2189596a87aab918 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Thu, 14 Nov 2024 10:49:32 +0100 Subject: [PATCH 10/20] Cache DP correction results --- .../interface/DoublePeakCorrection.h | 7 ++--- .../PPSTimingCalibrationPCLHarvester.cc | 12 ++++---- .../plugins/PPSTimingCalibrationPCLWorker.cc | 11 +++---- .../src/DoublePeakCorrection.cc | 29 ++++++++----------- 4 files changed, 24 insertions(+), 35 deletions(-) diff --git a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h index c75f675062276..b51122b1947fe 100644 --- a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -1,8 +1,7 @@ #ifndef CalibPPS_TimingCalibration_DoublePeakCorrection_h #define CalibPPS_TimingCalibration_DoublePeakCorrection_h -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "PlaneMap.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" @@ -14,13 +13,13 @@ class DoublePeakCorrection { public: void extractLsAndTimeOffset(const std::string&, const unsigned int, const std::vector&); - bool isCorrectionNeeded(const TH2F* tVsLs, const PlaneKey& planeKey) const; + 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& planeKey) const; private: const TH2F* getTVsLs(TFile&, const std::string&, const CTPPSDiamondDetId&); - void fillLsAndTimeOffset(const TH2F*, const PlaneKey&); std::tuple findLsAndTimePeaks(const TH2F* tVsLs, const PlaneKey& planeKey) const; double findTimeOffset(const TH2F*, const double, const double) const; double findGaussianMean(const std::unique_ptr& tProjection, const double estimatedMean) const; diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index 2e7c9e0ab5294..ef79ce0789a13 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -9,12 +9,9 @@ * ****************************************************************************/ -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -#include "../interface/DoublePeakCorrection.h" -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "../interface/PlaneMap.h" -// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h" -#include "../interface/TimingCalibrationHistograms.h" +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h" #include "CondCore/DBOutputService/interface/PoolDBOutputService.h" @@ -136,7 +133,8 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM calibParams[channelKey] = {0.0, 0.0, 0.0, 0.0}; calibTime[channelKey] = {defaultOffset, defaultResolution}; if (fetchWorkerHistograms(iGetter, histograms, detId, planeKey, channelId, channelName)) { - if (tVsLsFilename_.empty() && doublePeakCorrection_.isCorrectionNeeded(histograms.leadingTimeVsLs[planeKey]->getTH2F(), planeKey)) { + doublePeakCorrection_.fillLsAndTimeOffset(histograms.leadingTimeVsLs[planeKey]->getTH2F(), planeKey); + if (tVsLsFilename_.empty() && doublePeakCorrection_.isCorrectionNeeded(planeKey)) { calibTime[channelKey] = {0.0, 0.0}; } else { TProfile* const tVsTotProfile{ diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc index 34c8bed4f8204..9e929db02cab2 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc @@ -9,12 +9,9 @@ * ****************************************************************************/ -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -#include "../interface/DoublePeakCorrection.h" -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "../interface/PlaneMap.h" -// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" -#include "../interface/TimingCalibrationData.h" +#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" @@ -117,7 +114,7 @@ void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, if (!tcData.leadingTimeVsLs.contains(planeKey)) { detId.planeName(planeName); tcData.leadingTimeVsLs[planeKey] = - iBooker.book2D("tvsls_" + planeName, planeName + ";LS;t (ns)", 3000, 0.0, 3000.0, 500, 0.0, 20.0); + iBooker.book2D("tvsls_" + planeName, planeName + ";LS;t (ns)", 3000, 1.0, 3001.0, 500, 0.0, 20.0); } } } diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index ee6709f308baa..c9498778684e4 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -1,5 +1,4 @@ -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -#include "../interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" #include "FWCore/Utilities/interface/EDMException.h" @@ -15,9 +14,7 @@ void DoublePeakCorrection::extractLsAndTimeOffset(const std::string& tVsLsFilena const std::string runNumber{std::to_string(run)}; for (const auto& detId : detIds) { const PlaneKey planeKey{detId.arm(), detId.station(), detId.plane()}; - if (!lsAndTimeOffsets_.contains(planeKey)) { - fillLsAndTimeOffset(getTVsLs(tVsLsFile, runNumber, detId), planeKey); - } + fillLsAndTimeOffset(getTVsLs(tVsLsFile, runNumber, detId), planeKey); } } else { throw edm::Exception{edm::errors::FileOpenError} << "Can't open the file with t vs LS: " << tVsLsFilename << '.'; @@ -40,11 +37,12 @@ const TH2F* DoublePeakCorrection::getTVsLs(TFile& tVsLsFile, } void DoublePeakCorrection::fillLsAndTimeOffset(const TH2F* tVsLs, const PlaneKey& planeKey) { - const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); - - if (doublePeakLs != 1) { - lsAndTimeOffsets_[planeKey] = {doublePeakLs, - findTimeOffset(tVsLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; + if (!lsAndTimeOffsets_.contains(planeKey)) { + const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); + if (doublePeakLs != 1) { + lsAndTimeOffsets_[planeKey] = {doublePeakLs, + findTimeOffset(tVsLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; + } } } @@ -54,7 +52,6 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak auto numOfTBins = static_cast(tVsLs->GetNbinsY()); double firstPeakTWithMaxCount{0.0}; double secondPeakTWithMaxCount{0.0}; - double prevTWithMaxCount{0.0}; for (unsigned int lsBin{1}; lsBin <= numOfLsBins; ++lsBin) { double tMaxCount{0}; for (unsigned int tBin{1}; tBin <= numOfTBins; ++tBin) { @@ -67,12 +64,11 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak } } - if (prevTWithMaxCount != 0.0 && secondPeakTWithMaxCount - prevTWithMaxCount > TMaxDiff_) { + if (firstPeakTWithMaxCount != 0.0 && std::abs(secondPeakTWithMaxCount - firstPeakTWithMaxCount) > TMaxDiff_) { return {lsBin, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; } - firstPeakTWithMaxCount = prevTWithMaxCount; - prevTWithMaxCount = secondPeakTWithMaxCount; + firstPeakTWithMaxCount = secondPeakTWithMaxCount; } return {1, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; @@ -100,9 +96,8 @@ double DoublePeakCorrection::findGaussianMean(const std::unique_ptr& tProj throw edm::Exception{edm::errors::FatalRootError} << "Double peak Gaussian fit not valid."; } -bool DoublePeakCorrection::isCorrectionNeeded(const TH2F* tVsLs, const PlaneKey& planeKey) const { - const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); - return doublePeakLs != 1; +bool DoublePeakCorrection::isCorrectionNeeded(const PlaneKey& planeKey) const { + return lsAndTimeOffsets_.contains(planeKey); } double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, From ac354733e16c85dff318ecdd865a24094fae7354 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Mon, 18 Nov 2024 18:39:13 +0100 Subject: [PATCH 11/20] Separate fitting of two t projections for lower t differences --- .../interface/DoublePeakCorrection.h | 13 ++-- .../interface/TimingCalibrationData.h | 19 ++++- .../interface/TimingCalibrationHistograms.h | 19 ----- .../PPSTimingCalibrationPCLHarvester.cc | 74 +++++++++---------- .../plugins/PPSTimingCalibrationPCLWorker.cc | 9 ++- .../src/DoublePeakCorrection.cc | 23 +++--- 6 files changed, 78 insertions(+), 79 deletions(-) delete mode 100644 CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h diff --git a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h index b51122b1947fe..f1430dc8da10e 100644 --- a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -1,7 +1,8 @@ #ifndef CalibPPS_TimingCalibration_DoublePeakCorrection_h #define CalibPPS_TimingCalibration_DoublePeakCorrection_h -#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "../interface/PlaneMap.h" #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" @@ -16,17 +17,15 @@ class DoublePeakCorrection { 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& planeKey) const; + double getEncodedLsAndTimeOffset(const PlaneKey&) const; private: const TH2F* getTVsLs(TFile&, const std::string&, const CTPPSDiamondDetId&); - std::tuple findLsAndTimePeaks(const TH2F* tVsLs, const PlaneKey& planeKey) const; - double findTimeOffset(const TH2F*, const double, const double) const; - double findGaussianMean(const std::unique_ptr& tProjection, const double estimatedMean) const; + 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_; - - static constexpr double TMaxDiff_ = 1.5; }; #endif diff --git a/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h index 9d0626de1d04a..cce40f4a5b24b 100644 --- a/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h +++ b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h @@ -1,12 +1,23 @@ #ifndef CalibPPS_TimingCalibration_TimingCalibrationData_h #define CalibPPS_TimingCalibration_TimingCalibrationData_h -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" #include "DoublePeakCorrection.h" -// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h" -#include "TimingCalibrationHistograms.h" -struct TimingCalibrationData : TimingCalibrationHistograms { +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "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; }; diff --git a/CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h b/CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h deleted file mode 100644 index 564a3d0716da0..0000000000000 --- a/CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h +++ /dev/null @@ -1,19 +0,0 @@ -#ifndef CalibPPS_TimingCalibration_TimingCalibrationHistograms_h -#define CalibPPS_TimingCalibration_TimingCalibrationHistograms_h - -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "PlaneMap.h" - -#include "DQMServices/Core/interface/MonitorElement.h" - -struct TimingCalibrationHistograms { - using PlaneMonitorMap = std::unordered_map; - using ChannelMonitorMap = std::unordered_map; - - ChannelMonitorMap leadingTime; - ChannelMonitorMap toT; - ChannelMonitorMap leadingTimeVsToT; - PlaneMonitorMap leadingTimeVsLs; -}; - -#endif diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index ef79ce0789a13..b6f191dda0c14 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -9,9 +9,12 @@ * ****************************************************************************/ -#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "CalibPPS/TimingCalibration/interface/TimingCalibrationHistograms.h" +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" +#include "../interface/DoublePeakCorrection.h" +#include "../interface/PlaneMap.h" +#include "../interface/TimingCalibrationData.h" #include "CondCore/DBOutputService/interface/PoolDBOutputService.h" @@ -46,16 +49,12 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester { private: void dqmEndJob(DQMStore::IBooker&, DQMStore::IGetter&) override; - bool fetchWorkerHistograms(DQMStore::IGetter&, - TimingCalibrationHistograms&, - const CTPPSDiamondDetId&, - const PlaneKey&, - const uint32_t, - const std::string&) const; + 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(dqm::reco::MonitorElement*, const double, const double) const; - DoublePeakCorrection doublePeakCorrection_; + TimingCalibrationData timingCalibartionData_; const std::string dqmDir_; const std::string formula_; const std::string tVsLsFilename_; @@ -91,7 +90,7 @@ void PPSTimingCalibrationPCLHarvester::beginRun(const edm::Run& iRun, const edm: detIds_.push_back(detId); } } - doublePeakCorrection_.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds_); + timingCalibartionData_.doublePeakCorrection.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds_); } //------------------------------------------------------------------------------ @@ -120,30 +119,30 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM {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 histograms; 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())}; + 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, histograms, detId, planeKey, channelId, channelName)) { - doublePeakCorrection_.fillLsAndTimeOffset(histograms.leadingTimeVsLs[planeKey]->getTH2F(), planeKey); - if (tVsLsFilename_.empty() && doublePeakCorrection_.isCorrectionNeeded(planeKey)) { + 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, histograms.leadingTimeVsToT.at(channelId), channelName)}; + createTVsTotProfile(iBooker, timingCalibartionData_.leadingTimeVsToT.at(channelId), channelName)}; - const double defaultUpperLowerAsymptotesDiff = histograms.leadingTime.at(channelId)->getRMS(); - const double defaultCenterOfDistribution = histograms.toT.at(channelId)->getMean(); + const double defaultUpperLowerAsymptotesDiff = timingCalibartionData_.leadingTime.at(channelId)->getRMS(); + const double defaultCenterOfDistribution = timingCalibartionData_.toT.at(channelId)->getMean(); const double defaultLowerAsymptote = - histograms.leadingTime.at(channelId)->getMean() - defaultUpperLowerAsymptotesDiff; + timingCalibartionData_.leadingTime.at(channelId)->getMean() - defaultUpperLowerAsymptotesDiff; double bestChiSqDivNdf{std::numeric_limits::max()}; double bestLowerTotRange{0.0}; @@ -152,8 +151,8 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM for (const double lowerThresholdFractionOfMax : thresholds) { interp.SetParameters( defaultUpperLowerAsymptotesDiff, defaultCenterOfDistribution, defaultFitSlope, defaultLowerAsymptote); - const auto [lowerTotRange, upperTotRange] = - findFitRange(histograms.toT.at(channelId), lowerThresholdFractionOfMax, upperThresholdFractionOfMax); + const auto [lowerTotRange, upperTotRange] = findFitRange( + timingCalibartionData_.toT.at(channelId), lowerThresholdFractionOfMax, upperThresholdFractionOfMax); const TFitResultPtr& tVsTotProfileFitResult{ tVsTotProfile->Fit(&interp, "BNS", "", lowerTotRange, upperTotRange)}; @@ -172,7 +171,8 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM tVsTotProfile->Fit(&interp, "B", "", bestLowerTotRange, bestUpperTotRange); calibParams[channelKey] = { interp.GetParameter(0), interp.GetParameter(1), interp.GetParameter(2), interp.GetParameter(3)}; - calibTime[channelKey] = {doublePeakCorrection_.getEncodedLsAndTimeOffset(planeKey), defaultResolution}; + calibTime[channelKey] = {timingCalibartionData_.doublePeakCorrection.getEncodedLsAndTimeOffset(planeKey), + defaultResolution}; } else { edm::LogWarning{"PPSTimingCalibrationPCLHarvester:dqmEndJob"} << "Fit did not converge for channel (" << detId << ")."; @@ -192,34 +192,34 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM //------------------------------------------------------------------------------ bool PPSTimingCalibrationPCLHarvester::fetchWorkerHistograms(DQMStore::IGetter& iGetter, - TimingCalibrationHistograms& histograms, const CTPPSDiamondDetId& detId, const PlaneKey& planeKey, const uint32_t channelId, - const std::string& channelName) const { - histograms.leadingTime[channelId] = iGetter.get(dqmDir_ + "/t_" + channelName); - if (!histograms.leadingTime.at(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; } - histograms.toT[channelId] = iGetter.get(dqmDir_ + "/tot_" + channelName); - if (!histograms.toT.at(channelId)) { + 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; } - histograms.leadingTimeVsToT[channelId] = iGetter.get(dqmDir_ + "/tvstot_" + channelName); - if (!histograms.leadingTimeVsToT.at(channelId)) { + 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(histograms.leadingTimeVsToT.at(channelId)->getEntries()); + 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_ @@ -227,11 +227,11 @@ bool PPSTimingCalibrationPCLHarvester::fetchWorkerHistograms(DQMStore::IGetter& return false; } - if (!histograms.leadingTimeVsLs.contains(planeKey)) { + if (!timingCalibartionData_.leadingTimeVsLs.contains(planeKey)) { std::string planeName; detId.planeName(planeName); - histograms.leadingTimeVsLs[planeKey] = iGetter.get(dqmDir_ + "/tvsls_" + planeName); - if (!histograms.leadingTimeVsLs.at(planeKey)) { + 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."; diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc index 9e929db02cab2..a93d6739e293b 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc @@ -9,9 +9,12 @@ * ****************************************************************************/ -#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" +#include "../interface/DoublePeakCorrection.h" +#include "../interface/PlaneMap.h" +#include "../interface/TimingCalibrationData.h" #include "DataFormats/Common/interface/DetSetVector.h" #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index c9498778684e4..3bcc1fcb1b7e4 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -1,4 +1,5 @@ -#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "../interface/DoublePeakCorrection.h" #include "FWCore/Utilities/interface/EDMException.h" @@ -41,7 +42,7 @@ void DoublePeakCorrection::fillLsAndTimeOffset(const TH2F* tVsLs, const PlaneKey const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); if (doublePeakLs != 1) { lsAndTimeOffsets_[planeKey] = {doublePeakLs, - findTimeOffset(tVsLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; + findTimeOffset(tVsLs, doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; } } } @@ -64,7 +65,8 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak } } - if (firstPeakTWithMaxCount != 0.0 && std::abs(secondPeakTWithMaxCount - firstPeakTWithMaxCount) > TMaxDiff_) { + constexpr double tMinDiff{1.0}; + if (firstPeakTWithMaxCount != 0.0 && std::abs(secondPeakTWithMaxCount - firstPeakTWithMaxCount) > tMinDiff) { return {lsBin, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; } @@ -75,21 +77,24 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak } double DoublePeakCorrection::findTimeOffset(const TH2F* tVsLs, + const unsigned int doublePeakLs, const double firstPeakEstimatedMean, const double secondPeakEstimatedMean) const { - const std::unique_ptr tProjection{tVsLs->ProjectionY()}; - return findGaussianMean(tProjection, secondPeakEstimatedMean) - findGaussianMean(tProjection, firstPeakEstimatedMean); + 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& tProjection, +double DoublePeakCorrection::findGaussianMean(const std::unique_ptr& peak, const double estimatedMean) const { constexpr unsigned int meanParamIndex{1}; - const double fitLeftBound{estimatedMean - TMaxDiff_}; - const double fitRightBound{estimatedMean + TMaxDiff_}; + 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{tProjection->Fit(&fitFunction, "NS", "", fitLeftBound, fitRightBound)}; + const TFitResultPtr& peakFit{peak->Fit(&fitFunction, "NS", "", fitLeftBound, fitRightBound)}; if (peakFit->IsValid()) { return peakFit->Parameter(meanParamIndex); } From 940ef1e9c729cd0f39ed42915660a6a5c144cec9 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Mon, 18 Nov 2024 19:27:51 +0100 Subject: [PATCH 12/20] Remove relative includes --- .../interface/DoublePeakCorrection.h | 3 +- .../interface/TimingCalibrationData.h | 7 +--- .../PPSTimingCalibrationPCLHarvester.cc | 9 ++--- .../plugins/PPSTimingCalibrationPCLWorker.cc | 40 +++++++++---------- .../src/DoublePeakCorrection.cc | 10 ++--- 5 files changed, 29 insertions(+), 40 deletions(-) diff --git a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h index f1430dc8da10e..33b0573cb518e 100644 --- a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -1,8 +1,7 @@ #ifndef CalibPPS_TimingCalibration_DoublePeakCorrection_h #define CalibPPS_TimingCalibration_DoublePeakCorrection_h -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "../interface/PlaneMap.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" #include "DataFormats/CTPPSDetId/interface/CTPPSDiamondDetId.h" diff --git a/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h index cce40f4a5b24b..9cfc57bc18413 100644 --- a/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h +++ b/CalibPPS/TimingCalibration/interface/TimingCalibrationData.h @@ -1,11 +1,8 @@ #ifndef CalibPPS_TimingCalibration_TimingCalibrationData_h #define CalibPPS_TimingCalibration_TimingCalibrationData_h -#include "DoublePeakCorrection.h" - -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -#include "PlaneMap.h" - +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" #include "DQMServices/Core/interface/MonitorElement.h" diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index b6f191dda0c14..3d8e573b4e104 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -9,12 +9,9 @@ * ****************************************************************************/ -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" -#include "../interface/DoublePeakCorrection.h" -#include "../interface/PlaneMap.h" -#include "../interface/TimingCalibrationData.h" +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/PlaneMap.h" +#include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" #include "CondCore/DBOutputService/interface/PoolDBOutputService.h" diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc index a93d6739e293b..1d835a360ddd5 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLWorker.cc @@ -9,12 +9,9 @@ * ****************************************************************************/ -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -// #include "CalibPPS/TimingCalibration/interface/PlaneMap.h" -// #include "CalibPPS/TimingCalibration/interface/TimingCalibrationData.h" -#include "../interface/DoublePeakCorrection.h" -#include "../interface/PlaneMap.h" -#include "../interface/TimingCalibrationData.h" +#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" @@ -78,7 +75,7 @@ PPSTimingCalibrationPCLWorker::PPSTimingCalibrationPCLWorker(const edm::Paramete void PPSTimingCalibrationPCLWorker::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup, - TimingCalibrationData& tcData) const { + TimingCalibrationData& timingCalibrationData) const { std::vector detIds; const auto& geom = iSetup.getData(geomEsToken_); for (auto it = geom.beginSensor(); it != geom.endSensor(); ++it) { @@ -87,7 +84,7 @@ void PPSTimingCalibrationPCLWorker::dqmBeginRun(const edm::Run& iRun, detIds.push_back(detId); } } - tcData.doublePeakCorrection.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds); + timingCalibrationData.doublePeakCorrection.extractLsAndTimeOffset(tVsLsFilename_, iRun.run(), detIds); } //------------------------------------------------------------------------------ @@ -95,7 +92,7 @@ void PPSTimingCalibrationPCLWorker::dqmBeginRun(const edm::Run& iRun, void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, const edm::Run& iRun, const edm::EventSetup& iSetup, - TimingCalibrationData& tcData) const { + TimingCalibrationData& timingCalibrationData) const { iBooker.cd(); iBooker.setCurrentFolder(dqmDir_); @@ -109,14 +106,15 @@ void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, detId.channelName(channelName); const PlaneKey planeKey{detId.arm(), detId.station(), detId.plane()}; - tcData.leadingTime[channelId] = + timingCalibrationData.leadingTime[channelId] = iBooker.book1D("t_" + channelName, channelName + ";t (ns);Entries", 1200, -60.0, 60.0); - tcData.toT[channelId] = iBooker.book1D("tot_" + channelName, channelName + ";ToT (ns);Entries", 160, -20.0, 20.0); - tcData.leadingTimeVsToT[channelId] = + 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 (!tcData.leadingTimeVsLs.contains(planeKey)) { + if (!timingCalibrationData.leadingTimeVsLs.contains(planeKey)) { detId.planeName(planeName); - tcData.leadingTimeVsLs[planeKey] = + timingCalibrationData.leadingTimeVsLs[planeKey] = iBooker.book2D("tvsls_" + planeName, planeName + ";LS;t (ns)", 3000, 1.0, 3001.0, 500, 0.0, 20.0); } } @@ -127,7 +125,7 @@ void PPSTimingCalibrationPCLWorker::bookHistograms(DQMStore::IBooker& iBooker, void PPSTimingCalibrationPCLWorker::dqmAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup, - const TimingCalibrationData& tcData) const { + const TimingCalibrationData& timingCalibrationData) const { edm::Handle dsvRechits; // then extract the rechits information for later processing searchForProduct(iEvent, diamondRecHitTokens_, recHitTags_, dsvRechits); @@ -144,7 +142,7 @@ void PPSTimingCalibrationPCLWorker::dqmAnalyze(const edm::Event& iEvent, for (const auto& dsRechits : *dsvRechits) { const CTPPSDiamondDetId detId{dsRechits.detId()}; const uint32_t channelId{detId.rawId()}; - if (!tcData.leadingTimeVsToT.contains(channelId)) { + if (!timingCalibrationData.leadingTimeVsToT.contains(channelId)) { edm::LogWarning{"PPSTimingCalibrationPCLWorker:dqmAnalyze"} << "Pad with Detector ID =" << detId << " is not set to be monitored."; continue; @@ -156,11 +154,11 @@ void PPSTimingCalibrationPCLWorker::dqmAnalyze(const edm::Event& iEvent, // skip invalid rechits if (rechit.time() != 0.0 && rechit.toT() >= 0.0) { const double correctedLeadingTime{ - tcData.doublePeakCorrection.getCorrectedLeadingTime(rechit.time(), ls, planeKey)}; - tcData.leadingTime.at(channelId)->Fill(correctedLeadingTime); - tcData.toT.at(channelId)->Fill(rechit.toT()); - tcData.leadingTimeVsToT.at(channelId)->Fill(rechit.toT(), correctedLeadingTime); - tcData.leadingTimeVsLs.at(planeKey)->Fill(ls, 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); } } } diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index 3bcc1fcb1b7e4..6c4f4c772624a 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -1,5 +1,4 @@ -// #include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" -#include "../interface/DoublePeakCorrection.h" +#include "CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h" #include "FWCore/Utilities/interface/EDMException.h" @@ -41,8 +40,8 @@ void DoublePeakCorrection::fillLsAndTimeOffset(const TH2F* tVsLs, const PlaneKey if (!lsAndTimeOffsets_.contains(planeKey)) { const auto [doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount] = findLsAndTimePeaks(tVsLs, planeKey); if (doublePeakLs != 1) { - lsAndTimeOffsets_[planeKey] = {doublePeakLs, - findTimeOffset(tVsLs, doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; + lsAndTimeOffsets_[planeKey] = { + doublePeakLs, findTimeOffset(tVsLs, doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount)}; } } } @@ -85,8 +84,7 @@ double DoublePeakCorrection::findTimeOffset(const TH2F* tVsLs, return findGaussianMean(secondPeak, secondPeakEstimatedMean) - findGaussianMean(firstPeak, firstPeakEstimatedMean); } -double DoublePeakCorrection::findGaussianMean(const std::unique_ptr& peak, - const double estimatedMean) const { +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}; From 89a963bddc8ebe557d731979cdbbb9ca2c307de3 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Tue, 19 Nov 2024 18:49:01 +0100 Subject: [PATCH 13/20] Check for DP in sliding windows --- .../TimingCalibration/src/DoublePeakCorrection.cc | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index 6c4f4c772624a..c90f910edc572 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -52,6 +52,8 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak 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) { @@ -66,10 +68,17 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak constexpr double tMinDiff{1.0}; if (firstPeakTWithMaxCount != 0.0 && std::abs(secondPeakTWithMaxCount - firstPeakTWithMaxCount) > tMinDiff) { - return {lsBin, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; + constexpr unsigned int minTDiffCount{3}; + ++tDiffCount; + if (tDiffCount == 1) { + doublePeakLs = lsBin; + } else if (tDiffCount == minTDiffCount) { + return {doublePeakLs, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; + } + } else { + tDiffCount = 0; + firstPeakTWithMaxCount = secondPeakTWithMaxCount; } - - firstPeakTWithMaxCount = secondPeakTWithMaxCount; } return {1, firstPeakTWithMaxCount, secondPeakTWithMaxCount}; From 6b992b97b1729712ddf7e63b9efe9cbdffda26c1 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Sat, 7 Dec 2024 11:22:51 +0100 Subject: [PATCH 14/20] Compute max ToT bin only once per channel --- .../PPSTimingCalibrationPCLHarvester.cc | 44 ++++++++----------- 1 file changed, 19 insertions(+), 25 deletions(-) diff --git a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc index 3d8e573b4e104..5c6f320db455e 100644 --- a/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc +++ b/CalibPPS/TimingCalibration/plugins/PPSTimingCalibrationPCLHarvester.cc @@ -49,7 +49,7 @@ class PPSTimingCalibrationPCLHarvester : public DQMEDHarvester { 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(dqm::reco::MonitorElement*, const double, const double) const; + std::pair findFitRange(const TH1F*, const int, const double, const double) const; TimingCalibrationData timingCalibartionData_; const std::string dqmDir_; @@ -109,7 +109,6 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM iGetter.cd(); iGetter.setCurrentFolder(dqmDir_); - constexpr double defaultFitSlope{0.8}; constexpr double defaultOffset{0.0}; constexpr double defaultResolution{0.1}; constexpr std::array thresholds{ @@ -136,20 +135,23 @@ void PPSTimingCalibrationPCLHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQM TProfile* const tVsTotProfile{ createTVsTotProfile(iBooker, timingCalibartionData_.leadingTimeVsToT.at(channelId), channelName)}; - const double defaultUpperLowerAsymptotesDiff = timingCalibartionData_.leadingTime.at(channelId)->getRMS(); - const double defaultCenterOfDistribution = timingCalibartionData_.toT.at(channelId)->getMean(); - const double defaultLowerAsymptote = - timingCalibartionData_.leadingTime.at(channelId)->getMean() - defaultUpperLowerAsymptotesDiff; + 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( - timingCalibartionData_.toT.at(channelId), lowerThresholdFractionOfMax, upperThresholdFractionOfMax); + const auto [lowerTotRange, upperTotRange] = + findFitRange(tot, maxTotBin, lowerThresholdFractionOfMax, upperThresholdFractionOfMax); const TFitResultPtr& tVsTotProfileFitResult{ tVsTotProfile->Fit(&interp, "BNS", "", lowerTotRange, upperTotRange)}; @@ -258,28 +260,20 @@ TProfile* PPSTimingCalibrationPCLHarvester::createTVsTotProfile(DQMStore::IBooke return monitorTProfile; } -//------------------------------------------------------------------------------ - std::pair PPSTimingCalibrationPCLHarvester::findFitRange( - dqm::reco::MonitorElement* tot, + const TH1F* tot, + const int maxTotBin, const double lowerThresholdFractionOfMax, const double upperThresholdFractionOfMax) const { - constexpr double totUpperLimitMaxSearch{20.0}; - int maxTotBin{1}; - const int numOfToTBins{tot->getNbinsX()}; - const TAxis* const totXAxis{tot->getTH1()->GetXaxis()}; - for (int i{2}; i <= numOfToTBins || totXAxis->GetBinCenter(i) <= totUpperLimitMaxSearch; ++i) { - if (tot->getBinContent(i) > tot->getBinContent(maxTotBin)) { - maxTotBin = i; - } - } - 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 * tot->getBinContent(maxTotBin)}; - while (tot->getBinContent(lowerLimitPos) >= lowerThreshold && + const double lowerThreshold{lowerThresholdFractionOfMax * maxTotCount}; + while (tot->GetBinContent(lowerLimitPos) >= lowerThreshold && totXAxis->GetBinCenter(lowerLimitPos) > totLowerLimitRangeSearch) { --lowerLimitPos; } @@ -290,8 +284,8 @@ std::pair PPSTimingCalibrationPCLHarvester::findFitRange( double upperTotRange{15}; if (upperThresholdFractionOfMax != FixedFitBoundIndication_) { int upperLimitPos{maxTotBin}; - const double upperThreshold{upperThresholdFractionOfMax * tot->getBinContent(maxTotBin)}; - while (tot->getBinContent(upperLimitPos) >= upperThreshold && + const double upperThreshold{upperThresholdFractionOfMax * maxTotCount}; + while (tot->GetBinContent(upperLimitPos) >= upperThreshold && totXAxis->GetBinCenter(upperLimitPos) < totUpperLimitRangeSearch) { ++upperLimitPos; } From 09e0bb66ee584e554c87f5e9a11c1d6d4d35d1b0 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Sat, 7 Dec 2024 20:18:02 +0100 Subject: [PATCH 15/20] Increase sliiding window size to 5 --- CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index c90f910edc572..36eba263ba517 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -68,7 +68,7 @@ std::tuple DoublePeakCorrection::findLsAndTimePeak constexpr double tMinDiff{1.0}; if (firstPeakTWithMaxCount != 0.0 && std::abs(secondPeakTWithMaxCount - firstPeakTWithMaxCount) > tMinDiff) { - constexpr unsigned int minTDiffCount{3}; + constexpr unsigned int minTDiffCount{5}; ++tDiffCount; if (tDiffCount == 1) { doublePeakLs = lsBin; From 063c2fd33f68fe28890acea677c64bfcbcdf5725 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Wed, 11 Dec 2024 20:58:04 +0100 Subject: [PATCH 16/20] Encode LS and TO with the same abs value --- CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index 36eba263ba517..dd0077bf0558d 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -126,9 +126,12 @@ double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, double DoublePeakCorrection::getEncodedLsAndTimeOffset(const PlaneKey& planeKey) const { if (auto it = lsAndTimeOffsets_.find(planeKey); it != std::end(lsAndTimeOffsets_)) { - constexpr double encodingMultiple = 100'000.0; const auto [doublePeakLs, doublePeakTimeOffset] = it->second; - return doublePeakLs * encodingMultiple + doublePeakTimeOffset; + constexpr unsigned int lsEncodingMultiple{100}; + if (doublePeakTimeOffset >= 0.0) { + return doublePeakLs * lsEncodingMultiple + doublePeakTimeOffset; + } + return -(doublePeakLs * lsEncodingMultiple - doublePeakTimeOffset); } return 0.0; } From 3fca71ac9628a0acbcf504ff3a1b00233cceb7a2 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Wed, 11 Dec 2024 22:01:58 +0100 Subject: [PATCH 17/20] Correct DP in RecHit algorithm --- .../interface/DoublePeakCorrection.h | 4 ++++ .../src/DoublePeakCorrection.cc | 13 +++++++++++ RecoPPS/Local/BuildFile.xml | 1 + .../CTPPSDiamondRecHitProducerAlgorithm.h | 4 +++- .../interface/TimingRecHitProducerAlgorithm.h | 2 +- .../TotemT2RecHitProducerAlgorithm.h | 4 +++- .../TotemTimingRecHitProducerAlgorithm.h | 4 +++- .../plugins/CTPPSDiamondRecHitProducer.cc | 22 +++++++++++++------ .../Local/plugins/TotemT2RecHitProducer.cc | 2 +- .../plugins/TotemTimingRecHitProducer.cc | 2 +- .../CTPPSDiamondRecHitProducerAlgorithm.cc | 19 +++++++++++----- .../src/TotemT2RecHitProducerAlgorithm.cc | 4 +++- .../src/TotemTimingRecHitProducerAlgorithm.cc | 4 +++- 13 files changed, 65 insertions(+), 20 deletions(-) diff --git a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h index 33b0573cb518e..e8fccf95471ef 100644 --- a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -18,6 +18,8 @@ class DoublePeakCorrection { 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; @@ -25,6 +27,8 @@ class DoublePeakCorrection { 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/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index dd0077bf0558d..4502383d5792d 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -135,3 +135,16 @@ double DoublePeakCorrection::getEncodedLsAndTimeOffset(const PlaneKey& planeKey) } return 0.0; } + +double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, + const unsigned int ls, + const double encodedLsAndTimeOffset) { + const unsigned int doublePeakLs = 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/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..d4feccc346ce2 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.run(), iEvent.luminosityBlock()); } 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..b93b1205546dc 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 + const double toa{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()); From 5f7321780e6d74eafb6ee7819e6ea87a07d69ae3 Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Fri, 13 Dec 2024 11:45:50 +0100 Subject: [PATCH 18/20] Correct DP only from Run 3 onwards --- .../TimingCalibration/src/DoublePeakCorrection.cc | 14 ++++++++------ .../Local/plugins/CTPPSDiamondRecHitProducer.cc | 2 +- .../src/CTPPSDiamondRecHitProducerAlgorithm.cc | 4 ++-- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index 4502383d5792d..32377491a663f 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -139,12 +139,14 @@ double DoublePeakCorrection::getEncodedLsAndTimeOffset(const PlaneKey& planeKey) double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, const unsigned int ls, const double encodedLsAndTimeOffset) { - const unsigned int doublePeakLs = std::abs(encodedLsAndTimeOffset) / LsEncodingMultiple_; - if (ls >= doublePeakLs) { - const double doublePeakTimeOffset = encodedLsAndTimeOffset >= 0 - ? encodedLsAndTimeOffset - doublePeakLs * LsEncodingMultiple_ - : encodedLsAndTimeOffset + doublePeakLs * LsEncodingMultiple_; - return leadingTime - doublePeakTimeOffset; + if (encodedLsAndTimeOffset >= 0.0) { + const unsigned int doublePeakLs = 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/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc b/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc index d4feccc346ce2..d2a27c275ee3b 100644 --- a/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc +++ b/RecoPPS/Local/plugins/CTPPSDiamondRecHitProducer.cc @@ -87,7 +87,7 @@ void CTPPSDiamondRecHitProducer::produce(edm::Event& iEvent, const edm::EventSet algo_.setCalibration(iSetup.getData(timingCalibrationToken_), iSetup.getData(timingCalibrationLUTToken_)); // produce the rechits collection - algo_.build(iSetup.getData(geometryToken_), digis, *pOut, iEvent.run(), iEvent.luminosityBlock()); + algo_.build(iSetup.getData(geometryToken_), digis, *pOut, iEvent.luminosityBlock(), isRun2_); } iEvent.put(std::move(pOut)); diff --git a/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc b/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc index b93b1205546dc..17fd550382b5c 100644 --- a/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc +++ b/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc @@ -81,8 +81,8 @@ void CTPPSDiamondRecHitProducerAlgorithm::build(const CTPPSGeometry& geom, // 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 - const double toa{DoublePeakCorrection::getCorrectedLeadingTime(t0, ls, ch_t_offset)}; + // 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, From a6cc2a9c7a59d4076b0d9c05115a0129b251faff Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Tue, 17 Dec 2024 23:31:13 +0100 Subject: [PATCH 19/20] Change DP-corrected leading time function --- .../interface/DoublePeakCorrection.h | 2 +- .../TimingCalibration/src/DoublePeakCorrection.cc | 11 +++++------ .../Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc | 2 +- 3 files changed, 7 insertions(+), 8 deletions(-) diff --git a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h index e8fccf95471ef..c230bec63b193 100644 --- a/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h +++ b/CalibPPS/TimingCalibration/interface/DoublePeakCorrection.h @@ -18,7 +18,7 @@ class DoublePeakCorrection { 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); + static double GetCorrectedLeadingTime(const double, const unsigned int, const double); private: const TH2F* getTVsLs(TFile&, const std::string&, const CTPPSDiamondDetId&); diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index 32377491a663f..42dae03c97ace 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -127,19 +127,18 @@ double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, double DoublePeakCorrection::getEncodedLsAndTimeOffset(const PlaneKey& planeKey) const { if (auto it = lsAndTimeOffsets_.find(planeKey); it != std::end(lsAndTimeOffsets_)) { const auto [doublePeakLs, doublePeakTimeOffset] = it->second; - constexpr unsigned int lsEncodingMultiple{100}; - if (doublePeakTimeOffset >= 0.0) { - return doublePeakLs * lsEncodingMultiple + doublePeakTimeOffset; + if (doublePeakTimeOffset > 0.0) { + return doublePeakLs * LsEncodingMultiple_ + doublePeakTimeOffset; } - return -(doublePeakLs * lsEncodingMultiple - doublePeakTimeOffset); + return -(doublePeakLs * LsEncodingMultiple_ - doublePeakTimeOffset); } return 0.0; } -double DoublePeakCorrection::getCorrectedLeadingTime(const double leadingTime, +double DoublePeakCorrection::GetCorrectedLeadingTime(const double leadingTime, const unsigned int ls, const double encodedLsAndTimeOffset) { - if (encodedLsAndTimeOffset >= 0.0) { + if (encodedLsAndTimeOffset != 0.0) { const unsigned int doublePeakLs = std::abs(encodedLsAndTimeOffset) / LsEncodingMultiple_; if (ls >= doublePeakLs) { const double doublePeakTimeOffset = encodedLsAndTimeOffset >= 0 diff --git a/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc b/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc index 17fd550382b5c..f4627acb110b4 100644 --- a/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc +++ b/RecoPPS/Local/src/CTPPSDiamondRecHitProducerAlgorithm.cc @@ -82,7 +82,7 @@ void CTPPSDiamondRecHitProducerAlgorithm::build(const CTPPSGeometry& geom, // 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); + const double toa = is_run_2 ? t0 : DoublePeakCorrection::GetCorrectedLeadingTime(t0, ls, ch_t_offset); rec_hits.emplace_back( // spatial information x_pos, From 0813d122009c40f4e1c8c4e194c6ec907318839c Mon Sep 17 00:00:00 2001 From: Tomasz Ostafin Date: Sat, 4 Jan 2025 12:16:36 +0100 Subject: [PATCH 20/20] Don't correct missing leading edges --- CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc index 42dae03c97ace..c4da6018d0cb8 100644 --- a/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc +++ b/CalibPPS/TimingCalibration/src/DoublePeakCorrection.cc @@ -138,8 +138,8 @@ double DoublePeakCorrection::getEncodedLsAndTimeOffset(const PlaneKey& planeKey) double DoublePeakCorrection::GetCorrectedLeadingTime(const double leadingTime, const unsigned int ls, const double encodedLsAndTimeOffset) { - if (encodedLsAndTimeOffset != 0.0) { - const unsigned int doublePeakLs = std::abs(encodedLsAndTimeOffset) / LsEncodingMultiple_; + 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_