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