Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-12-08 23:57:35 

0001 /****************************************************************************

0002 * Authors: 

0003 *  Jan Kašpar (jan.kaspar@gmail.com) 

0004 *  Mateusz Kocot (mateuszkocot99@gmail.com)

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   // ------------ x alignment ------------

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   // ------------ x alignment relative ------------

0090   void xAlignmentRelative(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter, const PPSAlignmentConfiguration& cfg);
0091 
0092   // ------------ y alignment ------------

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   // ------------ other member data and methods ------------

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   // variables from parameters

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   // other class variables

0126   std::unique_ptr<TFile> debugFile_;
0127   std::ofstream textResultsFile_;
0128   int seqPos = 1;  // position in sequence_

0129 
0130   CTPPSRPAlignmentCorrectionsData xAliResults_;
0131 
0132   CTPPSRPAlignmentCorrectionsData xAliRelResults_;
0133   CTPPSRPAlignmentCorrectionsData xAliRelResultsSlopeFixed_;
0134 
0135   CTPPSRPAlignmentCorrectionsData yAliResults_;
0136   CTPPSRPAlignmentCorrectionsData yAliResultsSlopeFixed_;
0137 };
0138 
0139 // -------------------------------- DQMEDHarvester methods --------------------------------

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.)),  // um -> mm

0154       yCorrRange_(std::make_pair(iConfig.getParameter<double>("y_corr_min") / 1000.,
0155                                  iConfig.getParameter<double>("y_corr_max") / 1000.)),  // um -> mm

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     // print in um

0178     li << "* x_corr_min: " << std::fixed << xCorrRange_.first * 1000. << ", x_corr_max: " << xCorrRange_.second * 1000.
0179        << "\n";
0180     // print in um

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   // setting default sh_x values from config

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   // merge results from all the specified methods

0255   CTPPSRPAlignmentCorrectionsData finalResults;
0256   if (doXAli) {  // x alignment

0257     finalResults.addCorrections(xAliResults_);
0258     if (doXAliRel) {  // merge with x alignment relative

0259       for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0260         // extract shifts

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         // merge the results

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) {  // y alignment

0285     if (yAliFinalSlopeFixed_) {
0286       finalResults.addCorrections(yAliResultsSlopeFixed_);
0287     } else {
0288       finalResults.addCorrections(yAliResults_);
0289     }
0290   }
0291 
0292   // check if the results are within the reasonability ranges xCorrRange and yCorrRange

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   // print the text results

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   // if requested, store the results in a DB object

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   // if debug_, save nice-looking cut plots with the worker data in the debug ROOT file

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 // -------------------------------- x alignment methods --------------------------------

0357 
0358 // Builds graph from a vector of points (with errors).

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 // Builds a TGraphErrors from slice plots represented as MonitorElements.

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")  // find "h_y"

0384     {
0385       // retrieve parent directory

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       // collect "p_y_diffFN_vs_y" corresponding to found "h_y"

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 // Matches reference data with test data.

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   // print config

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   // make spline from g_ref

0445   auto s_ref = std::make_unique<TSpline3>("s_ref", g_ref->GetX(), g_ref->GetY(), g_ref->GetN());
0446 
0447   // book match-quality graphs

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   // optimalisation variables

0459   double S2_norm_best = 1E100;
0460 
0461   for (double sh = sh_min; sh <= sh_max; sh += cfg.x_ali_sh_step()) {
0462     // calculate chi^2

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     // update best result

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     // fill in graphs

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   // determine uncertainty

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   // print results

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     // save graphs

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     // save results

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     // save debug canvas

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 // method o

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       // require minimal number of points

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       // matching

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       // save the results

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       // update the shift

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 // -------------------------------- x alignment relative methods --------------------------------

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   // processing

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     // calculate the results without slope fixed

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     // calculate the results with slope fixed

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     // rebook the diffFN plot in the harvester

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   // write results

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 // -------------------------------- y alignment methods --------------------------------

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   // unreasonable fit?

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   // find mode for each bin

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   // processing

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       // calculate the results without slope fixed

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       // calculate the results with slope fixed

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   // write results

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 // -------------------------------- other methods --------------------------------

0989 
0990 // Creates a plot showing a cut applied by the worker. Used only for debug purposes.

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 // Points in TGraph should be sorted (TGraph::Sort())

1020 // if n, binWidth, or min is set to -1, method will find it on its own

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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!