File indexing completed on 2023-03-17 11:27:02
0001
0002
0003
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
0128 if (!rp_plots_.empty())
0129 return;
0130
0131
0132 const auto& opticalFunctions = iSetup.getData(opticsESToken_);
0133
0134
0135 if (opticalFunctions.empty())
0136 return;
0137
0138
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;
0150 const double th_ep = 1E-6;
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
0188 for (const auto& ap : arm_plots_) {
0189
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);