CMS 3D CMS Logo

CTPPSTrackDistributionPlotter.cc
Go to the documentation of this file.
1 /****************************************************************************
2  *
3  * This is a part of CTPPS validation software
4  * Authors:
5  * Jan Kašpar
6  * Laurent Forthomme
7  *
8  ****************************************************************************/
9 
12 
15 
17 
21 
22 #include "TFile.h"
23 #include "TH2D.h"
24 #include "TH1D.h"
25 #include "TProfile.h"
26 #include "TProfile2D.h"
27 #include "TGraph.h"
28 
29 #include <map>
30 #include <memory>
31 
32 //----------------------------------------------------------------------------------------------------
33 
35 public:
37 
39 
40 private:
41  void analyze(const edm::Event&, const edm::EventSetup&) override;
42  void endJob() override;
43 
45 
47 
49 
50  unsigned int events_total_;
51  std::map<unsigned int, unsigned int> events_per_arm_;
52 
53  struct RPPlots {
55 
56  std::unique_ptr<TH2D> h2_y_vs_x;
57  std::unique_ptr<TProfile> p_y_vs_x;
58  std::unique_ptr<TH1D> h_x;
59  std::unique_ptr<TH1D> h_y;
60  std::unique_ptr<TH1D> h_time;
61 
62  std::unique_ptr<TH2D> h2_de_x_vs_x; // "delta" refers to distance between two tracks in the same RP
63  std::unique_ptr<TH2D> h2_de_x_vs_y;
64  std::unique_ptr<TH2D> h2_de_y_vs_x;
65  std::unique_ptr<TH2D> h2_de_y_vs_y;
66 
68 
69  void init(bool pixel, double pitch) {
70  const double bin_size_x = (pixel) ? pitch * cos(18.4 / 180. * M_PI) : 100E-3;
71 
72  h2_y_vs_x = std::make_unique<TH2D>("", "", 300, -10., +70., 600, -30., +30.);
73  p_y_vs_x = std::make_unique<TProfile>("", "", 300, -10., +70.);
74 
75  int n_mi = ceil(10. / bin_size_x);
76  int n_pl = ceil(70. / bin_size_x);
77 
78  h_x = std::make_unique<TH1D>("", "", n_mi + n_pl, -n_mi * bin_size_x, +n_pl * bin_size_x);
79 
80  h_y = std::make_unique<TH1D>("", "", 300, -15., +15.);
81 
82  h_time = std::make_unique<TH1D>("", ";time", 500, -50., +50.);
83 
84  h2_de_x_vs_x =
85  std::make_unique<TH2D>("h2_de_x_vs_x", "h2_de_x_vs_x;x;distance in x axis", 300, -30., +30., 300, -3., +3.);
86  h2_de_x_vs_y =
87  std::make_unique<TH2D>("h2_de_x_vs_y", "h2_de_x_vs_y;y;distance in x axis", 300, -30., +30., 300, -3., +3.);
88  h2_de_y_vs_x =
89  std::make_unique<TH2D>("h2_de_y_vs_x", "h2_de_y_vs_x;x;distance in y axis", 300, -30., +30., 300, -3., +3.);
90  h2_de_y_vs_y =
91  std::make_unique<TH2D>("h2_de_y_vs_y", "h2_de_y_vs_y;y;distance in y axis", 300, -30., +30., 300, -3., +3.);
92 
93  initialized = true;
94  }
95 
96  void fillOneTrack(double x, double y, double time) {
97  h2_y_vs_x->Fill(x, y);
98  p_y_vs_x->Fill(x, y);
99  h_x->Fill(x);
100  h_y->Fill(y);
101  h_time->Fill(time);
102  }
103 
104  void write() const {
105  h2_y_vs_x->Write("h2_y_vs_x");
106  p_y_vs_x->Write("p_y_vs_x");
107  h_x->Write("h_x");
108  h_y->Write("h_y");
109  h_time->Write("h_time");
110 
111  h2_de_x_vs_x->Write("h2_de_x_vs_x");
112  h2_de_x_vs_y->Write("h2_de_x_vs_y");
113  h2_de_y_vs_x->Write("h2_de_y_vs_x");
114  h2_de_y_vs_y->Write("h2_de_y_vs_y");
115  }
116  };
117 
118  std::map<unsigned int, RPPlots> rpPlots;
119 
120  struct ArmPlots {
121  unsigned int rpId_N, rpId_F;
122 
123  std::unique_ptr<TH1D> h_de_x, h_de_y;
124  std::unique_ptr<TProfile> p_de_x_vs_x, p_de_y_vs_x;
125  std::unique_ptr<TProfile> p_de_x_vs_y, p_de_y_vs_y;
126  std::unique_ptr<TProfile2D> p2_de_x_vs_x_y, p2_de_y_vs_x_y;
127  std::unique_ptr<TH2D> h2_de_x_vs_x, h2_de_y_vs_x;
128  std::unique_ptr<TH2D> h2_de_y_vs_de_x;
129 
130  struct Stat {
131  unsigned int sN = 0, s1 = 0;
132  };
133 
134  std::map<unsigned int, std::map<unsigned int, Stat>> m_stat;
135 
137  : h_de_x(new TH1D("", ";x^{F} - x^{N}", 100, -1., +1.)),
138  h_de_y(new TH1D("", ";y^{F} - y^{N}", 100, -1., +1.)),
139 
140  p_de_x_vs_x(new TProfile("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., "s")),
141  p_de_y_vs_x(new TProfile("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., "s")),
142 
143  p_de_x_vs_y(new TProfile("", ";y^{N};x^{F} - x^{N}", 80, -20., +20., "s")),
144  p_de_y_vs_y(new TProfile("", ";y^{N};y^{F} - y^{N}", 80, -20., +20., "s")),
145 
146  p2_de_x_vs_x_y(new TProfile2D("", ";x;y;x^{F} - x^{N}", 80, 0., 40., 80, -20., +20., "s")),
147  p2_de_y_vs_x_y(new TProfile2D("", ";x;y;y^{F} - y^{N}", 80, 0., 40., 80, -20., +20., "s")),
148 
149  h2_de_x_vs_x(new TH2D("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., 100, -1., +1.)),
150  h2_de_y_vs_x(new TH2D("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., 100, -1., +1.)),
151  h2_de_y_vs_de_x(new TH2D("", ";x^{F} - x^{N};y^{F} - y^{N}", 100, -1., +1., 100, -1., +1.)) {}
152 
153  void fill(double x_N, double y_N, double x_F, double y_F) {
154  h_de_x->Fill(x_F - x_N);
155  h_de_y->Fill(y_F - y_N);
156 
157  p_de_x_vs_x->Fill(x_N, x_F - x_N);
158  p_de_y_vs_x->Fill(x_N, y_F - y_N);
159 
160  p_de_x_vs_y->Fill(y_N, x_F - x_N);
161  p_de_y_vs_y->Fill(y_N, y_F - y_N);
162 
163  p2_de_x_vs_x_y->Fill(x_N, y_N, x_F - x_N);
164  p2_de_y_vs_x_y->Fill(x_N, y_N, y_F - y_N);
165 
166  h2_de_x_vs_x->Fill(x_N, x_F - x_N);
167  h2_de_y_vs_x->Fill(x_N, y_F - y_N);
168 
169  h2_de_y_vs_de_x->Fill(x_F - x_N, y_F - y_N);
170  }
171 
172  void write() const {
173  h_de_x->Write("h_de_x");
174  h_de_y->Write("h_de_y");
175 
176  p_de_x_vs_x->Write("p_de_x_vs_x");
177  buildRMSHistogram(p_de_x_vs_x)->Write("h_rms_de_x_vs_x");
178  p_de_y_vs_x->Write("p_de_y_vs_x");
179  buildRMSHistogram(p_de_y_vs_x)->Write("h_rms_de_y_vs_x");
180 
181  p_de_x_vs_y->Write("p_de_x_vs_y");
182  buildRMSHistogram(p_de_x_vs_y)->Write("h_rms_de_x_vs_y");
183  p_de_y_vs_y->Write("p_de_y_vs_y");
184  buildRMSHistogram(p_de_y_vs_y)->Write("h_rms_de_y_vs_y");
185 
186  p2_de_x_vs_x_y->Write("p2_de_x_vs_x_y");
187  buildRMSHistogram(p2_de_x_vs_x_y)->Write("h2_rms_de_x_vs_x_y");
188  p2_de_y_vs_x_y->Write("p2_de_y_vs_x_y");
189  buildRMSHistogram(p2_de_y_vs_x_y)->Write("h2_rms_de_y_vs_x_y");
190 
191  h2_de_x_vs_x->Write("h2_de_x_vs_x");
192  h2_de_y_vs_x->Write("h2_de_y_vs_x");
193 
194  h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x");
195 
196  for (const auto& rp : m_stat) {
197  TGraph* g = new TGraph();
198 
199  char buf[100];
200  sprintf(buf, "g_mean_track_mult_run_%u", rp.first);
201  g->SetName(buf);
202 
203  for (const auto& lsp : rp.second) {
204  const int idx = g->GetN();
205  const double m = (lsp.second.s1 > 0) ? double(lsp.second.sN) / lsp.second.s1 : 0.;
206  g->SetPoint(idx, lsp.first, m);
207  }
208 
209  g->Write();
210  }
211  }
212 
213  static std::unique_ptr<TH1D> buildRMSHistogram(const std::unique_ptr<TProfile>& p) {
214  std::unique_ptr<TH1D> output =
215  std::make_unique<TH1D>("", p->GetTitle(), p->GetNbinsX(), p->GetXaxis()->GetXmin(), p->GetXaxis()->GetXmax());
216 
217  for (int bi = 1; bi <= output->GetNbinsX(); ++bi)
218  output->SetBinContent(bi, p->GetBinError(bi));
219 
220  return output;
221  }
222 
223  static std::unique_ptr<TH2D> buildRMSHistogram(const std::unique_ptr<TProfile2D>& p) {
224  std::unique_ptr<TH2D> output = std::make_unique<TH2D>("",
225  p->GetTitle(),
226  p->GetNbinsX(),
227  p->GetXaxis()->GetXmin(),
228  p->GetXaxis()->GetXmax(),
229  p->GetNbinsY(),
230  p->GetYaxis()->GetXmin(),
231  p->GetYaxis()->GetXmax());
232 
233  for (int bi_x = 1; bi_x <= output->GetNbinsX(); ++bi_x)
234  for (int bi_y = 1; bi_y <= output->GetNbinsY(); ++bi_y)
235  output->SetBinContent(bi_x, bi_y, p->GetBinError(bi_x, bi_y));
236 
237  return output;
238  }
239  };
240 
241  std::map<unsigned int, ArmPlots> armPlots;
242 };
243 
244 //----------------------------------------------------------------------------------------------------
245 
247  : tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
248  x_pitch_pixels_(iConfig.getUntrackedParameter<double>("x_pitch_pixels", 150E-3)),
249  outputFile_(iConfig.getParameter<std::string>("outputFile")),
250  events_total_(0) {
251  armPlots[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
252  armPlots[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
253 
254  armPlots[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
255  armPlots[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
256 }
257 
258 //----------------------------------------------------------------------------------------------------
259 
261  // get input
263  iEvent.getByToken(tracksToken_, tracks);
264 
265  // process tracks
266  std::map<unsigned int, unsigned int> m_mult;
267 
268  for (const auto& trk : *tracks) {
269  CTPPSDetId rpId(trk.rpId());
270  unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
271  bool rpPixel = (rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
272 
273  auto& pl = rpPlots[rpDecId];
274  if (!pl.initialized)
275  pl.init(rpPixel, x_pitch_pixels_);
276 
277  pl.fillOneTrack(trk.x(), trk.y(), trk.time());
278 
279  m_mult[rpId.arm()]++;
280  }
281 
282  for (unsigned int arm = 0; arm < 2; ++arm) {
283  auto& st = armPlots[arm].m_stat[iEvent.id().run()][iEvent.id().luminosityBlock()];
284  st.s1++;
285  st.sN += m_mult[arm];
286  }
287 
288  for (const auto& t1 : *tracks) {
289  const CTPPSDetId rpId1(t1.rpId());
290 
291  for (const auto& t2 : *tracks) {
292  const CTPPSDetId rpId2(t2.rpId());
293 
294  if (rpId1.arm() != rpId2.arm())
295  continue;
296 
297  auto& ap = armPlots[rpId1.arm()];
298 
299  const unsigned int rpDecId1 = rpId1.arm() * 100 + rpId1.station() * 10 + rpId1.rp();
300  const unsigned int rpDecId2 = rpId2.arm() * 100 + rpId2.station() * 10 + rpId2.rp();
301 
302  if (rpDecId1 == ap.rpId_N && rpDecId2 == ap.rpId_F)
303  ap.fill(t1.x(), t1.y(), t2.x(), t2.y());
304  }
305  }
306 
307  for (unsigned int i = 0; i < tracks->size(); ++i) {
308  for (unsigned int j = 0; j < tracks->size(); ++j) {
309  if (i == j)
310  continue;
311 
312  const auto& tr_i = tracks->at(i);
313  const auto& tr_j = tracks->at(j);
314 
315  if (tr_i.rpId() != tr_j.rpId())
316  continue;
317 
318  CTPPSDetId rpId(tr_i.rpId());
319  unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
320 
321  auto& pl = rpPlots[rpDecId];
322 
323  pl.h2_de_x_vs_x->Fill(tr_i.x(), tr_j.x() - tr_i.x());
324  pl.h2_de_x_vs_y->Fill(tr_i.y(), tr_j.x() - tr_i.x());
325  pl.h2_de_y_vs_x->Fill(tr_i.x(), tr_j.y() - tr_i.y());
326  pl.h2_de_y_vs_y->Fill(tr_i.y(), tr_j.y() - tr_i.y());
327  }
328  }
329 
330  // update counters
331  events_total_++;
332 
333  if (m_mult[0] > 0)
334  events_per_arm_[0]++;
335  if (m_mult[1] > 0)
336  events_per_arm_[1]++;
337 }
338 
339 //----------------------------------------------------------------------------------------------------
340 
342  edm::LogInfo("CTPPSTrackDistributionPlotter")
343  << " events processed: " << events_total_ << " (" << std::scientific << double(events_total_) << ")\n"
344  << " events with tracks in sector 45: " << events_per_arm_[0] << " (" << double(events_per_arm_[0]) << ")\n"
345  << " events with tracks in sector 56: " << events_per_arm_[1] << " (" << double(events_per_arm_[1]) << ")";
346 
347  auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
348 
349  for (const auto& it : rpPlots) {
350  gDirectory = f_out->mkdir(Form("RP %u", it.first));
351  it.second.write();
352  }
353 
354  for (const auto& it : armPlots) {
355  gDirectory = f_out->mkdir(Form("arm %u", it.first));
356  it.second.write();
357  }
358 }
359 
360 //----------------------------------------------------------------------------------------------------
361 
constexpr int32_t ceil(float num)
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
std::map< unsigned int, unsigned int > events_per_arm_
static std::unique_ptr< TH2D > buildRMSHistogram(const std::unique_ptr< TProfile2D > &p)
uint32_t arm() const
Definition: CTPPSDetId.h:51
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
void analyze(const edm::Event &, const edm::EventSetup &) override
std::map< unsigned int, ArmPlots > armPlots
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
std::map< unsigned int, std::map< unsigned int, Stat > > m_stat
int iEvent
Definition: GenABIO.cc:224
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
edm::EDGetTokenT< CTPPSLocalTrackLiteCollection > tracksToken_
static std::unique_ptr< TH1D > buildRMSHistogram(const std::unique_ptr< TProfile > &p)
std::map< unsigned int, RPPlots > rpPlots
#define M_PI
Log< level::Info, false > LogInfo
void fillOneTrack(double x, double y, double time)
auto const & tracks
cannot be loose
std::vector< CTPPSLocalTrackLite > CTPPSLocalTrackLiteCollection
Collection of CTPPSLocalTrackLite objects.
Base class for CTPPS detector IDs.
Definition: CTPPSDetId.h:32
HLT enums.
CTPPSTrackDistributionPlotter(const edm::ParameterSet &)
void fill(double x_N, double y_N, double x_F, double y_F)