Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-09-11 22:29:07 

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

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

0093   void xAlignmentRelative(DQMStore::IBooker& iBooker, DQMStore::IGetter& iGetter, const PPSAlignmentConfiguration& cfg);
0094 
0095   // ------------ y alignment ------------

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

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

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

0133   std::unique_ptr<TFile> debugFile_;
0134   std::ofstream textResultsFile_;
0135   int seqPos = 1;  // position in sequence_

0136 
0137   CTPPSRPAlignmentCorrectionsData xAliResults_;
0138 
0139   CTPPSRPAlignmentCorrectionsData xAliRelResults_;
0140   CTPPSRPAlignmentCorrectionsData xAliRelResultsSlopeFixed_;
0141 
0142   CTPPSRPAlignmentCorrectionsData yAliResults_;
0143   CTPPSRPAlignmentCorrectionsData yAliResultsSlopeFixed_;
0144 };
0145 
0146 // -------------------------------- DQMEDHarvester methods --------------------------------

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

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

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

0187     li << "* x_corr_min: " << std::fixed << xCorrRange_.first * 1000. << ", x_corr_max: " << xCorrRange_.second * 1000.
0188        << "\n";
0189     // print in um

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", true);
0213   desc.add<bool>("y_ali_final_slope_fixed", true);
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   // setting default sh_x values from config

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

0269   CTPPSRPAlignmentCorrectionsData finalResults;
0270   if (doXAli) {  // x alignment

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

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

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

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

0299     if (yAliFinalSlopeFixed_) {
0300       finalResults.addCorrections(yAliResultsSlopeFixed_);
0301     } else {
0302       finalResults.addCorrections(yAliResults_);
0303     }
0304   }
0305 
0306   // check if the results are within the reasonability ranges xCorrRange and yCorrRange

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

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

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

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

0376 
0377 // Builds graph from a vector of points (with errors).

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

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

0403     {
0404       // retrieve parent directory

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(), '_', '.');  // replace _ with .

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

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

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

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

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

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

0479   double S2_norm_best = 1E100;
0480 
0481   for (double sh = sh_min; sh <= sh_max; sh += cfg.x_ali_sh_step()) {
0482     // calculate chi^2

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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     CTPPSRPAlignmentCorrectionData rpResult_N(+b / 2., b_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0743     xAliRelResults_.setRPCorrection(sc.rp_N_.id_, rpResult_N);
0744     CTPPSRPAlignmentCorrectionData rpResult_F(-b / 2., b_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0745     xAliRelResults_.setRPCorrection(sc.rp_F_.id_, rpResult_F);
0746 
0747     // calculate the results with slope fixed

0748     ff_sl_fix->SetParameters(0., slope, 0.);
0749     ff_sl_fix->FixParameter(1, slope);
0750     ff_sl_fix->FixParameter(2, -sh_x_N);
0751     ff_sl_fix->SetLineColor(4);
0752     p_x_diffFN_vs_x_N->Fit(ff_sl_fix.get(), "Q+", "", x_min, x_max);
0753 
0754     const double b_fs = ff_sl_fix->GetParameter(0), b_fs_unc = ff_sl_fix->GetParError(0);
0755 
0756     CTPPSRPAlignmentCorrectionData rpResult_sl_fix_N(+b_fs / 2., b_fs_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0757     xAliRelResultsSlopeFixed_.setRPCorrection(sc.rp_N_.id_, rpResult_sl_fix_N);
0758     CTPPSRPAlignmentCorrectionData rpResult_sl_fix_F(-b_fs / 2., b_fs_unc / 2., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.);
0759     xAliRelResultsSlopeFixed_.setRPCorrection(sc.rp_F_.id_, rpResult_sl_fix_F);
0760 
0761     edm::LogInfo("PPSAlignmentHarvester")
0762         << "[x_alignment_relative] " << std::fixed << std::setprecision(3) << "ff: " << ff->GetParameter(0) << " + "
0763         << ff->GetParameter(1) << " * (x - " << ff->GetParameter(2) << "), ff_sl_fix: " << ff_sl_fix->GetParameter(0)
0764         << " + " << ff_sl_fix->GetParameter(1) << " * (x - " << ff_sl_fix->GetParameter(2) << ")";
0765 
0766     // rebook the diffFN plot in the harvester

0767     iBooker.setCurrentFolder(dqmDir_ + "/harvester/x_alignment_relative/" + sc.name_);
0768     iBooker.bookProfile("p_x_diffFN_vs_x_N", p_x_diffFN_vs_x_N);
0769 
0770     if (debug_) {
0771       p_x_diffFN_vs_x_N->Write("p_x_diffFN_vs_x_N");
0772 
0773       auto g_results = std::make_unique<TGraph>();
0774       g_results->SetPoint(0, sh_x_N, 0.);
0775       g_results->SetPoint(1, a, a_unc);
0776       g_results->SetPoint(2, b, b_unc);
0777       g_results->SetPoint(3, b_fs, b_fs_unc);
0778       g_results->Write("g_results");
0779     }
0780   }
0781 
0782   // write results

0783   edm::LogInfo("PPSAlignmentHarvester") << seqPos << ": x_alignment_relative:\n"
0784                                         << xAliRelResults_ << seqPos + 1 << ": x_alignment_relative_sl_fix:\n"
0785                                         << xAliRelResultsSlopeFixed_;
0786 
0787   if (textResultsFile_.is_open()) {
0788     textResultsFile_ << seqPos << ": x_alignment_relative:\n" << xAliRelResults_ << "\n";
0789     textResultsFile_ << seqPos << ": x_alignment_relative_sl_fix:\n" << xAliRelResultsSlopeFixed_ << "\n\n";
0790   }
0791 }
0792 
0793 // -------------------------------- y alignment methods --------------------------------

0794 
0795 double PPSAlignmentHarvester::findMax(const TF1* ff_fit) {
0796   const double mu = ff_fit->GetParameter(1);
0797   const double si = ff_fit->GetParameter(2);
0798 
0799   // unreasonable fit?

0800   if (si > 25. || std::fabs(mu) > 100.)
0801     return 1E100;
0802 
0803   double x_max = 1E100;
0804   double y_max = -1E100;
0805   for (double x = mu - si; x <= mu + si; x += 0.001) {
0806     double y = ff_fit->Eval(x);
0807     if (y > y_max) {
0808       x_max = x;
0809       y_max = y;
0810     }
0811   }
0812 
0813   return x_max;
0814 }
0815 
0816 TH1D* PPSAlignmentHarvester::buildModeGraph(DQMStore::IBooker& iBooker,
0817                                             const MonitorElement* h2_y_vs_x,
0818                                             const PPSAlignmentConfiguration& cfg,
0819                                             const PPSAlignmentConfiguration::RPConfig& rpc) {
0820   TDirectory* d_top = nullptr;
0821   if (debug_)
0822     d_top = gDirectory;
0823 
0824   auto ff_fit = std::make_unique<TF1>("ff_fit", "[0] * exp(-(x-[1])*(x-[1])/2./[2]/[2]) + [3] + [4]*x");
0825 
0826   int h_n = h2_y_vs_x->getNbinsX();
0827   double diff = h2_y_vs_x->getTH2D()->GetXaxis()->GetBinWidth(1) / 2.;
0828   auto h_mode = iBooker.book1DD("mode",
0829                                 ";x (mm); mode of y (mm)",
0830                                 h_n,
0831                                 h2_y_vs_x->getTH2D()->GetXaxis()->GetBinCenter(1) - diff,
0832                                 h2_y_vs_x->getTH2D()->GetXaxis()->GetBinCenter(h_n) + diff);
0833 
0834   // find mode for each bin

0835   for (int bix = 1; bix <= h_n; bix++) {
0836     const double x = h2_y_vs_x->getTH2D()->GetXaxis()->GetBinCenter(bix);
0837 
0838     char buf[100];
0839     sprintf(buf, "h_y_x=%.3f", x);
0840     TH1D* h_y = h2_y_vs_x->getTH2D()->ProjectionY(buf, bix, bix);
0841 
0842     if (h_y->GetEntries() < cfg.multSelProjYMinEntries())
0843       continue;
0844 
0845     if (debug_) {
0846       sprintf(buf, "x=%.3f", x);
0847       gDirectory = d_top->mkdir(buf);
0848     }
0849 
0850     double conMax = -1.;
0851     double conMax_x = 0.;
0852     for (int biy = 1; biy < h_y->GetNbinsX(); biy++) {
0853       if (h_y->GetBinContent(biy) > conMax) {
0854         conMax = h_y->GetBinContent(biy);
0855         conMax_x = h_y->GetBinCenter(biy);
0856       }
0857     }
0858 
0859     ff_fit->SetParameters(conMax, conMax_x, h_y->GetRMS() * 0.75, 0., 0.);
0860     ff_fit->FixParameter(4, 0.);
0861 
0862     double x_min = rpc.x_min_fit_mode_, x_max = rpc.x_max_fit_mode_;
0863     h_y->Fit(ff_fit.get(), "Q", "", x_min, x_max);
0864 
0865     ff_fit->ReleaseParameter(4);
0866     double w = std::min(4., 2. * ff_fit->GetParameter(2));
0867     x_min = ff_fit->GetParameter(1) - w;
0868     x_max = std::min(rpc.y_max_fit_mode_, ff_fit->GetParameter(1) + w);
0869 
0870     h_y->Fit(ff_fit.get(), "Q", "", x_min, x_max);
0871 
0872     if (debug_)
0873       h_y->Write("h_y");
0874 
0875     double y_mode = findMax(ff_fit.get());
0876     const double y_mode_fit_unc = ff_fit->GetParameter(2) / 10;
0877     const double y_mode_sys_unc = cfg.y_mode_sys_unc();
0878     double y_mode_unc = std::sqrt(y_mode_fit_unc * y_mode_fit_unc + y_mode_sys_unc * y_mode_sys_unc);
0879 
0880     const double chiSqThreshold = cfg.chiSqThreshold();
0881 
0882     const bool valid =
0883         !(std::fabs(y_mode_unc) > cfg.y_mode_unc_max_valid() || std::fabs(y_mode) > cfg.y_mode_max_valid() ||
0884           ff_fit->GetChisquare() / ff_fit->GetNDF() > chiSqThreshold);
0885 
0886     if (debug_) {
0887       auto g_data = std::make_unique<TGraph>();
0888       g_data->SetPoint(0, y_mode, y_mode_unc);
0889       g_data->SetPoint(1, ff_fit->GetChisquare(), ff_fit->GetNDF());
0890       g_data->SetPoint(2, valid, 0.);
0891       g_data->Write("g_data");
0892     }
0893 
0894     if (!valid)
0895       continue;
0896 
0897     h_mode->Fill(x, y_mode);
0898     h_mode->setBinError(bix, y_mode_unc);
0899   }
0900 
0901   return h_mode->getTH1D();
0902 }
0903 
0904 void PPSAlignmentHarvester::yAlignment(DQMStore::IBooker& iBooker,
0905                                        DQMStore::IGetter& iGetter,
0906                                        const PPSAlignmentConfiguration& cfg) {
0907   TDirectory* yAliDir = nullptr;
0908   if (debug_)
0909     yAliDir = debugFile_->mkdir((std::to_string(seqPos) + ": y_alignment").c_str());
0910 
0911   auto ff = std::make_unique<TF1>("ff", "[0] + [1]*(x - [2])");
0912   auto ff_sl_fix = std::make_unique<TF1>("ff_sl_fix", "[0] + [1]*(x - [2])");
0913 
0914   // processing

0915   for (const auto& sc : {cfg.sectorConfig45(), cfg.sectorConfig56()}) {
0916     for (const auto& rpc : {sc.rp_F_, sc.rp_N_}) {
0917       TDirectory* rpDir = nullptr;
0918       if (debug_) {
0919         rpDir = yAliDir->mkdir(rpc.name_.c_str());
0920         gDirectory = rpDir->mkdir("x");
0921       }
0922 
0923       auto* h2_y_vs_x =
0924           iGetter.get(dqmDir_ + "/worker/" + sc.name_ + "/multiplicity selection/" + rpc.name_ + "/h2_y_vs_x");
0925 
0926       if (h2_y_vs_x == nullptr) {
0927         edm::LogWarning("PPSAlignmentHarvester") << "[y_alignment] " << rpc.name_ << ": cannot load data, skipping";
0928         continue;
0929       }
0930 
0931       iBooker.setCurrentFolder(dqmDir_ + "/harvester/y alignment/" + rpc.name_);
0932       auto* h_y_cen_vs_x = buildModeGraph(iBooker, h2_y_vs_x, cfg, rpc);
0933 
0934       if ((unsigned int)h_y_cen_vs_x->GetEntries() < cfg.modeGraphMinN()) {
0935         edm::LogWarning("PPSAlignmentHarvester")
0936             << "[y_alignment] " << rpc.name_ << ": insufficient data, skipping (mode graph "
0937             << h_y_cen_vs_x->GetEntries() << "/" << cfg.modeGraphMinN() << ")";
0938         continue;
0939       }
0940 
0941       const double x_min = cfg.alignment_y_ranges().at(rpc.id_).x_min_;
0942       const double x_max = cfg.alignment_y_ranges().at(rpc.id_).x_max_;
0943 
0944       const double sh_x = sh_x_map_[rpc.id_];
0945       double slope = rpc.slope_;
0946 
0947       // calculate the results without slope fixed

0948       ff->SetParameters(0., 0., 0.);
0949       ff->FixParameter(2, -sh_x);
0950       ff->SetLineColor(2);
0951       h_y_cen_vs_x->Fit(ff.get(), "Q", "", x_min, x_max);
0952 
0953       const double a = ff->GetParameter(1), a_unc = ff->GetParError(1);
0954       const double b = ff->GetParameter(0), b_unc = ff->GetParError(0);
0955 
0956       edm::LogInfo("PPSAlignmentHarvester")
0957           << "[y_alignment] " << rpc.name_ << ":\n"
0958           << std::fixed << std::setprecision(3) << "    x_min = " << x_min << ", x_max = " << x_max << "\n"
0959           << "    sh_x = " << sh_x << ", slope (fix) = " << slope << ", slope (fitted) = " << a;
0960 
0961       CTPPSRPAlignmentCorrectionData rpResult(0., 0., b, b_unc, 0., 0., 0., 0., 0., 0., 0., 0.);
0962       yAliResults_.setRPCorrection(rpc.id_, rpResult);
0963 
0964       // calculate the results with slope fixed

0965       ff_sl_fix->SetParameters(0., 0., 0.);
0966       ff_sl_fix->FixParameter(1, slope);
0967       ff_sl_fix->FixParameter(2, -sh_x);
0968       ff_sl_fix->SetLineColor(4);
0969       h_y_cen_vs_x->Fit(ff_sl_fix.get(), "Q+", "", x_min, x_max);
0970 
0971       const double b_fs = ff_sl_fix->GetParameter(0), b_fs_unc = ff_sl_fix->GetParError(0);
0972 
0973       CTPPSRPAlignmentCorrectionData rpResult_sl_fix(0., 0., b_fs, b_fs_unc, 0., 0., 0., 0., 0., 0., 0., 0.);
0974       yAliResultsSlopeFixed_.setRPCorrection(rpc.id_, rpResult_sl_fix);
0975 
0976       edm::LogInfo("PPSAlignmentHarvester")
0977           << "[y_alignment] " << std::fixed << std::setprecision(3) << "ff: " << ff->GetParameter(0) << " + "
0978           << ff->GetParameter(1) << " * (x - " << ff->GetParameter(2) << "), ff_sl_fix: " << ff_sl_fix->GetParameter(0)
0979           << " + " << ff_sl_fix->GetParameter(1) << " * (x - " << ff_sl_fix->GetParameter(2) << ")";
0980 
0981       if (debug_) {
0982         gDirectory = rpDir;
0983 
0984         h_y_cen_vs_x->SetTitle(";x (mm); mode of y (mm)");
0985         h_y_cen_vs_x->Write("h_y_cen_vs_x");
0986 
0987         auto g_results = std::make_unique<TGraph>();
0988         g_results->SetPoint(0, sh_x, 0.);
0989         g_results->SetPoint(1, a, a_unc);
0990         g_results->SetPoint(2, b, b_unc);
0991         g_results->SetPoint(3, b_fs, b_fs_unc);
0992         g_results->Write("g_results");
0993       }
0994     }
0995   }
0996 
0997   // write results

0998   edm::LogInfo("PPSAlignmentHarvester") << seqPos << ": y_alignment:\n"
0999                                         << yAliResults_ << seqPos << ": y_alignment_sl_fix:\n"
1000                                         << yAliResultsSlopeFixed_;
1001 
1002   if (textResultsFile_.is_open()) {
1003     textResultsFile_ << seqPos << ": y_alignment:\n" << yAliResults_ << "\n";
1004     textResultsFile_ << seqPos << ": y_alignment_sl_fix:\n" << yAliResultsSlopeFixed_ << "\n\n";
1005   }
1006 }
1007 
1008 // -------------------------------- other methods --------------------------------

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

1011 void PPSAlignmentHarvester::writeCutPlot(
1012     TH2D* h, const double a, const double c, const double n_si, const double si, const std::string& label) {
1013   auto canvas = std::make_unique<TCanvas>();
1014   canvas->SetName(label.c_str());
1015   canvas->SetLogz(1);
1016 
1017   h->Draw("colz");
1018 
1019   double x_min = -30.;
1020   double x_max = 30.;
1021 
1022   auto g_up = std::make_unique<TGraph>();
1023   g_up->SetName("g_up");
1024   g_up->SetPoint(0, x_min, -a * x_min - c + n_si * si);
1025   g_up->SetPoint(1, x_max, -a * x_max - c + n_si * si);
1026   g_up->SetLineColor(1);
1027   g_up->Draw("l");
1028 
1029   auto g_down = std::make_unique<TGraph>();
1030   g_down->SetName("g_down");
1031   g_down->SetPoint(0, x_min, -a * x_min - c - n_si * si);
1032   g_down->SetPoint(1, x_max, -a * x_max - c - n_si * si);
1033   g_down->SetLineColor(1);
1034   g_down->Draw("l");
1035 
1036   canvas->Write();
1037 }
1038 
1039 // Points in TGraph should be sorted (TGraph::Sort())

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

1041 std::unique_ptr<TH1D> PPSAlignmentHarvester::getTH1DFromTGraphErrors(
1042     TGraphErrors* graph, const std::string& title, const std::string& labels, int n, double binWidth, double min) {
1043   std::unique_ptr<TH1D> hist;
1044   if (n == 0) {
1045     hist = std::make_unique<TH1D>(title.c_str(), labels.c_str(), 0, -10., 10.);
1046   } else if (n == 1) {
1047     hist = std::make_unique<TH1D>(title.c_str(), labels.c_str(), 1, graph->GetPointX(0) - 5., graph->GetPointX(0) + 5.);
1048   } else {
1049     n = n == -1 ? graph->GetN() : n;
1050     binWidth = binWidth == -1 ? graph->GetPointX(1) - graph->GetPointX(0) : binWidth;
1051     double diff = binWidth / 2.;
1052     min = min == -1 ? graph->GetPointX(0) - diff : min;
1053     double max = min + n * binWidth;
1054     hist = std::make_unique<TH1D>(title.c_str(), labels.c_str(), n, min, max);
1055   }
1056 
1057   for (int i = 0; i < graph->GetN(); i++) {
1058     double x, y;
1059     graph->GetPoint(i, x, y);
1060     hist->Fill(x, y);
1061     hist->SetBinError(hist->GetXaxis()->FindBin(x), graph->GetErrorY(i));
1062   }
1063   return hist;
1064 }
1065 
1066 // Get Long 32-bit detector ID from short 3-digit ID

1067 CTPPSRPAlignmentCorrectionsData PPSAlignmentHarvester::getLongIdResults(CTPPSRPAlignmentCorrectionsData shortIdResults) {
1068   CTPPSRPAlignmentCorrectionsData longIdResults;
1069   for (const auto& [shortId, correction] : shortIdResults.getRPMap()) {
1070     unsigned int arm = shortId / 100;
1071     unsigned int station = (shortId / 10) % 10;
1072     unsigned int rp = shortId % 10;
1073 
1074     uint32_t longDetId = detectorId_ << 28 | subdetectorId_ << 25 | arm << 24 | station << 22 | rp << 19;
1075 
1076     longIdResults.addRPCorrection(longDetId, correction);
1077   }
1078 
1079   return longIdResults;
1080 }
1081 
1082 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!