File indexing completed on 2024-04-06 11:58:32
0001
0002
0003
0004
0005
0006
0007 #include "DQMServices/Core/interface/DQMEDHarvester.h"
0008 #include "DQMServices/Core/interface/DQMStore.h"
0009 #include "DQMServices/Core/interface/MonitorElement.h"
0010
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016
0017 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0018
0019 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0020
0021 #include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionData.h"
0022 #include "CondFormats/PPSObjects/interface/CTPPSRPAlignmentCorrectionsData.h"
0023 #include "CondFormats/DataRecord/interface/CTPPSRPAlignmentCorrectionsDataRcd.h"
0024
0025 #include "CondFormats/PPSObjects/interface/PPSAlignmentConfiguration.h"
0026 #include "CondFormats/DataRecord/interface/PPSAlignmentConfigurationRcd.h"
0027
0028 #include "CalibPPS/AlignmentGlobal/interface/utils.h"
0029
0030 #include <memory>
0031 #include <map>
0032 #include <vector>
0033 #include <string>
0034 #include <fstream>
0035 #include <iomanip>
0036 #include <cmath>
0037 #include <utility>
0038 #include <algorithm>
0039
0040 #include "TH1D.h"
0041 #include "TH2D.h"
0042 #include "TGraph.h"
0043 #include "TGraphErrors.h"
0044 #include "TF1.h"
0045 #include "TProfile.h"
0046 #include "TFile.h"
0047 #include "TKey.h"
0048 #include "TSpline.h"
0049 #include "TCanvas.h"
0050
0051
0052
0053 class PPSAlignmentHarvester : public DQMEDHarvester {
0054 public:
0055 PPSAlignmentHarvester(const edm::ParameterSet& iConfig);
0056 ~PPSAlignmentHarvester() override;
0057
0058 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0059
0060 private:
0061 void dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) override;
0062 void dqmEndRun(DQMStore::IBooker& iBooker,
0063 DQMStore::IGetter& iGetter,
0064 edm::Run const& iRun,
0065 edm::EventSetup const& iSetup) override;
0066
0067
0068 std::unique_ptr<TGraphErrors> buildGraphFromVector(const std::vector<PPSAlignmentConfiguration::PointErrors>& pv);
0069 std::unique_ptr<TGraphErrors> buildGraphFromMonitorElements(DQMStore::IGetter& iGetter,
0070 const PPSAlignmentConfiguration::RPConfig& rpc,
0071 const std::vector<MonitorElement*>& mes,
0072 const unsigned int fitProfileMinBinEntries,
0073 const unsigned int fitProfileMinNReasonable);
0074 void doMatch(DQMStore::IBooker& iBooker,
0075 const PPSAlignmentConfiguration& cfg,
0076 const PPSAlignmentConfiguration::RPConfig& rpc,
0077 TGraphErrors* g_ref,
0078 TGraphErrors* g_test,
0079 const PPSAlignmentConfiguration::SelectionRange& range_ref,
0080 const double sh_min,
0081 const double sh_max,
0082 double& sh_best,
0083 double& sh_best_unc);
0084
0085 void xAlignment(DQMStore::IBooker& iBooker,
0086 DQMStore::IGetter& iGetter,
0087 const PPSAlignmentConfiguration& cfg,
0088 const PPSAlignmentConfiguration& cfg_ref);
0089
0090 std::map<unsigned int, double> sh_x_map_;
0091
0092
0093 void xAlignmentRelative(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter, const PPSAlignmentConfiguration& cfg);
0094
0095
0096 static double findMax(const TF1* ff_fit);
0097 TH1D* buildModeGraph(DQMStore::IBooker& iBooker,
0098 const MonitorElement* h2_y_vs_x,
0099 const PPSAlignmentConfiguration& cfg,
0100 const PPSAlignmentConfiguration::RPConfig& rpc);
0101
0102 void yAlignment(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter, const PPSAlignmentConfiguration& cfg);
0103
0104
0105 static void writeCutPlot(
0106 TH2D* h, const double a, const double c, const double si, const double n_si, const std::string& label);
0107 static std::unique_ptr<TH1D> getTH1DFromTGraphErrors(TGraphErrors* graph,
0108 const std::string& title = "",
0109 const std::string& labels = "",
0110 int n = -1,
0111 double binWidth = -1.,
0112 double min = -1.);
0113
0114 CTPPSRPAlignmentCorrectionsData getLongIdResults(CTPPSRPAlignmentCorrectionsData finalResults);
0115
0116 edm::ESGetToken<PPSAlignmentConfiguration, PPSAlignmentConfigurationRcd> esTokenTest_;
0117 edm::ESGetToken<PPSAlignmentConfiguration, PPSAlignmentConfigurationRcd> esTokenReference_;
0118
0119
0120 const std::string dqmDir_;
0121 const std::vector<std::string> sequence_;
0122 const bool overwriteShX_;
0123 const bool writeSQLiteResults_;
0124 const bool xAliRelFinalSlopeFixed_;
0125 const bool yAliFinalSlopeFixed_;
0126 const std::pair<double, double> xCorrRange_;
0127 const std::pair<double, double> yCorrRange_;
0128 const unsigned int detectorId_;
0129 const unsigned int subdetectorId_;
0130 const bool debug_;
0131
0132
0133 std::unique_ptr<TFile> debugFile_;
0134 std::ofstream textResultsFile_;
0135 int seqPos = 1;
0136
0137 CTPPSRPAlignmentCorrectionsData xAliResults_;
0138
0139 CTPPSRPAlignmentCorrectionsData xAliRelResults_;
0140 CTPPSRPAlignmentCorrectionsData xAliRelResultsSlopeFixed_;
0141
0142 CTPPSRPAlignmentCorrectionsData yAliResults_;
0143 CTPPSRPAlignmentCorrectionsData yAliResultsSlopeFixed_;
0144 };
0145
0146
0147
0148 PPSAlignmentHarvester::PPSAlignmentHarvester(const edm::ParameterSet& iConfig)
0149 : esTokenTest_(esConsumes<PPSAlignmentConfiguration, PPSAlignmentConfigurationRcd, edm::Transition::EndRun>(
0150 edm::ESInputTag("", ""))),
0151 esTokenReference_(esConsumes<PPSAlignmentConfiguration, PPSAlignmentConfigurationRcd, edm::Transition::EndRun>(
0152 edm::ESInputTag("", "reference"))),
0153 dqmDir_(iConfig.getParameter<std::string>("dqm_dir")),
0154 sequence_(iConfig.getParameter<std::vector<std::string>>("sequence")),
0155 overwriteShX_(iConfig.getParameter<bool>("overwrite_sh_x")),
0156 writeSQLiteResults_(iConfig.getParameter<bool>("write_sqlite_results")),
0157 xAliRelFinalSlopeFixed_(iConfig.getParameter<bool>("x_ali_rel_final_slope_fixed")),
0158 yAliFinalSlopeFixed_(iConfig.getParameter<bool>("y_ali_final_slope_fixed")),
0159 xCorrRange_(std::make_pair(iConfig.getParameter<double>("x_corr_min") / 1000.,
0160 iConfig.getParameter<double>("x_corr_max") / 1000.)),
0161 yCorrRange_(std::make_pair(iConfig.getParameter<double>("y_corr_min") / 1000.,
0162 iConfig.getParameter<double>("y_corr_max") / 1000.)),
0163 detectorId_(iConfig.getParameter<unsigned int>("detector_id")),
0164 subdetectorId_(iConfig.getParameter<unsigned int>("subdetector_id")),
0165 debug_(iConfig.getParameter<bool>("debug")) {
0166 auto textResultsPath = iConfig.getParameter<std::string>("text_results_path");
0167 if (!textResultsPath.empty()) {
0168 textResultsFile_.open(textResultsPath, std::ios::out | std::ios::trunc);
0169 }
0170 if (debug_) {
0171 debugFile_ = std::make_unique<TFile>("debug_harvester.root", "recreate");
0172 }
0173
0174 edm::LogInfo("PPSAlignmentHarvester").log([&](auto& li) {
0175 li << "parameters:\n";
0176 li << "* dqm_dir: " << dqmDir_ << "\n";
0177 li << "* sequence:\n";
0178 for (unsigned int i = 0; i < sequence_.size(); i++) {
0179 li << " " << i + 1 << ": " << sequence_[i] << "\n";
0180 }
0181 li << "* overwrite_sh_x: " << std::boolalpha << overwriteShX_ << "\n";
0182 li << "* text_results_path: " << textResultsPath << "\n";
0183 li << "* write_sqlite_results: " << std::boolalpha << writeSQLiteResults_ << "\n";
0184 li << "* x_ali_rel_final_slope_fixed: " << std::boolalpha << xAliRelFinalSlopeFixed_ << "\n";
0185 li << "* y_ali_final_slope_fixed: " << std::boolalpha << yAliFinalSlopeFixed_ << "\n";
0186
0187 li << "* x_corr_min: " << std::fixed << xCorrRange_.first * 1000. << ", x_corr_max: " << xCorrRange_.second * 1000.
0188 << "\n";
0189
0190 li << "* y_corr_min: " << std::fixed << yCorrRange_.first * 1000. << ", y_corr_max: " << yCorrRange_.second * 1000.
0191 << "\n";
0192 li << "* detector_id: " << detectorId_ << "\n";
0193 li << "* subdetector_id: " << subdetectorId_ << "\n";
0194 li << "* debug: " << std::boolalpha << debug_;
0195 });
0196 }
0197
0198 PPSAlignmentHarvester::~PPSAlignmentHarvester() {
0199 if (textResultsFile_.is_open()) {
0200 textResultsFile_.close();
0201 }
0202 }
0203
0204 void PPSAlignmentHarvester::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0205 edm::ParameterSetDescription desc;
0206
0207 desc.add<std::string>("dqm_dir", "AlCaReco/PPSAlignment");
0208 desc.add<std::vector<std::string>>("sequence", {"x_alignment", "x_alignment_relative", "y_alignment"});
0209 desc.add<bool>("overwrite_sh_x", true);
0210 desc.add<std::string>("text_results_path", "./alignment_results.txt");
0211 desc.add<bool>("write_sqlite_results", false);
0212 desc.add<bool>("x_ali_rel_final_slope_fixed", false);
0213 desc.add<bool>("y_ali_final_slope_fixed", false);
0214 desc.add<double>("x_corr_min", -1'000'000.);
0215 desc.add<double>("x_corr_max", 1'000'000.);
0216 desc.add<double>("y_corr_min", -1'000'000.);
0217 desc.add<double>("y_corr_max", 1'000'000.);
0218 desc.add<unsigned int>("detector_id", 7);
0219 desc.add<unsigned int>("subdetector_id", 4);
0220 desc.add<bool>("debug", false);
0221
0222 descriptions.addWithDefaultLabel(desc);
0223 }
0224
0225 void PPSAlignmentHarvester::dqmEndJob(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter) {}
0226
0227 void PPSAlignmentHarvester::dqmEndRun(DQMStore::IBooker& iBooker,
0228 DQMStore::IGetter& iGetter,
0229 edm::Run const& iRun,
0230 edm::EventSetup const& iSetup) {
0231 const auto& cfg = iSetup.getData(esTokenTest_);
0232
0233 const auto& cfg_ref = iSetup.getData(esTokenReference_);
0234
0235
0236 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0237 for (const auto& rpc : {sc.rp_N_, sc.rp_F_}) {
0238 sh_x_map_[rpc.id_] = rpc.sh_x_;
0239 }
0240 }
0241 edm::LogInfo("PPSAlignmentHarvester").log([&](auto& li) {
0242 li << "Setting sh_x from config of:\n";
0243 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0244 for (const auto& rpc : {sc.rp_N_, sc.rp_F_}) {
0245 li << " " << rpc.name_ << " to " << std::fixed << std::setprecision(3) << rpc.sh_x_;
0246 if (rpc.name_ != "R_2_F")
0247 li << "\n";
0248 }
0249 }
0250 });
0251
0252 bool doXAli = false, doXAliRel = false, doYAli = false;
0253 for (const std::string& aliMethod : sequence_) {
0254 if (aliMethod == "x_alignment") {
0255 xAlignment(iBooker, iGetter, cfg, cfg_ref);
0256 doXAli = true;
0257 } else if (aliMethod == "x_alignment_relative") {
0258 xAlignmentRelative(iBooker, iGetter, cfg);
0259 doXAliRel = true;
0260 } else if (aliMethod == "y_alignment") {
0261 yAlignment(iBooker, iGetter, cfg);
0262 doYAli = true;
0263 } else
0264 edm::LogError("PPSAlignmentHarvester") << aliMethod << " is a wrong method name.";
0265 seqPos++;
0266 }
0267
0268
0269 CTPPSRPAlignmentCorrectionsData finalResults;
0270 if (doXAli) {
0271 finalResults.addCorrections(xAliResults_);
0272 if (doXAliRel) {
0273 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0274
0275 double d_x_N = xAliResults_.getRPCorrection(sc.rp_N_.id_).getShX();
0276 double d_x_F = xAliResults_.getRPCorrection(sc.rp_F_.id_).getShX();
0277
0278 double d_x_rel_N, d_x_rel_F;
0279 if (xAliRelFinalSlopeFixed_) {
0280 d_x_rel_N = xAliRelResultsSlopeFixed_.getRPCorrection(sc.rp_N_.id_).getShX();
0281 d_x_rel_F = xAliRelResultsSlopeFixed_.getRPCorrection(sc.rp_F_.id_).getShX();
0282 } else {
0283 d_x_rel_N = xAliRelResults_.getRPCorrection(sc.rp_N_.id_).getShX();
0284 d_x_rel_F = xAliRelResults_.getRPCorrection(sc.rp_F_.id_).getShX();
0285 }
0286
0287
0288 double b = d_x_rel_N - d_x_rel_F;
0289 double xCorrRel = b + d_x_F - d_x_N;
0290
0291 CTPPSRPAlignmentCorrectionData corrRelN(xCorrRel / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0292 finalResults.addRPCorrection(sc.rp_N_.id_, corrRelN);
0293 CTPPSRPAlignmentCorrectionData corrRelF(-xCorrRel / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0294 finalResults.addRPCorrection(sc.rp_F_.id_, corrRelF);
0295 }
0296 }
0297 }
0298 if (doYAli) {
0299 if (yAliFinalSlopeFixed_) {
0300 finalResults.addCorrections(yAliResultsSlopeFixed_);
0301 } else {
0302 finalResults.addCorrections(yAliResults_);
0303 }
0304 }
0305
0306
0307 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0308 for (const auto& rpc : {sc.rp_F_, sc.rp_N_}) {
0309 auto& rpResults = finalResults.getRPCorrection(rpc.id_);
0310
0311 if (!(xCorrRange_.first <= rpResults.getShX() && rpResults.getShX() <= xCorrRange_.second)) {
0312 edm::LogWarning("PPSAlignmentHarvester")
0313 << "The horizontal shift of " << rpc.name_ << " (" << std::fixed << std::setw(9) << std::setprecision(1)
0314 << rpResults.getShX() * 1000. << " um) outside of the reasonability range. Setting it to 0.";
0315 rpResults.setShX(0.);
0316 rpResults.setShXUnc(0.);
0317 }
0318
0319 if (!(yCorrRange_.first <= rpResults.getShY() && rpResults.getShY() <= yCorrRange_.second)) {
0320 edm::LogWarning("PPSAlignmentHarvester")
0321 << "The vertical shift of " << rpc.name_ << " (" << std::fixed << std::setw(9) << std::setprecision(1)
0322 << rpResults.getShY() * 1000. << " um) outside of the reasonability range. Setting it to 0.";
0323 rpResults.setShY(0.);
0324 rpResults.setShYUnc(0.);
0325 }
0326 }
0327 }
0328
0329
0330 edm::LogInfo("PPSAlignmentHarvester") << "final merged results:\n" << finalResults;
0331
0332 if (textResultsFile_.is_open()) {
0333 textResultsFile_ << "final merged results:\n" << finalResults;
0334 }
0335
0336
0337 if (writeSQLiteResults_) {
0338 CTPPSRPAlignmentCorrectionsData longIdFinalResults = getLongIdResults(finalResults);
0339 edm::LogInfo("PPSAlignmentHarvester") << "trying to store final merged results with long ids:\n"
0340 << longIdFinalResults;
0341
0342 edm::Service<cond::service::PoolDBOutputService> poolDbService;
0343 if (poolDbService.isAvailable()) {
0344 poolDbService->writeOneIOV(
0345 longIdFinalResults, poolDbService->currentTime(), "CTPPSRPAlignmentCorrectionsDataRcd");
0346 } else {
0347 edm::LogWarning("PPSAlignmentHarvester")
0348 << "Could not store the results in a DB object. PoolDBService not available.";
0349 }
0350 }
0351
0352
0353 if (debug_) {
0354 TDirectory* cutsDir = debugFile_->mkdir("cuts");
0355 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0356 TDirectory* sectorDir = cutsDir->mkdir(sc.name_.c_str());
0357
0358 gDirectory = sectorDir->mkdir("cut_h");
0359 auto* h2_cut_h_bef_monitor = iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/cuts/cut_h/h2_cut_h_bef");
0360 auto* h2_cut_h_aft_monitor = iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/cuts/cut_h/h2_cut_h_aft");
0361 writeCutPlot(
0362 h2_cut_h_bef_monitor->getTH2D(), sc.cut_h_a_, sc.cut_h_c_, cfg.n_si(), sc.cut_h_si_, "canvas_before");
0363 writeCutPlot(h2_cut_h_aft_monitor->getTH2D(), sc.cut_h_a_, sc.cut_h_c_, cfg.n_si(), sc.cut_h_si_, "canvas_after");
0364
0365 gDirectory = sectorDir->mkdir("cut_v");
0366 auto* h2_cut_v_bef_monitor = iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/cuts/cut_v/h2_cut_v_bef");
0367 auto* h2_cut_v_aft_monitor = iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/cuts/cut_v/h2_cut_v_aft");
0368 writeCutPlot(
0369 h2_cut_v_bef_monitor->getTH2D(), sc.cut_v_a_, sc.cut_v_c_, cfg.n_si(), sc.cut_v_si_, "canvas_before");
0370 writeCutPlot(h2_cut_v_aft_monitor->getTH2D(), sc.cut_v_a_, sc.cut_v_c_, cfg.n_si(), sc.cut_v_si_, "canvas_after");
0371 }
0372 }
0373 }
0374
0375
0376
0377
0378 std::unique_ptr<TGraphErrors> PPSAlignmentHarvester::buildGraphFromVector(
0379 const std::vector<PPSAlignmentConfiguration::PointErrors>& pv) {
0380 auto g = std::make_unique<TGraphErrors>();
0381
0382 for (unsigned int i = 0; i < pv.size(); i++) {
0383 const auto& p = pv[i];
0384 g->SetPoint(i, p.x_, p.y_);
0385 g->SetPointError(i, p.ex_, p.ey_);
0386 }
0387 g->Sort();
0388
0389 return g;
0390 }
0391
0392
0393 std::unique_ptr<TGraphErrors> PPSAlignmentHarvester::buildGraphFromMonitorElements(
0394 DQMStore::IGetter& iGetter,
0395 const PPSAlignmentConfiguration::RPConfig& rpc,
0396 const std::vector<MonitorElement*>& mes,
0397 const unsigned int fitProfileMinBinEntries,
0398 const unsigned int fitProfileMinNReasonable) {
0399 auto g = std::make_unique<TGraphErrors>();
0400
0401 for (auto* me : mes) {
0402 if (me->getName() == "h_y")
0403 {
0404
0405 std::string parentPath = me->getPathname();
0406 size_t parentPos = parentPath.substr(0, parentPath.size() - 1).find_last_of('/') + 1;
0407 std::string parentName = parentPath.substr(parentPos);
0408 std::replace(parentName.begin(), parentName.end(), '_', '.');
0409 size_t d = parentName.find('-');
0410 const double x_min = std::stod(parentName.substr(0, d));
0411 const double x_max = std::stod(parentName.substr(d + 1));
0412
0413 TH1D* h_y = me->getTH1D();
0414
0415
0416 auto* p_y_diffFN_vs_y_monitor = iGetter.get(parentPath + "p_y_diffFN_vs_y");
0417 if (p_y_diffFN_vs_y_monitor == nullptr) {
0418 edm::LogWarning("PPSAlignmentHarvester") << "[x_alignment] could not find p_y_diffFN_vs_y in: " << parentPath;
0419 continue;
0420 }
0421 TProfile* p_y_diffFN_vs_y = p_y_diffFN_vs_y_monitor->getTProfile();
0422
0423 double y_cen = h_y->GetMean() + rpc.y_cen_add_;
0424 double y_width = h_y->GetRMS() * rpc.y_width_mult_;
0425
0426 double sl, sl_unc;
0427 int fr = alig_utils::fitProfile(
0428 p_y_diffFN_vs_y, y_cen, y_width, fitProfileMinBinEntries, fitProfileMinNReasonable, sl, sl_unc);
0429 if (fr != 0)
0430 continue;
0431
0432 if (debug_)
0433 p_y_diffFN_vs_y->Write(parentName.c_str());
0434
0435 int idx = g->GetN();
0436 g->SetPoint(idx, (x_max + x_min) / 2., sl);
0437 g->SetPointError(idx, (x_max - x_min) / 2., sl_unc);
0438 }
0439 }
0440 g->Sort();
0441
0442 return g;
0443 }
0444
0445
0446 void PPSAlignmentHarvester::doMatch(DQMStore::IBooker& iBooker,
0447 const PPSAlignmentConfiguration& cfg,
0448 const PPSAlignmentConfiguration::RPConfig& rpc,
0449 TGraphErrors* g_ref,
0450 TGraphErrors* g_test,
0451 const PPSAlignmentConfiguration::SelectionRange& range_ref,
0452 const double sh_min,
0453 const double sh_max,
0454 double& sh_best,
0455 double& sh_best_unc) {
0456 const auto range_test = cfg.alignment_x_meth_o_ranges().at(rpc.id_);
0457
0458
0459 edm::LogInfo("PPSAlignmentHarvester") << std::fixed << std::setprecision(3) << "[x_alignment] "
0460 << "ref: x_min = " << range_ref.x_min_ << ", x_max = " << range_ref.x_max_
0461 << "\n"
0462 << "test: x_min = " << range_test.x_min_ << ", x_max = " << range_test.x_max_;
0463
0464
0465 auto s_ref = std::make_unique<TSpline3>("s_ref", g_ref->GetX(), g_ref->GetY(), g_ref->GetN());
0466
0467
0468 auto g_n_points = std::make_unique<TGraph>();
0469 g_n_points->SetName("g_n_points");
0470 g_n_points->SetTitle(";sh;N");
0471 auto g_chi_sq = std::make_unique<TGraph>();
0472 g_chi_sq->SetName("g_chi_sq");
0473 g_chi_sq->SetTitle(";sh;S2");
0474 auto g_chi_sq_norm = std::make_unique<TGraph>();
0475 g_chi_sq_norm->SetName("g_chi_sq_norm");
0476 g_chi_sq_norm->SetTitle(";sh;S2 / N");
0477
0478
0479 double S2_norm_best = 1E100;
0480
0481 for (double sh = sh_min; sh <= sh_max; sh += cfg.x_ali_sh_step()) {
0482
0483 int n_points = 0;
0484 double S2 = 0.;
0485
0486 for (int i = 0; i < g_test->GetN(); ++i) {
0487 const double x_test = g_test->GetX()[i];
0488 const double y_test = g_test->GetY()[i];
0489 const double y_test_unc = g_test->GetErrorY(i);
0490
0491 const double x_ref = x_test + sh;
0492
0493 if (x_ref < range_ref.x_min_ || x_ref > range_ref.x_max_ || x_test < range_test.x_min_ ||
0494 x_test > range_test.x_max_)
0495 continue;
0496
0497 const double y_ref = s_ref->Eval(x_ref);
0498
0499 int js = -1, jg = -1;
0500 double xs = -1E100, xg = +1E100;
0501 for (int j = 0; j < g_ref->GetN(); ++j) {
0502 const double x = g_ref->GetX()[j];
0503 if (x < x_ref && x > xs) {
0504 xs = x;
0505 js = j;
0506 }
0507 if (x > x_ref && x < xg) {
0508 xg = x;
0509 jg = j;
0510 }
0511 }
0512 if (jg == -1)
0513 jg = js;
0514
0515 const double y_ref_unc = (g_ref->GetErrorY(js) + g_ref->GetErrorY(jg)) / 2.;
0516
0517 n_points++;
0518 const double S2_inc = pow(y_test - y_ref, 2.) / (y_ref_unc * y_ref_unc + y_test_unc * y_test_unc);
0519 S2 += S2_inc;
0520 }
0521
0522
0523 double S2_norm = S2 / n_points;
0524
0525 if (S2_norm < S2_norm_best) {
0526 S2_norm_best = S2_norm;
0527 sh_best = sh;
0528 }
0529
0530
0531 int idx = g_n_points->GetN();
0532 g_n_points->SetPoint(idx, sh, n_points);
0533 g_chi_sq->SetPoint(idx, sh, S2);
0534 g_chi_sq_norm->SetPoint(idx, sh, S2_norm);
0535 }
0536
0537 auto ff_pol2 = std::make_unique<TF1>("ff_pol2", "[0] + [1]*x + [2]*x*x");
0538
0539
0540 double fit_range = cfg.methOUncFitRange();
0541 g_chi_sq->Fit(ff_pol2.get(), "Q", "", sh_best - fit_range, sh_best + fit_range);
0542 sh_best_unc = 1. / sqrt(ff_pol2->GetParameter(2));
0543
0544
0545 edm::LogInfo("PPSAlignmentHarvester") << std::fixed << std::setprecision(3) << "[x_alignment] "
0546 << "sh_best = (" << sh_best << " +- " << sh_best_unc << ") mm";
0547
0548 auto g_test_shifted = std::make_unique<TGraphErrors>(*g_test);
0549 for (int i = 0; i < g_test_shifted->GetN(); ++i) {
0550 g_test_shifted->GetX()[i] += sh_best;
0551 }
0552
0553 std::unique_ptr<TH1D> histPtr = getTH1DFromTGraphErrors(
0554 g_test_shifted.get(), "test_shifted", ";x (mm);S", rpc.x_slice_n_, rpc.x_slice_w_, rpc.x_slice_min_ + sh_best);
0555 iBooker.book1DD("h_test_shifted", histPtr.get());
0556
0557 if (debug_) {
0558
0559 g_n_points->Write();
0560 g_chi_sq->Write();
0561 g_chi_sq_norm->Write();
0562 g_test_shifted->SetTitle(";x (mm);S");
0563 g_test_shifted->Write("g_test_shifted");
0564
0565
0566 auto g_results = std::make_unique<TGraph>();
0567 g_results->SetName("g_results");
0568 g_results->SetPoint(0, sh_best, sh_best_unc);
0569 g_results->SetPoint(1, range_ref.x_min_, range_ref.x_max_);
0570 g_results->SetPoint(2, range_test.x_min_, range_test.x_max_);
0571 g_results->Write();
0572
0573
0574 auto c_cmp = std::make_unique<TCanvas>("c_cmp");
0575 g_ref->SetLineColor(1);
0576 g_ref->SetName("g_ref");
0577 g_ref->Draw("apl");
0578
0579 g_test->SetLineColor(4);
0580 g_test->SetName("g_test");
0581 g_test->Draw("pl");
0582
0583 g_test_shifted->SetLineColor(2);
0584 g_test_shifted->SetName("g_test_shifted");
0585
0586 g_test_shifted->Draw("pl");
0587 c_cmp->Write();
0588 }
0589 }
0590
0591
0592 void PPSAlignmentHarvester::xAlignment(DQMStore::IBooker& iBooker,
0593 DQMStore::IGetter& iGetter,
0594 const PPSAlignmentConfiguration& cfg,
0595 const PPSAlignmentConfiguration& cfg_ref) {
0596 TDirectory* xAliDir = nullptr;
0597 if (debug_)
0598 xAliDir = debugFile_->mkdir((std::to_string(seqPos) + ": x alignment").c_str());
0599
0600 for (const auto& [sc, sc_ref] : {std::make_pair(cfg.sectorConfig45(), cfg_ref.sectorConfig45()),
0601 std::make_pair(cfg.sectorConfig56(), cfg_ref.sectorConfig56())}) {
0602 for (const auto& [rpc, rpc_ref] :
0603 {std::make_pair(sc.rp_F_, sc_ref.rp_F_), std::make_pair(sc.rp_N_, sc_ref.rp_N_)}) {
0604 auto mes_test = iGetter.getAllContents(dqmDir_ + "/worker/" + sc.name_ + "/near_far/x slices " + rpc.position_);
0605 if (mes_test.empty()) {
0606 edm::LogWarning("PPSAlignmentHarvester") << "[x_alignment] " << rpc.name_ << ": could not load mes_test";
0607 continue;
0608 }
0609
0610 TDirectory* rpDir = nullptr;
0611 if (debug_)
0612 rpDir = xAliDir->mkdir(rpc.name_.c_str());
0613
0614 auto vec_ref = cfg_ref.matchingReferencePoints().at(rpc.id_);
0615 if (vec_ref.empty()) {
0616 edm::LogInfo("PPSAlignmentHarvester") << "[x_alignment] " << rpc.name_ << ": reference points vector is empty";
0617 continue;
0618 }
0619
0620 std::unique_ptr<TGraphErrors> g_ref = buildGraphFromVector(vec_ref);
0621
0622 if (debug_)
0623 gDirectory = rpDir->mkdir("fits_test");
0624 std::unique_ptr<TGraphErrors> g_test = buildGraphFromMonitorElements(
0625 iGetter, rpc, mes_test, cfg.fitProfileMinBinEntries(), cfg.fitProfileMinNReasonable());
0626
0627
0628 if (g_ref->GetN() < (int)cfg.methOGraphMinN() || g_test->GetN() < (int)cfg.methOGraphMinN()) {
0629 edm::LogWarning("PPSAlignmentHarvester")
0630 << "[x_alignment] " << rpc.name_ << ": insufficient data, skipping (g_ref " << g_ref->GetN() << "/"
0631 << cfg.methOGraphMinN() << ", g_test " << g_test->GetN() << "/" << cfg.methOGraphMinN() << ")";
0632 continue;
0633 }
0634
0635 iBooker.setCurrentFolder(dqmDir_ + "/harvester/x alignment/" + rpc.name_);
0636
0637 std::unique_ptr<TH1D> histPtr = getTH1DFromTGraphErrors(
0638 g_ref.get(), "ref", ";x (mm);S", rpc_ref.x_slice_n_, rpc_ref.x_slice_w_, rpc_ref.x_slice_min_);
0639 iBooker.book1DD("h_ref", histPtr.get());
0640
0641 histPtr =
0642 getTH1DFromTGraphErrors(g_test.get(), "test", ";x (mm);S", rpc.x_slice_n_, rpc.x_slice_w_, rpc.x_slice_min_);
0643 iBooker.book1DD("h_test", histPtr.get());
0644
0645 if (debug_) {
0646 gDirectory = rpDir;
0647 g_ref->SetTitle(";x (mm);S");
0648 g_ref->Write("g_ref");
0649 g_test->SetTitle(";x (mm);S");
0650 g_test->Write("g_test");
0651 }
0652
0653 const auto& shiftRange = cfg.matchingShiftRanges().at(rpc.id_);
0654 double sh = 0., sh_unc = 0.;
0655
0656
0657 doMatch(iBooker,
0658 cfg,
0659 rpc,
0660 g_ref.get(),
0661 g_test.get(),
0662 cfg_ref.alignment_x_meth_o_ranges().at(rpc.id_),
0663 shiftRange.x_min_,
0664 shiftRange.x_max_,
0665 sh,
0666 sh_unc);
0667
0668
0669 CTPPSRPAlignmentCorrectionData rpResult(sh, sh_unc, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0670 xAliResults_.setRPCorrection(rpc.id_, rpResult);
0671 edm::LogInfo("PPSAlignmentHarvester") << std::fixed << std::setprecision(3) << "[x_alignment] "
0672 << "Setting sh_x of " << rpc.name_ << " to " << sh;
0673
0674
0675 if (overwriteShX_) {
0676 sh_x_map_[rpc.id_] = sh;
0677 }
0678 }
0679 }
0680
0681 edm::LogInfo("PPSAlignmentHarvester") << seqPos << ": x_alignment:\n" << xAliResults_;
0682
0683 if (textResultsFile_.is_open())
0684 textResultsFile_ << seqPos << ": x_alignment:\n" << xAliResults_ << "\n\n";
0685 }
0686
0687
0688
0689 void PPSAlignmentHarvester::xAlignmentRelative(DQMStore::IBooker& iBooker,
0690 DQMStore::IGetter& iGetter,
0691 const PPSAlignmentConfiguration& cfg) {
0692 TDirectory* xAliRelDir = nullptr;
0693 if (debug_)
0694 xAliRelDir = debugFile_->mkdir((std::to_string(seqPos) + ": x_alignment_relative").c_str());
0695
0696 auto ff = std::make_unique<TF1>("ff", "[0] + [1]*(x - [2])");
0697 auto ff_sl_fix = std::make_unique<TF1>("ff_sl_fix", "[0] + [1]*(x - [2])");
0698
0699
0700 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0701 TDirectory* sectorDir = nullptr;
0702 if (debug_) {
0703 sectorDir = xAliRelDir->mkdir(sc.name_.c_str());
0704 gDirectory = sectorDir;
0705 }
0706
0707 auto* p_x_diffFN_vs_x_N_monitor = iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/near_far/p_x_diffFN_vs_x_N");
0708 if (p_x_diffFN_vs_x_N_monitor == nullptr) {
0709 edm::LogWarning("PPSAlignmentHarvester")
0710 << "[x_alignment_relative] " << sc.name_ << ": cannot load data, skipping";
0711 continue;
0712 }
0713 TProfile* p_x_diffFN_vs_x_N = p_x_diffFN_vs_x_N_monitor->getTProfile();
0714
0715 if (p_x_diffFN_vs_x_N->GetEntries() < cfg.nearFarMinEntries()) {
0716 edm::LogWarning("PPSAlignmentHarvester")
0717 << "[x_alignment_relative] " << sc.name_ << ": insufficient data, skipping (near_far "
0718 << p_x_diffFN_vs_x_N->GetEntries() << "/" << cfg.nearFarMinEntries() << ")";
0719 continue;
0720 }
0721
0722 const double x_min = cfg.alignment_x_relative_ranges().at(sc.rp_N_.id_).x_min_;
0723 const double x_max = cfg.alignment_x_relative_ranges().at(sc.rp_N_.id_).x_max_;
0724
0725 const double sh_x_N = sh_x_map_[sc.rp_N_.id_];
0726 double slope = sc.slope_;
0727
0728
0729 ff->SetParameters(0., slope, 0.);
0730 ff->FixParameter(2, -sh_x_N);
0731 ff->SetLineColor(2);
0732 p_x_diffFN_vs_x_N->Fit(ff.get(), "Q", "", x_min, x_max);
0733
0734 const double a = ff->GetParameter(1), a_unc = ff->GetParError(1);
0735 const double b = ff->GetParameter(0), b_unc = ff->GetParError(0);
0736
0737 edm::LogInfo("PPSAlignmentHarvester")
0738 << "[x_alignment_relative] " << sc.name_ << ":\n"
0739 << std::fixed << std::setprecision(3) << " x_min = " << x_min << ", x_max = " << x_max << "\n"
0740 << " sh_x_N = " << sh_x_N << ", slope (fix) = " << slope << ", slope (fitted) = " << a;
0741
0742
0743 CTPPSRPAlignmentCorrectionData rpResult_N(+b / 2., b_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0744 xAliRelResults_.setRPCorrection(sc.rp_N_.id_, rpResult_N);
0745 CTPPSRPAlignmentCorrectionData rpResult_F(-b / 2., b_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0746 xAliRelResults_.setRPCorrection(sc.rp_F_.id_, rpResult_F);
0747
0748
0749 ff_sl_fix->SetParameters(0., slope, 0.);
0750 ff_sl_fix->FixParameter(1, slope);
0751 ff_sl_fix->FixParameter(2, -sh_x_N);
0752 ff_sl_fix->SetLineColor(4);
0753 p_x_diffFN_vs_x_N->Fit(ff_sl_fix.get(), "Q+", "", x_min, x_max);
0754
0755 const double b_fs = ff_sl_fix->GetParameter(0), b_fs_unc = ff_sl_fix->GetParError(0);
0756
0757
0758 CTPPSRPAlignmentCorrectionData rpResult_sl_fix_N(+b_fs / 2., b_fs_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0759 xAliRelResultsSlopeFixed_.setRPCorrection(sc.rp_N_.id_, rpResult_sl_fix_N);
0760 CTPPSRPAlignmentCorrectionData rpResult_sl_fix_F(-b_fs / 2., b_fs_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0761 xAliRelResultsSlopeFixed_.setRPCorrection(sc.rp_F_.id_, rpResult_sl_fix_F);
0762
0763 edm::LogInfo("PPSAlignmentHarvester")
0764 << "[x_alignment_relative] " << std::fixed << std::setprecision(3) << "ff: " << ff->GetParameter(0) << " + "
0765 << ff->GetParameter(1) << " * (x - " << ff->GetParameter(2) << "), ff_sl_fix: " << ff_sl_fix->GetParameter(0)
0766 << " + " << ff_sl_fix->GetParameter(1) << " * (x - " << ff_sl_fix->GetParameter(2) << ")";
0767
0768
0769 iBooker.setCurrentFolder(dqmDir_ + "/harvester/x_alignment_relative/" + sc.name_);
0770 iBooker.bookProfile("p_x_diffFN_vs_x_N", p_x_diffFN_vs_x_N);
0771
0772 if (debug_) {
0773 p_x_diffFN_vs_x_N->Write("p_x_diffFN_vs_x_N");
0774
0775 auto g_results = std::make_unique<TGraph>();
0776 g_results->SetPoint(0, sh_x_N, 0.);
0777 g_results->SetPoint(1, a, a_unc);
0778 g_results->SetPoint(2, b, b_unc);
0779 g_results->SetPoint(3, b_fs, b_fs_unc);
0780 g_results->Write("g_results");
0781 }
0782 }
0783
0784
0785 edm::LogInfo("PPSAlignmentHarvester") << seqPos << ": x_alignment_relative:\n"
0786 << xAliRelResults_ << seqPos + 1 << ": x_alignment_relative_sl_fix:\n"
0787 << xAliRelResultsSlopeFixed_;
0788
0789 if (textResultsFile_.is_open()) {
0790 textResultsFile_ << seqPos << ": x_alignment_relative:\n" << xAliRelResults_ << "\n";
0791 textResultsFile_ << seqPos << ": x_alignment_relative_sl_fix:\n" << xAliRelResultsSlopeFixed_ << "\n\n";
0792 }
0793 }
0794
0795
0796
0797 double PPSAlignmentHarvester::findMax(const TF1* ff_fit) {
0798 const double mu = ff_fit->GetParameter(1);
0799 const double si = ff_fit->GetParameter(2);
0800
0801
0802 if (si > 25. || std::fabs(mu) > 100.)
0803 return 1E100;
0804
0805 double x_max = 1E100;
0806 double y_max = -1E100;
0807 for (double x = mu - si; x <= mu + si; x += 0.001) {
0808 double y = ff_fit->Eval(x);
0809 if (y > y_max) {
0810 x_max = x;
0811 y_max = y;
0812 }
0813 }
0814
0815 return x_max;
0816 }
0817
0818 TH1D* PPSAlignmentHarvester::buildModeGraph(DQMStore::IBooker& iBooker,
0819 const MonitorElement* h2_y_vs_x,
0820 const PPSAlignmentConfiguration& cfg,
0821 const PPSAlignmentConfiguration::RPConfig& rpc) {
0822 TDirectory* d_top = nullptr;
0823 if (debug_)
0824 d_top = gDirectory;
0825
0826 auto ff_fit = std::make_unique<TF1>("ff_fit", "[0] * exp(-(x-[1])*(x-[1])/2./[2]/[2]) + [3] + [4]*x");
0827
0828 int h_n = h2_y_vs_x->getNbinsX();
0829 double diff = h2_y_vs_x->getTH2D()->GetXaxis()->GetBinWidth(1) / 2.;
0830 auto h_mode = iBooker.book1DD("mode",
0831 ";x (mm); mode of y (mm)",
0832 h_n,
0833 h2_y_vs_x->getTH2D()->GetXaxis()->GetBinCenter(1) - diff,
0834 h2_y_vs_x->getTH2D()->GetXaxis()->GetBinCenter(h_n) + diff);
0835
0836
0837 for (int bix = 1; bix <= h_n; bix++) {
0838 const double x = h2_y_vs_x->getTH2D()->GetXaxis()->GetBinCenter(bix);
0839
0840 char buf[100];
0841 sprintf(buf, "h_y_x=%.3f", x);
0842 TH1D* h_y = h2_y_vs_x->getTH2D()->ProjectionY(buf, bix, bix);
0843
0844 if (h_y->GetEntries() < cfg.multSelProjYMinEntries())
0845 continue;
0846
0847 if (debug_) {
0848 sprintf(buf, "x=%.3f", x);
0849 gDirectory = d_top->mkdir(buf);
0850 }
0851
0852 double conMax = -1.;
0853 double conMax_x = 0.;
0854 for (int biy = 1; biy < h_y->GetNbinsX(); biy++) {
0855 if (h_y->GetBinContent(biy) > conMax) {
0856 conMax = h_y->GetBinContent(biy);
0857 conMax_x = h_y->GetBinCenter(biy);
0858 }
0859 }
0860
0861 ff_fit->SetParameters(conMax, conMax_x, h_y->GetRMS() * 0.75, 0., 0.);
0862 ff_fit->FixParameter(4, 0.);
0863
0864 double x_min = rpc.x_min_fit_mode_, x_max = rpc.x_max_fit_mode_;
0865 h_y->Fit(ff_fit.get(), "Q", "", x_min, x_max);
0866
0867 ff_fit->ReleaseParameter(4);
0868 double w = std::min(4., 2. * ff_fit->GetParameter(2));
0869 x_min = ff_fit->GetParameter(1) - w;
0870 x_max = std::min(rpc.y_max_fit_mode_, ff_fit->GetParameter(1) + w);
0871
0872 h_y->Fit(ff_fit.get(), "Q", "", x_min, x_max);
0873
0874 if (debug_)
0875 h_y->Write("h_y");
0876
0877 double y_mode = findMax(ff_fit.get());
0878 const double y_mode_fit_unc = ff_fit->GetParameter(2) / 10;
0879 const double y_mode_sys_unc = cfg.y_mode_sys_unc();
0880 double y_mode_unc = std::sqrt(y_mode_fit_unc * y_mode_fit_unc + y_mode_sys_unc * y_mode_sys_unc);
0881
0882 const double chiSqThreshold = cfg.chiSqThreshold();
0883
0884 const bool valid =
0885 !(std::fabs(y_mode_unc) > cfg.y_mode_unc_max_valid() || std::fabs(y_mode) > cfg.y_mode_max_valid() ||
0886 ff_fit->GetChisquare() / ff_fit->GetNDF() > chiSqThreshold);
0887
0888 if (debug_) {
0889 auto g_data = std::make_unique<TGraph>();
0890 g_data->SetPoint(0, y_mode, y_mode_unc);
0891 g_data->SetPoint(1, ff_fit->GetChisquare(), ff_fit->GetNDF());
0892 g_data->SetPoint(2, valid, 0.);
0893 g_data->Write("g_data");
0894 }
0895
0896 if (!valid)
0897 continue;
0898
0899 h_mode->Fill(x, y_mode);
0900 h_mode->setBinError(bix, y_mode_unc);
0901 }
0902
0903 return h_mode->getTH1D();
0904 }
0905
0906 void PPSAlignmentHarvester::yAlignment(DQMStore::IBooker& iBooker,
0907 DQMStore::IGetter& iGetter,
0908 const PPSAlignmentConfiguration& cfg) {
0909 TDirectory* yAliDir = nullptr;
0910 if (debug_)
0911 yAliDir = debugFile_->mkdir((std::to_string(seqPos) + ": y_alignment").c_str());
0912
0913 auto ff = std::make_unique<TF1>("ff", "[0] + [1]*(x - [2])");
0914 auto ff_sl_fix = std::make_unique<TF1>("ff_sl_fix", "[0] + [1]*(x - [2])");
0915
0916
0917 for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0918 for (const auto& rpc : {sc.rp_F_, sc.rp_N_}) {
0919 TDirectory* rpDir = nullptr;
0920 if (debug_) {
0921 rpDir = yAliDir->mkdir(rpc.name_.c_str());
0922 gDirectory = rpDir->mkdir("x");
0923 }
0924
0925 auto* h2_y_vs_x =
0926 iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/multiplicity selection/" + rpc.name_ + "/h2_y_vs_x");
0927
0928 if (h2_y_vs_x == nullptr) {
0929 edm::LogWarning("PPSAlignmentHarvester") << "[y_alignment] " << rpc.name_ << ": cannot load data, skipping";
0930 continue;
0931 }
0932
0933 iBooker.setCurrentFolder(dqmDir_ + "/harvester/y alignment/" + rpc.name_);
0934 auto* h_y_cen_vs_x = buildModeGraph(iBooker, h2_y_vs_x, cfg, rpc);
0935
0936 if ((unsigned int)h_y_cen_vs_x->GetEntries() < cfg.modeGraphMinN()) {
0937 edm::LogWarning("PPSAlignmentHarvester")
0938 << "[y_alignment] " << rpc.name_ << ": insufficient data, skipping (mode graph "
0939 << h_y_cen_vs_x->GetEntries() << "/" << cfg.modeGraphMinN() << ")";
0940 continue;
0941 }
0942
0943 const double x_min = cfg.alignment_y_ranges().at(rpc.id_).x_min_;
0944 const double x_max = cfg.alignment_y_ranges().at(rpc.id_).x_max_;
0945
0946 const double sh_x = sh_x_map_[rpc.id_];
0947 double slope = rpc.slope_;
0948
0949
0950 ff->SetParameters(0., 0., 0.);
0951 ff->FixParameter(2, -sh_x);
0952 ff->SetLineColor(2);
0953 h_y_cen_vs_x->Fit(ff.get(), "Q", "", x_min, x_max);
0954
0955 const double a = ff->GetParameter(1), a_unc = ff->GetParError(1);
0956 const double b = ff->GetParameter(0), b_unc = ff->GetParError(0);
0957
0958 edm::LogInfo("PPSAlignmentHarvester")
0959 << "[y_alignment] " << rpc.name_ << ":\n"
0960 << std::fixed << std::setprecision(3) << " x_min = " << x_min << ", x_max = " << x_max << "\n"
0961 << " sh_x = " << sh_x << ", slope (fix) = " << slope << ", slope (fitted) = " << a;
0962
0963
0964 CTPPSRPAlignmentCorrectionData rpResult(0., 0., -b, b_unc, 0., 0., 0., 0., 0., 0., 0., 0.);
0965 yAliResults_.setRPCorrection(rpc.id_, rpResult);
0966
0967
0968 ff_sl_fix->SetParameters(0., 0., 0.);
0969 ff_sl_fix->FixParameter(1, slope);
0970 ff_sl_fix->FixParameter(2, -sh_x);
0971 ff_sl_fix->SetLineColor(4);
0972 h_y_cen_vs_x->Fit(ff_sl_fix.get(), "Q+", "", x_min, x_max);
0973
0974 const double b_fs = ff_sl_fix->GetParameter(0), b_fs_unc = ff_sl_fix->GetParError(0);
0975
0976 CTPPSRPAlignmentCorrectionData rpResult_sl_fix(0., 0., -b_fs, b_fs_unc, 0., 0., 0., 0., 0., 0., 0., 0.);
0977 yAliResultsSlopeFixed_.setRPCorrection(rpc.id_, rpResult_sl_fix);
0978
0979 edm::LogInfo("PPSAlignmentHarvester")
0980 << "[y_alignment] " << std::fixed << std::setprecision(3) << "ff: " << ff->GetParameter(0) << " + "
0981 << ff->GetParameter(1) << " * (x - " << ff->GetParameter(2) << "), ff_sl_fix: " << ff_sl_fix->GetParameter(0)
0982 << " + " << ff_sl_fix->GetParameter(1) << " * (x - " << ff_sl_fix->GetParameter(2) << ")";
0983
0984 if (debug_) {
0985 gDirectory = rpDir;
0986
0987 h_y_cen_vs_x->SetTitle(";x (mm); mode of y (mm)");
0988 h_y_cen_vs_x->Write("h_y_cen_vs_x");
0989
0990 auto g_results = std::make_unique<TGraph>();
0991 g_results->SetPoint(0, sh_x, 0.);
0992 g_results->SetPoint(1, a, a_unc);
0993 g_results->SetPoint(2, -b, b_unc);
0994 g_results->SetPoint(3, -b_fs, b_fs_unc);
0995 g_results->Write("g_results");
0996 }
0997 }
0998 }
0999
1000
1001 edm::LogInfo("PPSAlignmentHarvester") << seqPos << ": y_alignment:\n"
1002 << yAliResults_ << seqPos << ": y_alignment_sl_fix:\n"
1003 << yAliResultsSlopeFixed_;
1004
1005 if (textResultsFile_.is_open()) {
1006 textResultsFile_ << seqPos << ": y_alignment:\n" << yAliResults_ << "\n";
1007 textResultsFile_ << seqPos << ": y_alignment_sl_fix:\n" << yAliResultsSlopeFixed_ << "\n\n";
1008 }
1009 }
1010
1011
1012
1013
1014 void PPSAlignmentHarvester::writeCutPlot(
1015 TH2D* h, const double a, const double c, const double n_si, const double si, const std::string& label) {
1016 auto canvas = std::make_unique<TCanvas>();
1017 canvas->SetName(label.c_str());
1018 canvas->SetLogz(1);
1019
1020 h->Draw("colz");
1021
1022 double x_min = -30.;
1023 double x_max = 30.;
1024
1025 auto g_up = std::make_unique<TGraph>();
1026 g_up->SetName("g_up");
1027 g_up->SetPoint(0, x_min, -a * x_min - c + n_si * si);
1028 g_up->SetPoint(1, x_max, -a * x_max - c + n_si * si);
1029 g_up->SetLineColor(1);
1030 g_up->Draw("l");
1031
1032 auto g_down = std::make_unique<TGraph>();
1033 g_down->SetName("g_down");
1034 g_down->SetPoint(0, x_min, -a * x_min - c - n_si * si);
1035 g_down->SetPoint(1, x_max, -a * x_max - c - n_si * si);
1036 g_down->SetLineColor(1);
1037 g_down->Draw("l");
1038
1039 canvas->Write();
1040 }
1041
1042
1043
1044 std::unique_ptr<TH1D> PPSAlignmentHarvester::getTH1DFromTGraphErrors(
1045 TGraphErrors* graph, const std::string& title, const std::string& labels, int n, double binWidth, double min) {
1046 std::unique_ptr<TH1D> hist;
1047 if (n == 0) {
1048 hist = std::make_unique<TH1D>(title.c_str(), labels.c_str(), 0, -10., 10.);
1049 } else if (n == 1) {
1050 hist = std::make_unique<TH1D>(title.c_str(), labels.c_str(), 1, graph->GetPointX(0) - 5., graph->GetPointX(0) + 5.);
1051 } else {
1052 n = n == -1 ? graph->GetN() : n;
1053 binWidth = binWidth == -1 ? graph->GetPointX(1) - graph->GetPointX(0) : binWidth;
1054 double diff = binWidth / 2.;
1055 min = min == -1 ? graph->GetPointX(0) - diff : min;
1056 double max = min + n * binWidth;
1057 hist = std::make_unique<TH1D>(title.c_str(), labels.c_str(), n, min, max);
1058 }
1059
1060 for (int i = 0; i < graph->GetN(); i++) {
1061 double x, y;
1062 graph->GetPoint(i, x, y);
1063 hist->Fill(x, y);
1064 hist->SetBinError(hist->GetXaxis()->FindBin(x), graph->GetErrorY(i));
1065 }
1066 return hist;
1067 }
1068
1069
1070 CTPPSRPAlignmentCorrectionsData PPSAlignmentHarvester::getLongIdResults(CTPPSRPAlignmentCorrectionsData shortIdResults) {
1071 CTPPSRPAlignmentCorrectionsData longIdResults;
1072 for (const auto& [shortId, correction] : shortIdResults.getRPMap()) {
1073 unsigned int arm = shortId / 100;
1074 unsigned int station = (shortId / 10) % 10;
1075 unsigned int rp = shortId % 10;
1076
1077 uint32_t longDetId = detectorId_ << 28 | subdetectorId_ << 25 | arm << 24 | station << 22 | rp << 19;
1078
1079 longIdResults.addRPCorrection(longDetId, correction);
1080 }
1081
1082 return longIdResults;
1083 }
1084
1085 DEFINE_FWK_MODULE(PPSAlignmentHarvester);