File indexing completed on 2023-03-17 11:27:04
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0012
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017
0018 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0019 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0020 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLiteFwd.h"
0021
0022 #include "TFile.h"
0023 #include "TH2D.h"
0024 #include "TH1D.h"
0025 #include "TProfile.h"
0026 #include "TProfile2D.h"
0027 #include "TGraph.h"
0028
0029 #include <map>
0030 #include <memory>
0031
0032
0033
0034 class CTPPSTrackDistributionPlotter : public edm::one::EDAnalyzer<> {
0035 public:
0036 explicit CTPPSTrackDistributionPlotter(const edm::ParameterSet&);
0037
0038 ~CTPPSTrackDistributionPlotter() override {}
0039
0040 private:
0041 void analyze(const edm::Event&, const edm::EventSetup&) override;
0042 void endJob() override;
0043
0044 edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tracksToken_;
0045
0046 double x_pitch_pixels_;
0047
0048 std::string outputFile_;
0049
0050 unsigned int events_total_;
0051 std::map<unsigned int, unsigned int> events_per_arm_;
0052
0053 struct RPPlots {
0054 bool initialized;
0055
0056 std::unique_ptr<TH2D> h2_y_vs_x;
0057 std::unique_ptr<TProfile> p_y_vs_x;
0058 std::unique_ptr<TH1D> h_x;
0059 std::unique_ptr<TH1D> h_y;
0060 std::unique_ptr<TH1D> h_time;
0061
0062 std::unique_ptr<TH2D> h2_de_x_vs_x;
0063 std::unique_ptr<TH2D> h2_de_x_vs_y;
0064 std::unique_ptr<TH2D> h2_de_y_vs_x;
0065 std::unique_ptr<TH2D> h2_de_y_vs_y;
0066
0067 RPPlots() : initialized(false) {}
0068
0069 void init(bool pixel, double pitch) {
0070 const double bin_size_x = (pixel) ? pitch * cos(18.4 / 180. * M_PI) : 100E-3;
0071
0072 h2_y_vs_x = std::make_unique<TH2D>("", "", 300, -10., +70., 600, -30., +30.);
0073 p_y_vs_x = std::make_unique<TProfile>("", "", 300, -10., +70.);
0074
0075 int n_mi = ceil(10. / bin_size_x);
0076 int n_pl = ceil(70. / bin_size_x);
0077
0078 h_x = std::make_unique<TH1D>("", "", n_mi + n_pl, -n_mi * bin_size_x, +n_pl * bin_size_x);
0079
0080 h_y = std::make_unique<TH1D>("", "", 300, -15., +15.);
0081
0082 h_time = std::make_unique<TH1D>("", ";time", 500, -50., +50.);
0083
0084 h2_de_x_vs_x =
0085 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.);
0086 h2_de_x_vs_y =
0087 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.);
0088 h2_de_y_vs_x =
0089 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.);
0090 h2_de_y_vs_y =
0091 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.);
0092
0093 initialized = true;
0094 }
0095
0096 void fillOneTrack(double x, double y, double time) {
0097 h2_y_vs_x->Fill(x, y);
0098 p_y_vs_x->Fill(x, y);
0099 h_x->Fill(x);
0100 h_y->Fill(y);
0101 h_time->Fill(time);
0102 }
0103
0104 void write() const {
0105 h2_y_vs_x->Write("h2_y_vs_x");
0106 p_y_vs_x->Write("p_y_vs_x");
0107 h_x->Write("h_x");
0108 h_y->Write("h_y");
0109 h_time->Write("h_time");
0110
0111 h2_de_x_vs_x->Write("h2_de_x_vs_x");
0112 h2_de_x_vs_y->Write("h2_de_x_vs_y");
0113 h2_de_y_vs_x->Write("h2_de_y_vs_x");
0114 h2_de_y_vs_y->Write("h2_de_y_vs_y");
0115 }
0116 };
0117
0118 std::map<unsigned int, RPPlots> rpPlots;
0119
0120 struct ArmPlots {
0121 unsigned int rpId_N, rpId_F;
0122
0123 std::unique_ptr<TH1D> h_de_x, h_de_y;
0124 std::unique_ptr<TProfile> p_de_x_vs_x, p_de_y_vs_x;
0125 std::unique_ptr<TProfile> p_de_x_vs_y, p_de_y_vs_y;
0126 std::unique_ptr<TProfile2D> p2_de_x_vs_x_y, p2_de_y_vs_x_y;
0127 std::unique_ptr<TH2D> h2_de_x_vs_x, h2_de_y_vs_x;
0128 std::unique_ptr<TH2D> h2_de_y_vs_de_x;
0129
0130 struct Stat {
0131 unsigned int sN = 0, s1 = 0;
0132 };
0133
0134 std::map<unsigned int, std::map<unsigned int, Stat>> m_stat;
0135
0136 ArmPlots()
0137 : h_de_x(new TH1D("", ";x^{F} - x^{N}", 100, -1., +1.)),
0138 h_de_y(new TH1D("", ";y^{F} - y^{N}", 100, -1., +1.)),
0139
0140 p_de_x_vs_x(new TProfile("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., "s")),
0141 p_de_y_vs_x(new TProfile("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., "s")),
0142
0143 p_de_x_vs_y(new TProfile("", ";y^{N};x^{F} - x^{N}", 80, -20., +20., "s")),
0144 p_de_y_vs_y(new TProfile("", ";y^{N};y^{F} - y^{N}", 80, -20., +20., "s")),
0145
0146 p2_de_x_vs_x_y(new TProfile2D("", ";x;y;x^{F} - x^{N}", 80, 0., 40., 80, -20., +20., "s")),
0147 p2_de_y_vs_x_y(new TProfile2D("", ";x;y;y^{F} - y^{N}", 80, 0., 40., 80, -20., +20., "s")),
0148
0149 h2_de_x_vs_x(new TH2D("", ";x^{N};x^{F} - x^{N}", 80, 0., 40., 100, -1., +1.)),
0150 h2_de_y_vs_x(new TH2D("", ";x^{N};y^{F} - y^{N}", 80, 0., 40., 100, -1., +1.)),
0151 h2_de_y_vs_de_x(new TH2D("", ";x^{F} - x^{N};y^{F} - y^{N}", 100, -1., +1., 100, -1., +1.)) {}
0152
0153 void fill(double x_N, double y_N, double x_F, double y_F) {
0154 h_de_x->Fill(x_F - x_N);
0155 h_de_y->Fill(y_F - y_N);
0156
0157 p_de_x_vs_x->Fill(x_N, x_F - x_N);
0158 p_de_y_vs_x->Fill(x_N, y_F - y_N);
0159
0160 p_de_x_vs_y->Fill(y_N, x_F - x_N);
0161 p_de_y_vs_y->Fill(y_N, y_F - y_N);
0162
0163 p2_de_x_vs_x_y->Fill(x_N, y_N, x_F - x_N);
0164 p2_de_y_vs_x_y->Fill(x_N, y_N, y_F - y_N);
0165
0166 h2_de_x_vs_x->Fill(x_N, x_F - x_N);
0167 h2_de_y_vs_x->Fill(x_N, y_F - y_N);
0168
0169 h2_de_y_vs_de_x->Fill(x_F - x_N, y_F - y_N);
0170 }
0171
0172 void write() const {
0173 h_de_x->Write("h_de_x");
0174 h_de_y->Write("h_de_y");
0175
0176 p_de_x_vs_x->Write("p_de_x_vs_x");
0177 buildRMSHistogram(p_de_x_vs_x)->Write("h_rms_de_x_vs_x");
0178 p_de_y_vs_x->Write("p_de_y_vs_x");
0179 buildRMSHistogram(p_de_y_vs_x)->Write("h_rms_de_y_vs_x");
0180
0181 p_de_x_vs_y->Write("p_de_x_vs_y");
0182 buildRMSHistogram(p_de_x_vs_y)->Write("h_rms_de_x_vs_y");
0183 p_de_y_vs_y->Write("p_de_y_vs_y");
0184 buildRMSHistogram(p_de_y_vs_y)->Write("h_rms_de_y_vs_y");
0185
0186 p2_de_x_vs_x_y->Write("p2_de_x_vs_x_y");
0187 buildRMSHistogram(p2_de_x_vs_x_y)->Write("h2_rms_de_x_vs_x_y");
0188 p2_de_y_vs_x_y->Write("p2_de_y_vs_x_y");
0189 buildRMSHistogram(p2_de_y_vs_x_y)->Write("h2_rms_de_y_vs_x_y");
0190
0191 h2_de_x_vs_x->Write("h2_de_x_vs_x");
0192 h2_de_y_vs_x->Write("h2_de_y_vs_x");
0193
0194 h2_de_y_vs_de_x->Write("h2_de_y_vs_de_x");
0195
0196 for (const auto& rp : m_stat) {
0197 TGraph* g = new TGraph();
0198
0199 char buf[100];
0200 sprintf(buf, "g_mean_track_mult_run_%u", rp.first);
0201 g->SetName(buf);
0202
0203 for (const auto& lsp : rp.second) {
0204 const int idx = g->GetN();
0205 const double m = (lsp.second.s1 > 0) ? double(lsp.second.sN) / lsp.second.s1 : 0.;
0206 g->SetPoint(idx, lsp.first, m);
0207 }
0208
0209 g->Write();
0210 }
0211 }
0212
0213 static std::unique_ptr<TH1D> buildRMSHistogram(const std::unique_ptr<TProfile>& p) {
0214 std::unique_ptr<TH1D> output =
0215 std::make_unique<TH1D>("", p->GetTitle(), p->GetNbinsX(), p->GetXaxis()->GetXmin(), p->GetXaxis()->GetXmax());
0216
0217 for (int bi = 1; bi <= output->GetNbinsX(); ++bi)
0218 output->SetBinContent(bi, p->GetBinError(bi));
0219
0220 return output;
0221 }
0222
0223 static std::unique_ptr<TH2D> buildRMSHistogram(const std::unique_ptr<TProfile2D>& p) {
0224 std::unique_ptr<TH2D> output = std::make_unique<TH2D>("",
0225 p->GetTitle(),
0226 p->GetNbinsX(),
0227 p->GetXaxis()->GetXmin(),
0228 p->GetXaxis()->GetXmax(),
0229 p->GetNbinsY(),
0230 p->GetYaxis()->GetXmin(),
0231 p->GetYaxis()->GetXmax());
0232
0233 for (int bi_x = 1; bi_x <= output->GetNbinsX(); ++bi_x)
0234 for (int bi_y = 1; bi_y <= output->GetNbinsY(); ++bi_y)
0235 output->SetBinContent(bi_x, bi_y, p->GetBinError(bi_x, bi_y));
0236
0237 return output;
0238 }
0239 };
0240
0241 std::map<unsigned int, ArmPlots> armPlots;
0242 };
0243
0244
0245
0246 CTPPSTrackDistributionPlotter::CTPPSTrackDistributionPlotter(const edm::ParameterSet& iConfig)
0247 : tracksToken_(consumes<CTPPSLocalTrackLiteCollection>(iConfig.getParameter<edm::InputTag>("tagTracks"))),
0248 x_pitch_pixels_(iConfig.getUntrackedParameter<double>("x_pitch_pixels", 150E-3)),
0249 outputFile_(iConfig.getParameter<std::string>("outputFile")),
0250 events_total_(0) {
0251 armPlots[0].rpId_N = iConfig.getParameter<unsigned int>("rpId_45_N");
0252 armPlots[0].rpId_F = iConfig.getParameter<unsigned int>("rpId_45_F");
0253
0254 armPlots[1].rpId_N = iConfig.getParameter<unsigned int>("rpId_56_N");
0255 armPlots[1].rpId_F = iConfig.getParameter<unsigned int>("rpId_56_F");
0256 }
0257
0258
0259
0260 void CTPPSTrackDistributionPlotter::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0261
0262 edm::Handle<CTPPSLocalTrackLiteCollection> tracks;
0263 iEvent.getByToken(tracksToken_, tracks);
0264
0265
0266 std::map<unsigned int, unsigned int> m_mult;
0267
0268 for (const auto& trk : *tracks) {
0269 CTPPSDetId rpId(trk.rpId());
0270 unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0271 bool rpPixel = (rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
0272
0273 auto& pl = rpPlots[rpDecId];
0274 if (!pl.initialized)
0275 pl.init(rpPixel, x_pitch_pixels_);
0276
0277 pl.fillOneTrack(trk.x(), trk.y(), trk.time());
0278
0279 m_mult[rpId.arm()]++;
0280 }
0281
0282 for (unsigned int arm = 0; arm < 2; ++arm) {
0283 auto& st = armPlots[arm].m_stat[iEvent.id().run()][iEvent.id().luminosityBlock()];
0284 st.s1++;
0285 st.sN += m_mult[arm];
0286 }
0287
0288 for (const auto& t1 : *tracks) {
0289 const CTPPSDetId rpId1(t1.rpId());
0290
0291 for (const auto& t2 : *tracks) {
0292 const CTPPSDetId rpId2(t2.rpId());
0293
0294 if (rpId1.arm() != rpId2.arm())
0295 continue;
0296
0297 auto& ap = armPlots[rpId1.arm()];
0298
0299 const unsigned int rpDecId1 = rpId1.arm() * 100 + rpId1.station() * 10 + rpId1.rp();
0300 const unsigned int rpDecId2 = rpId2.arm() * 100 + rpId2.station() * 10 + rpId2.rp();
0301
0302 if (rpDecId1 == ap.rpId_N && rpDecId2 == ap.rpId_F)
0303 ap.fill(t1.x(), t1.y(), t2.x(), t2.y());
0304 }
0305 }
0306
0307 for (unsigned int i = 0; i < tracks->size(); ++i) {
0308 for (unsigned int j = 0; j < tracks->size(); ++j) {
0309 if (i == j)
0310 continue;
0311
0312 const auto& tr_i = tracks->at(i);
0313 const auto& tr_j = tracks->at(j);
0314
0315 if (tr_i.rpId() != tr_j.rpId())
0316 continue;
0317
0318 CTPPSDetId rpId(tr_i.rpId());
0319 unsigned int rpDecId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0320
0321 auto& pl = rpPlots[rpDecId];
0322
0323 pl.h2_de_x_vs_x->Fill(tr_i.x(), tr_j.x() - tr_i.x());
0324 pl.h2_de_x_vs_y->Fill(tr_i.y(), tr_j.x() - tr_i.x());
0325 pl.h2_de_y_vs_x->Fill(tr_i.x(), tr_j.y() - tr_i.y());
0326 pl.h2_de_y_vs_y->Fill(tr_i.y(), tr_j.y() - tr_i.y());
0327 }
0328 }
0329
0330
0331 events_total_++;
0332
0333 if (m_mult[0] > 0)
0334 events_per_arm_[0]++;
0335 if (m_mult[1] > 0)
0336 events_per_arm_[1]++;
0337 }
0338
0339
0340
0341 void CTPPSTrackDistributionPlotter::endJob() {
0342 edm::LogInfo("CTPPSTrackDistributionPlotter")
0343 << " events processed: " << events_total_ << " (" << std::scientific << double(events_total_) << ")\n"
0344 << " events with tracks in sector 45: " << events_per_arm_[0] << " (" << double(events_per_arm_[0]) << ")\n"
0345 << " events with tracks in sector 56: " << events_per_arm_[1] << " (" << double(events_per_arm_[1]) << ")";
0346
0347 auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0348
0349 for (const auto& it : rpPlots) {
0350 gDirectory = f_out->mkdir(Form("RP %u", it.first));
0351 it.second.write();
0352 }
0353
0354 for (const auto& it : armPlots) {
0355 gDirectory = f_out->mkdir(Form("arm %u", it.first));
0356 it.second.write();
0357 }
0358 }
0359
0360
0361
0362 DEFINE_FWK_MODULE(CTPPSTrackDistributionPlotter);
