CMS 3D CMS Logo

CTPPSOpticsPlotter.cc
Go to the documentation of this file.
1 /****************************************************************************
2  * Authors:
3  * Jan Kašpar
4  ****************************************************************************/
5 
10 
13 
15 
18 
19 #include "TFile.h"
20 #include "TGraph.h"
21 
22 #include <map>
23 #include <string>
24 
25 //----------------------------------------------------------------------------------------------------
26 
28 public:
29  explicit CTPPSOpticsPlotter(const edm::ParameterSet&);
30 
31 private:
32  void analyze(const edm::Event&, const edm::EventSetup&) override;
33  void endJob() override;
34 
36 
37  unsigned int rpId_45_N_, rpId_45_F_;
38  unsigned int rpId_56_N_, rpId_56_F_;
39 
41 
42  struct RPPlots {
43  std::unique_ptr<TGraph> g_v_x_vs_xi, g_L_x_vs_xi, g_x_D_vs_xi;
44  std::unique_ptr<TGraph> g_v_y_vs_xi, g_L_y_vs_xi, g_y_D_vs_xi;
45  std::unique_ptr<TGraph> h_y_vs_x_disp;
46 
48  : g_v_x_vs_xi(new TGraph),
49  g_L_x_vs_xi(new TGraph),
50  g_x_D_vs_xi(new TGraph),
51  g_v_y_vs_xi(new TGraph),
52  g_L_y_vs_xi(new TGraph),
53  g_y_D_vs_xi(new TGraph),
54  h_y_vs_x_disp(new TGraph) {}
55 
56  void write() const {
57  g_v_x_vs_xi->SetTitle(";xi;v_{x}");
58  g_v_x_vs_xi->Write("g_v_x_vs_xi");
59 
60  g_L_x_vs_xi->SetTitle(";xi;L_{x} (cm)");
61  g_L_x_vs_xi->Write("g_L_x_vs_xi");
62 
63  g_x_D_vs_xi->SetTitle(";xi;x_{D} (cm)");
64  g_x_D_vs_xi->Write("g_x_D_vs_xi");
65 
66  g_v_y_vs_xi->SetTitle(";xi;v_{y}");
67  g_v_y_vs_xi->Write("g_v_y_vs_xi");
68 
69  g_L_y_vs_xi->SetTitle(";xi;L_{y} (cm)");
70  g_L_y_vs_xi->Write("g_L_y_vs_xi");
71 
72  g_y_D_vs_xi->SetTitle(";xi;y_{D} (cm)");
73  g_y_D_vs_xi->Write("g_y_D_vs_xi");
74 
75  h_y_vs_x_disp->SetTitle(";x (cm);y (cm)");
76  h_y_vs_x_disp->Write("h_y_vs_x_disp");
77  }
78  };
79 
80  std::map<unsigned int, RPPlots> rp_plots_;
81 
82  struct ArmPlots {
83  unsigned int id_N, id_F;
84 
85  std::unique_ptr<TGraph> g_de_x_vs_x_disp, g_de_y_vs_x_disp;
86 
87  ArmPlots() : g_de_x_vs_x_disp(new TGraph), g_de_y_vs_x_disp(new TGraph) {}
88 
89  void write() const {
90  g_de_x_vs_x_disp->SetTitle(";x_N (cm);x_F - x_N (cm)");
91  g_de_x_vs_x_disp->Write("g_de_x_vs_x_disp");
92 
93  g_de_y_vs_x_disp->SetTitle(";x_N (cm);y_F - y_N (cm)");
94  g_de_y_vs_x_disp->Write("g_de_y_vs_x_disp");
95  }
96  };
97 
98  std::map<unsigned int, ArmPlots> arm_plots_;
99 };
100 
101 //----------------------------------------------------------------------------------------------------
102 
103 using namespace std;
104 using namespace edm;
105 
106 //----------------------------------------------------------------------------------------------------
107 
109  : opticsESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("opticsLabel")))),
110 
111  rpId_45_N_(iConfig.getParameter<unsigned int>("rpId_45_N")),
112  rpId_45_F_(iConfig.getParameter<unsigned int>("rpId_45_F")),
113  rpId_56_N_(iConfig.getParameter<unsigned int>("rpId_56_N")),
114  rpId_56_F_(iConfig.getParameter<unsigned int>("rpId_56_F")),
115 
116  outputFile_(iConfig.getParameter<string>("outputFile")) {
117  arm_plots_[0].id_N = rpId_45_N_;
118  arm_plots_[0].id_F = rpId_45_F_;
119 
120  arm_plots_[1].id_N = rpId_56_N_;
121  arm_plots_[1].id_F = rpId_56_F_;
122 }
123 
124 //----------------------------------------------------------------------------------------------------
125 
127  // stop if plots already made
128  if (!rp_plots_.empty())
129  return;
130 
131  // get conditions
132  const auto& opticalFunctions = iSetup.getData(opticsESToken_);
133 
134  // stop if conditions invalid
135  if (opticalFunctions.empty())
136  return;
137 
138  // make per-RP plots
139  for (const auto& it : opticalFunctions) {
140  CTPPSDetId rpId(it.first);
141  unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
142 
143  auto& pl = rp_plots_[rpDecId];
144 
145  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
147  it.second.transport(k_in_beam, k_out_beam);
148 
149  const double vtx_ep = 1E-4; // cm
150  const double th_ep = 1E-6; // rad
151 
152  for (double xi = 0.; xi < 0.30001; xi += 0.001) {
153  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
155  it.second.transport(k_in_xi, k_out_xi);
156 
157  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_vtx_x = {vtx_ep, 0., 0., 0., xi};
159  it.second.transport(k_in_xi_vtx_x, k_out_xi_vtx_x);
160 
161  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_th_x = {0., th_ep, 0., 0., xi};
163  it.second.transport(k_in_xi_th_x, k_out_xi_th_x);
164 
165  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_vtx_y = {0., 0., vtx_ep, 0., xi};
167  it.second.transport(k_in_xi_vtx_y, k_out_xi_vtx_y);
168 
169  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi_th_y = {0., 0., 0., th_ep, xi};
171  it.second.transport(k_in_xi_th_y, k_out_xi_th_y);
172 
173  int idx = pl.g_v_x_vs_xi->GetN();
174 
175  pl.g_v_x_vs_xi->SetPoint(idx, xi, (k_out_xi_vtx_x.x - k_out_xi.x) / vtx_ep);
176  pl.g_L_x_vs_xi->SetPoint(idx, xi, (k_out_xi_th_x.x - k_out_xi.x) / th_ep);
177  pl.g_x_D_vs_xi->SetPoint(idx, xi, k_out_xi.x - k_out_beam.x);
178 
179  pl.g_v_y_vs_xi->SetPoint(idx, xi, (k_out_xi_vtx_y.y - k_out_xi.y) / vtx_ep);
180  pl.g_L_y_vs_xi->SetPoint(idx, xi, (k_out_xi_th_y.y - k_out_xi.y) / th_ep);
181  pl.g_y_D_vs_xi->SetPoint(idx, xi, k_out_xi.y - k_out_beam.y);
182 
183  pl.h_y_vs_x_disp->SetPoint(idx, k_out_xi.x - k_out_beam.x, k_out_xi.y - k_out_beam.y);
184  }
185  }
186 
187  // make per-arm plots
188  for (const auto& ap : arm_plots_) {
189  // find optics objects
190  const LHCInterpolatedOpticalFunctionsSet *opt_N = nullptr, *opt_F = nullptr;
191 
192  for (const auto& it : opticalFunctions) {
193  CTPPSDetId rpId(it.first);
194  unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
195 
196  if (rpDecId == ap.second.id_N)
197  opt_N = &it.second;
198  if (rpDecId == ap.second.id_F)
199  opt_F = &it.second;
200  }
201 
202  if (!opt_N || !opt_F) {
203  edm::LogError("CTPPSOpticsPlotter::analyze") << "Cannot find optics objects for arm " << ap.first;
204  continue;
205  }
206 
207  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_beam = {0., 0., 0., 0., 0.};
208 
209  LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_beam_N, k_out_beam_F;
210  opt_N->transport(k_in_beam, k_out_beam_N);
211  opt_F->transport(k_in_beam, k_out_beam_F);
212 
213  for (double xi = 0.; xi < 0.30001; xi += 0.001) {
214  LHCInterpolatedOpticalFunctionsSet::Kinematics k_in_xi = {0., 0., 0., 0., xi};
215 
216  LHCInterpolatedOpticalFunctionsSet::Kinematics k_out_xi_N, k_out_xi_F;
217  opt_N->transport(k_in_xi, k_out_xi_N);
218  opt_F->transport(k_in_xi, k_out_xi_F);
219 
220  int idx = ap.second.g_de_x_vs_x_disp->GetN();
221 
222  ap.second.g_de_x_vs_x_disp->SetPoint(
223  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));
224  ap.second.g_de_y_vs_x_disp->SetPoint(
225  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));
226  }
227  }
228 }
229 
230 //----------------------------------------------------------------------------------------------------
231 
233  auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
234 
235  for (const auto& p : rp_plots_) {
236  gDirectory = f_out->mkdir(Form("%u", p.first));
237  p.second.write();
238  }
239 
240  for (const auto& p : arm_plots_) {
241  gDirectory = f_out->mkdir(Form("arm %u", p.first));
242  p.second.write();
243  }
244 }
245 
246 //----------------------------------------------------------------------------------------------------
247 
void endJob() override
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
void analyze(const edm::Event &, const edm::EventSetup &) override
std::unique_ptr< TGraph > g_L_x_vs_xi
CTPPSOpticsPlotter(const edm::ParameterSet &)
std::unique_ptr< TGraph > g_x_D_vs_xi
uint32_t arm() const
Definition: CTPPSDetId.h:57
Log< level::Error, false > LogError
std::map< unsigned int, ArmPlots > arm_plots_
std::unique_ptr< TGraph > g_L_y_vs_xi
std::unique_ptr< TGraph > g_de_x_vs_x_disp
int iEvent
Definition: GenABIO.cc:224
std::unique_ptr< TGraph > g_de_y_vs_x_disp
void transport(const Kinematics &input, Kinematics &output, bool calculateAngles=false) const
transports proton according to the splines
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
bool getData(T &iHolder) const
Definition: EventSetup.h:122
std::unique_ptr< TGraph > g_v_y_vs_xi
std::map< unsigned int, RPPlots > rp_plots_
std::unique_ptr< TGraph > g_v_x_vs_xi
std::unique_ptr< TGraph > h_y_vs_x_disp
uint32_t rp() const
Definition: CTPPSDetId.h:71
std::unique_ptr< TGraph > g_y_D_vs_xi
edm::ESGetToken< LHCInterpolatedOpticalFunctionsSetCollection, CTPPSInterpolatedOpticsRcd > opticsESToken_
uint32_t station() const
Definition: CTPPSDetId.h:64
Set of optical functions corresponding to one scoring plane along LHC, including splines for interpol...
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:32
HLT enums.