CTPPS/plugins/CTPPSOpticsPlotter.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013
0014 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0015
0016 #include "CondFormats/DataRecord/interface/CTPPSInterpolatedOpticsRcd.h"
0017 #include "CondFormats/PPSObjects/interface/LHCInterpolatedOpticalFunctionsSetCollection.h"
0018
0019 #include "TFile.h"
0020 #include "TGraph.h"
0021
0022 #include <map>
0023 #include <string>
0024
0025 //----------------------------------------------------------------------------------------------------
0026
0027 class CTPPSOpticsPlotter : public edm::one::EDAnalyzer<> {
0028 public:
0029   explicit CTPPSOpticsPlotter(const edm::ParameterSet&);
0030
0031 private:
0032   void analyze(const edm::Event&, const edm::EventSetup&) override;
0033   void endJob() override;
0034
0035   edm::ESGetToken<LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd> opticsESToken_;
0036
0037   unsigned int rpId_45_N_, rpId_45_F_;
0038   unsigned int rpId_56_N_, rpId_56_F_;
0039
0040   std::string outputFile_;
0041
0042   struct RPPlots {
0043     std::unique_ptr<TGraph> g_v_x_vs_xi, g_L_x_vs_xi, g_x_D_vs_xi;
0044     std::unique_ptr<TGraph> g_v_y_vs_xi, g_L_y_vs_xi, g_y_D_vs_xi;
0045     std::unique_ptr<TGraph> h_y_vs_x_disp;
0046
0047     RPPlots()
0048         : g_v_x_vs_xi(new TGraph),
0049           g_L_x_vs_xi(new TGraph),
0050           g_x_D_vs_xi(new TGraph),
0051           g_v_y_vs_xi(new TGraph),
0052           g_L_y_vs_xi(new TGraph),
0053           g_y_D_vs_xi(new TGraph),
0054           h_y_vs_x_disp(new TGraph) {}
0055
0056     void write() const {
0057       g_v_x_vs_xi->SetTitle(";xi;v_{x}");
0058       g_v_x_vs_xi->Write("g_v_x_vs_xi");
0059
0060       g_L_x_vs_xi->SetTitle(";xi;L_{x}   (cm)");
0061       g_L_x_vs_xi->Write("g_L_x_vs_xi");
0062
0063       g_x_D_vs_xi->SetTitle(";xi;x_{D}   (cm)");
0064       g_x_D_vs_xi->Write("g_x_D_vs_xi");
0065
0066       g_v_y_vs_xi->SetTitle(";xi;v_{y}");
0067       g_v_y_vs_xi->Write("g_v_y_vs_xi");
0068
0069       g_L_y_vs_xi->SetTitle(";xi;L_{y}   (cm)");
0070       g_L_y_vs_xi->Write("g_L_y_vs_xi");
0071
0072       g_y_D_vs_xi->SetTitle(";xi;y_{D}   (cm)");
0073       g_y_D_vs_xi->Write("g_y_D_vs_xi");
0074
0075       h_y_vs_x_disp->SetTitle(";x   (cm);y   (cm)");
0076       h_y_vs_x_disp->Write("h_y_vs_x_disp");
0077     }
0078   };
0079
0080   std::map<unsigned int, RPPlots> rp_plots_;
0081
0082   struct ArmPlots {
0083     unsigned int id_N, id_F;
0084
0085     std::unique_ptr<TGraph> g_de_x_vs_x_disp, g_de_y_vs_x_disp;
0086
0087     ArmPlots() : g_de_x_vs_x_disp(new TGraph), g_de_y_vs_x_disp(new TGraph) {}
0088
0089     void write() const {
0090       g_de_x_vs_x_disp->SetTitle(";x_N   (cm);x_F - x_N   (cm)");
0091       g_de_x_vs_x_disp->Write("g_de_x_vs_x_disp");
0092
0093       g_de_y_vs_x_disp->SetTitle(";x_N   (cm);y_F - y_N   (cm)");
0094       g_de_y_vs_x_disp->Write("g_de_y_vs_x_disp");
0095     }
0096   };
0097
0098   std::map<unsigned int, ArmPlots> arm_plots_;
0099 };
0100
0101 //----------------------------------------------------------------------------------------------------
0102
0103 using namespace std;
0104 using namespace edm;
0105
0106 //----------------------------------------------------------------------------------------------------
0107
0108 CTPPSOpticsPlotter::CTPPSOpticsPlotter(const edm::ParameterSet& iConfig)
0109     : opticsESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
0110
0111       rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
0112       rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
0113       rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
0114       rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
0115
0116       outputFile_(iConfig.getParameter<string>("outputFile")) {
0117   arm_plots_[0].id_N = rpId_45_N_;
0118   arm_plots_[0].id_F = rpId_45_F_;
0119
0120   arm_plots_[1].id_N = rpId_56_N_;
0121   arm_plots_[1].id_F = rpId_56_F_;
0122 }
0123
0124 //----------------------------------------------------------------------------------------------------
0125
0126 void CTPPSOpticsPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0127   // stop if plots already made
0128   if (!rp_plots_.empty())
0129     return;
0130
0131   // get conditions
0132   const auto& opticalFunctions = iSetup.getData(opticsESToken_);
0133
0134   // stop if conditions invalid
0135   if (opticalFunctions.empty())
0136     return;
0137
0138   // make per-RP plots
0139   for (const auto& it : opticalFunctions) {
0140     CTPPSDetId rpId(it.first);
0141     unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0142
0143     auto& pl = rp_plots_[rpDecId];
0144
0145     LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
0146     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam;
0147     it.second.transport(k_in_beam, k_out_beam);
0148
0149     const double vtx_ep = 1E-4;  // cm
0150     const double th_ep = 1E-6;   // rad
0151
0152     for (double xi = 0.; xi < 0.30001; xi += 0.001) {
0153       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
0154       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi;
0155       it.second.transport(k_in_xi, k_out_xi);
0156
0157       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_vtx_x = {vtx_ep, 0., 0., 0., xi};
0158       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_vtx_x;
0159       it.second.transport(k_in_xi_vtx_x, k_out_xi_vtx_x);
0160
0161       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_th_x = {0., th_ep, 0., 0., xi};
0162       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_th_x;
0163       it.second.transport(k_in_xi_th_x, k_out_xi_th_x);
0164
0165       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_vtx_y = {0., 0., vtx_ep, 0., xi};
0166       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_vtx_y;
0167       it.second.transport(k_in_xi_vtx_y, k_out_xi_vtx_y);
0168
0169       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_th_y = {0., 0., 0., th_ep, xi};
0170       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_th_y;
0171       it.second.transport(k_in_xi_th_y, k_out_xi_th_y);
0172
0173       int idx = pl.g_v_x_vs_xi->GetN();
0174
0175       pl.g_v_x_vs_xi->SetPoint(idx, xi, (k_out_xi_vtx_x.x - k_out_xi.x) / vtx_ep);
0176       pl.g_L_x_vs_xi->SetPoint(idx, xi, (k_out_xi_th_x.x - k_out_xi.x) / th_ep);
0177       pl.g_x_D_vs_xi->SetPoint(idx, xi, k_out_xi.x - k_out_beam.x);
0178
0179       pl.g_v_y_vs_xi->SetPoint(idx, xi, (k_out_xi_vtx_y.y - k_out_xi.y) / vtx_ep);
0180       pl.g_L_y_vs_xi->SetPoint(idx, xi, (k_out_xi_th_y.y - k_out_xi.y) / th_ep);
0181       pl.g_y_D_vs_xi->SetPoint(idx, xi, k_out_xi.y - k_out_beam.y);
0182
0183       pl.h_y_vs_x_disp->SetPoint(idx, k_out_xi.x - k_out_beam.x, k_out_xi.y - k_out_beam.y);
0184     }
0185   }
0186
0187   // make per-arm plots
0188   for (const auto& ap : arm_plots_) {
0189     // find optics objects
0190     const LHCInterpolatedOpticalFunctionsSet *opt_N = nullptr, *opt_F = nullptr;
0191
0192     for (const auto& it : opticalFunctions) {
0193       CTPPSDetId rpId(it.first);
0194       unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0195
0196       if (rpDecId == ap.second.id_N)
0197         opt_N = &it.second;
0198       if (rpDecId == ap.second.id_F)
0199         opt_F = &it.second;
0200     }
0201
0202     if (!opt_N || !opt_F) {
0203       edm::LogError("CTPPSOpticsPlotter::analyze") << "Cannot find optics objects for arm " << ap.first;
0204       continue;
0205     }
0206
0207     LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
0208
0209     LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam_N, k_out_beam_F;
0210     opt_N->transport(k_in_beam, k_out_beam_N);
0211     opt_F->transport(k_in_beam, k_out_beam_F);
0212
0213     for (double xi = 0.; xi < 0.30001; xi += 0.001) {
0214       LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
0215
0216       LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_N, k_out_xi_F;
0217       opt_N->transport(k_in_xi, k_out_xi_N);
0218       opt_F->transport(k_in_xi, k_out_xi_F);
0219
0220       int idx = ap.second.g_de_x_vs_x_disp->GetN();
0221
0222       ap.second.g_de_x_vs_x_disp->SetPoint(
0223           idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.x - k_out_beam_F.x) - (k_out_xi_N.x - k_out_beam_N.x));
0224       ap.second.g_de_y_vs_x_disp->SetPoint(
0225           idx, k_out_xi_N.x - k_out_beam_N.x, (k_out_xi_F.y - k_out_beam_F.y) - (k_out_xi_N.y - k_out_beam_N.y));
0226     }
0227   }
0228 }
0229
0230 //----------------------------------------------------------------------------------------------------
0231
0232 void CTPPSOpticsPlotter::endJob() {
0233   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0234
0235   for (const auto& p : rp_plots_) {
0236     gDirectory = f_out->mkdir(Form("%u", p.first));
0237     p.second.write();
0238   }
0239
0240   for (const auto& p : arm_plots_) {
0241     gDirectory = f_out->mkdir(Form("arm %u", p.first));
0242     p.second.write();
0243   }
0244 }
0245
0246 //----------------------------------------------------------------------------------------------------
0247
0248 DEFINE_FWK_MODULE(CTPPSOpticsPlotter);