File indexing completed on 2024-09-24 22:51:44
0001
0002
0003
0004
0005
0006 #include <memory>
0007
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016
0017 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0018
0019 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0020 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0021 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0022
0023 #include "Geometry/Records/interface/VeryForwardRealGeometryRecord.h"
0024 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0025 #include "CondFormats/DataRecord/interface/PPSAssociationCutsRcd.h"
0026 #include "CondFormats/PPSObjects/interface/PPSAssociationCuts.h"
0027
0028 #include "TFile.h"
0029 #include "TGraphErrors.h"
0030 #include "TH1D.h"
0031 #include "TH2D.h"
0032 #include "TProfile.h"
0033
0034
0035
0036 class CTPPSProtonReconstructionPlotter : public edm::one::EDAnalyzer<> {
0037 public:
0038 explicit CTPPSProtonReconstructionPlotter(const edm::ParameterSet &);
0039 ~CTPPSProtonReconstructionPlotter() override {}
0040
0041 private:
0042 void analyze(const edm::Event &, const edm::EventSetup &) override;
0043
0044 void endJob() override;
0045
0046 edm::EDGetTokenT<CTPPSLocalTrackLiteCollection> tokenTracks_;
0047 edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsSingleRP_;
0048 edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsMultiRP_;
0049 edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> geometryESToken_;
0050 edm::ESGetToken<PPSAssociationCuts, PPSAssociationCutsRcd> ppsAssociationCutsToken_;
0051
0052 unsigned int rpId_45_N_, rpId_45_F_;
0053 unsigned int rpId_56_N_, rpId_56_F_;
0054
0055 std::string outputFile_;
0056
0057 signed int maxNonEmptyEvents_;
0058
0059 static void profileToRMSGraph(TProfile *p, TGraphErrors *g) {
0060 for (int bi = 1; bi <= p->GetNbinsX(); ++bi) {
0061 double c = p->GetBinCenter(bi);
0062
0063 double N = p->GetBinEntries(bi);
0064 double Sy = p->GetBinContent(bi) * N;
0065 double Syy = p->GetSumw2()->At(bi);
0066
0067 double si_sq = Syy / N - Sy * Sy / N / N;
0068 double si = (si_sq >= 0.) ? sqrt(si_sq) : 0.;
0069 double si_unc_sq = si_sq / 2. / N;
0070 double si_unc = (si_unc_sq >= 0.) ? sqrt(si_unc_sq) : 0.;
0071
0072 int idx = g->GetN();
0073 g->SetPoint(idx, c, si);
0074 g->SetPointError(idx, 0., si_unc);
0075 }
0076 }
0077
0078 void CalculateTimingTrackingDistance(const reco::ForwardProton &proton,
0079 const CTPPSLocalTrackLite &tr,
0080 const CTPPSGeometry &geometry,
0081 double &x_tr,
0082 double &x_ti,
0083 double &de_x,
0084 double &de_x_unc);
0085
0086 struct SingleRPPlots {
0087 std::unique_ptr<TH1D> h_multiplicity;
0088 std::unique_ptr<TH1D> h_xi;
0089 std::unique_ptr<TH2D> h2_th_y_vs_xi;
0090 std::unique_ptr<TProfile> p_th_y_vs_xi;
0091
0092 std::map<unsigned int, TH1D *> m_h_xi_nTracks;
0093 std::unique_ptr<TH1D> h_xi_n1f1;
0094
0095 SingleRPPlots()
0096 : h_multiplicity(new TH1D("", ";reconstructed protons", 11, -0.5, 10.5)),
0097 h_xi(new TH1D("", ";#xi", 100, 0., 0.3)),
0098 h2_th_y_vs_xi(new TH2D("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)),
0099 p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3)),
0100 h_xi_n1f1(new TH1D("", ";#xi", 100, 0., 0.3)) {
0101 for (unsigned int n = 2; n <= 10; ++n)
0102 m_h_xi_nTracks[n] = new TH1D(*h_xi);
0103 }
0104
0105 void fill(const reco::ForwardProton &p, unsigned int nTracks, bool n1f1) {
0106 if (p.validFit()) {
0107 h_xi->Fill(p.xi());
0108
0109 const double th_y = p.thetaY();
0110 h2_th_y_vs_xi->Fill(p.xi(), th_y);
0111 p_th_y_vs_xi->Fill(p.xi(), th_y);
0112
0113 auto it = m_h_xi_nTracks.find(nTracks);
0114 if (it != m_h_xi_nTracks.end())
0115 it->second->Fill(p.xi());
0116
0117 if (n1f1)
0118 h_xi_n1f1->Fill(p.xi());
0119 }
0120 }
0121
0122 void write() const {
0123 h_multiplicity->Write("h_multiplicity");
0124 h_xi->Write("h_xi");
0125
0126 h2_th_y_vs_xi->Write("h2_th_y_vs_xi");
0127 p_th_y_vs_xi->Write("p_th_y_vs_xi");
0128
0129 TDirectory *d_top = gDirectory;
0130
0131 gDirectory = d_top->mkdir("h_xi_nTracks");
0132 for (const auto p : m_h_xi_nTracks) {
0133 char buf[100];
0134 sprintf(buf, "h_xi_nTracks_%u", p.first);
0135 p.second->Write(buf);
0136 }
0137
0138 gDirectory = d_top;
0139
0140 h_xi_n1f1->Write("h_xi_n1f1");
0141 }
0142 };
0143
0144 std::map<unsigned int, SingleRPPlots> singleRPPlots_;
0145
0146 struct MultiRPPlots {
0147 std::unique_ptr<TH1D> h_multiplicity;
0148 std::unique_ptr<TH1D> h_xi, h_th_x, h_th_y, h_vtx_y, h_t_unif, h_t, h_chi_sq, h_log_chi_sq, h_chi_sq_norm;
0149 std::unique_ptr<TH1D> h_t_xi_range1, h_t_xi_range2, h_t_xi_range3;
0150 std::unique_ptr<TH1D> h_time, h_time_unc;
0151 std::unique_ptr<TProfile> p_time_unc_vs_x_ClCo, p_time_unc_vs_xi;
0152 std::unique_ptr<TH1D> h_n_contrib_tracking_tracks, h_n_contrib_timing_tracks;
0153 std::unique_ptr<TH2D> h2_th_x_vs_xi, h2_th_y_vs_xi, h2_vtx_y_vs_xi, h2_t_vs_xi;
0154 std::unique_ptr<TProfile> p_th_x_vs_xi, p_th_y_vs_xi, p_vtx_y_vs_xi;
0155
0156 std::unique_ptr<TH2D> h2_timing_tracks_vs_prot_mult;
0157
0158 std::map<unsigned int, TH1D *> m_h_xi_nTracks;
0159 std::unique_ptr<TH1D> h_xi_n1f1;
0160
0161 std::unique_ptr<TH2D> h2_x_timing_vs_x_tracking_ClCo;
0162
0163 std::unique_ptr<TH1D> h_de_x_timing_vs_tracking, h_de_x_rel_timing_vs_tracking, h_de_x_match_timing_vs_tracking;
0164 std::unique_ptr<TH1D> h_de_x_timing_vs_tracking_ClCo, h_de_x_rel_timing_vs_tracking_ClCo,
0165 h_de_x_match_timing_vs_tracking_ClCo;
0166
0167 std::unique_ptr<TH2D> h2_y_vs_x_tt0_ClCo, h2_y_vs_x_tt1_ClCo, h2_y_vs_x_ttm_ClCo;
0168
0169 MultiRPPlots()
0170 : h_multiplicity(new TH1D("", ";reconstructed protons per event", 11, -0.5, 10.5)),
0171 h_xi(new TH1D("", ";#xi", 100, 0., 0.3)),
0172 h_th_x(new TH1D("", ";#theta_{x} (rad)", 250, -500E-6, +500E-6)),
0173 h_th_y(new TH1D("", ";#theta_{y} (rad)", 500, -1000E-6, +1000E-6)),
0174 h_vtx_y(new TH1D("", ";vtx_{y} (cm)", 100, -100E-3, +100E-3)),
0175 h_chi_sq(new TH1D("", ";#chi^{2}", 100, 0., 10.)),
0176 h_log_chi_sq(new TH1D("", ";log_{10} #chi^{2}", 100, -20., 5.)),
0177 h_chi_sq_norm(new TH1D("", ";#chi^{2}/ndf", 100, 0., 5.)),
0178 h_time(new TH1D("", ";time (ns)", 100, -2., +2.)),
0179 h_time_unc(new TH1D("", ";time unc (ns)", 100, -1., +1.)),
0180 p_time_unc_vs_x_ClCo(new TProfile("", ";x_tracking (mm);time unc (ns)", 100, 0., 30.)),
0181 p_time_unc_vs_xi(new TProfile("", ";xi;time unc (ns)", 100, 0., 0.3)),
0182 h_n_contrib_tracking_tracks(new TH1D("", ";n of contrib. tracking tracks per reco proton", 4, -0.5, +3.5)),
0183 h_n_contrib_timing_tracks(new TH1D("", ";n of contrib. timing tracks per reco proton", 4, -0.5, +3.5)),
0184 h2_th_x_vs_xi(new TH2D("", ";#xi;#theta_{x} (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)),
0185 h2_th_y_vs_xi(new TH2D("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3, 100, -500E-6, +500E-6)),
0186 h2_vtx_y_vs_xi(new TH2D("", ";#xi;vtx_{y} (cm)", 100, 0., 0.3, 100, -100E-3, +100E-3)),
0187 p_th_x_vs_xi(new TProfile("", ";#xi;#theta_{x} (rad)", 100, 0., 0.3)),
0188 p_th_y_vs_xi(new TProfile("", ";#xi;#theta_{y} (rad)", 100, 0., 0.3)),
0189 p_vtx_y_vs_xi(new TProfile("", ";#xi;vtx_{y} (cm)", 100, 0., 0.3)),
0190 h2_timing_tracks_vs_prot_mult(
0191 new TH2D("", ";reco protons per event;timing tracks per event", 11, -0.5, 10.5, 11, -0.5, 10.5)),
0192 h_xi_n1f1(new TH1D("", ";#xi", 100, 0., 0.3)),
0193
0194 h2_x_timing_vs_x_tracking_ClCo(
0195 new TH2D("", ";x_tracking (mm);x_timing (mm)", 100, 0., 20., 100, 0., 20.)),
0196 h_de_x_timing_vs_tracking(new TH1D("", ";#Delta x (mm)", 200, -1., +1.)),
0197 h_de_x_rel_timing_vs_tracking(new TH1D("", ";#Delta x / #sigma(x)", 200, -20., +20.)),
0198 h_de_x_match_timing_vs_tracking(new TH1D("", ";match between tracking and timing tracks", 2, -0.5, +1.5)),
0199 h_de_x_timing_vs_tracking_ClCo(new TH1D("", ";#Delta x (mm)", 200, -1., +1.)),
0200 h_de_x_rel_timing_vs_tracking_ClCo(new TH1D("", ";#Delta x / #sigma(x)", 200, -20., +20.)),
0201 h_de_x_match_timing_vs_tracking_ClCo(
0202 new TH1D("", ";match between tracking and timing tracks", 2, -0.5, +1.5)),
0203
0204 h2_y_vs_x_tt0_ClCo(new TH2D("", ";x (mm);y (mm)", 100, -5., 25., 100, -15., +15.)),
0205 h2_y_vs_x_tt1_ClCo(new TH2D("", ";x (mm);y (mm)", 100, -5., 25., 100, -15., +15.)),
0206 h2_y_vs_x_ttm_ClCo(new TH2D("", ";x (mm);y (mm)", 100, -5., 25., 100, -15., +15.)) {
0207 std::vector<double> v_t_bin_edges;
0208 for (double t = 0; t <= 5.;) {
0209 v_t_bin_edges.push_back(t);
0210 const double de_t = 0.05 + 0.09 * t + 0.02 * t * t;
0211 t += de_t;
0212 }
0213 h_t_unif = std::make_unique<TH1D>("", ";|t| (GeV^2)", 100, 0., 5.);
0214 h_t = std::make_unique<TH1D>("", ";|t| (GeV^2)", v_t_bin_edges.size() - 1, v_t_bin_edges.data());
0215 h_t_xi_range1 = std::make_unique<TH1D>("", ";|t| (GeV^2)", v_t_bin_edges.size() - 1, v_t_bin_edges.data());
0216 h_t_xi_range2 = std::make_unique<TH1D>("", ";|t| (GeV^2)", v_t_bin_edges.size() - 1, v_t_bin_edges.data());
0217 h_t_xi_range3 = std::make_unique<TH1D>("", ";|t| (GeV^2)", v_t_bin_edges.size() - 1, v_t_bin_edges.data());
0218 h2_t_vs_xi = std::make_unique<TH2D>(
0219 "", ";#xi;|t| (GeV^2)", 100, 0., 0.3, v_t_bin_edges.size() - 1, v_t_bin_edges.data());
0220
0221 for (unsigned int n = 2; n <= 10; ++n)
0222 m_h_xi_nTracks[n] = new TH1D(*h_xi);
0223 }
0224
0225 void fill(const reco::ForwardProton &p, unsigned int nTracks, bool n1f1) {
0226 if (!p.validFit())
0227 return;
0228
0229 unsigned int n_contrib_tracking_tracks = 0, n_contrib_timing_tracks = 0;
0230 for (const auto &tr : p.contributingLocalTracks()) {
0231 CTPPSDetId detId(tr->rpId());
0232 if (detId.subdetId() == CTPPSDetId::sdTrackingStrip || detId.subdetId() == CTPPSDetId::sdTrackingPixel)
0233 n_contrib_tracking_tracks++;
0234 if (detId.subdetId() == CTPPSDetId::sdTimingDiamond || detId.subdetId() == CTPPSDetId::sdTimingFastSilicon)
0235 n_contrib_timing_tracks++;
0236 }
0237
0238 const double th_x = p.thetaX();
0239 const double th_y = p.thetaY();
0240 const double mt = -p.t();
0241
0242 h_chi_sq->Fill(p.chi2());
0243 h_log_chi_sq->Fill(log10(p.chi2()));
0244 if (p.ndof() > 0)
0245 h_chi_sq_norm->Fill(p.normalizedChi2());
0246
0247 h_n_contrib_tracking_tracks->Fill(n_contrib_tracking_tracks);
0248 h_n_contrib_timing_tracks->Fill(n_contrib_timing_tracks);
0249
0250 h_xi->Fill(p.xi());
0251
0252 h_th_x->Fill(th_x);
0253 h_th_y->Fill(th_y);
0254
0255 h_vtx_y->Fill(p.vertex().y());
0256
0257 h_t_unif->Fill(mt);
0258 h_t->Fill(mt);
0259 if (p.xi() > 0.04 && p.xi() < 0.07)
0260 h_t_xi_range1->Fill(mt);
0261 if (p.xi() > 0.07 && p.xi() < 0.10)
0262 h_t_xi_range2->Fill(mt);
0263 if (p.xi() > 0.10 && p.xi() < 0.13)
0264 h_t_xi_range3->Fill(mt);
0265
0266 if (p.timeError() > 0.) {
0267 h_time->Fill(p.time());
0268 h_time_unc->Fill(p.timeError());
0269
0270 p_time_unc_vs_xi->Fill(p.xi(), p.timeError());
0271 }
0272
0273 h2_th_x_vs_xi->Fill(p.xi(), th_x);
0274 h2_th_y_vs_xi->Fill(p.xi(), th_y);
0275 h2_vtx_y_vs_xi->Fill(p.xi(), p.vertex().y());
0276 h2_t_vs_xi->Fill(p.xi(), mt);
0277
0278 p_th_x_vs_xi->Fill(p.xi(), th_x);
0279 p_th_y_vs_xi->Fill(p.xi(), th_y);
0280 p_vtx_y_vs_xi->Fill(p.xi(), p.vertex().y());
0281
0282 auto it = m_h_xi_nTracks.find(nTracks);
0283 if (it != m_h_xi_nTracks.end())
0284 it->second->Fill(p.xi());
0285
0286 if (n1f1)
0287 h_xi_n1f1->Fill(p.xi());
0288 }
0289
0290 void write() const {
0291 h_multiplicity->Write("h_multiplicity");
0292
0293 h_chi_sq->Write("h_chi_sq");
0294 h_log_chi_sq->Write("h_log_chi_sq");
0295 h_chi_sq_norm->Write("h_chi_sq_norm");
0296
0297 h_n_contrib_tracking_tracks->Write("h_n_contrib_tracking_tracks");
0298 h_n_contrib_timing_tracks->Write("h_n_contrib_timing_tracks");
0299
0300 h2_timing_tracks_vs_prot_mult->Write("h2_timing_tracks_vs_prot_mult");
0301
0302 h_xi->Write("h_xi");
0303
0304 h_th_x->Write("h_th_x");
0305 h2_th_x_vs_xi->Write("h2_th_x_vs_xi");
0306 p_th_x_vs_xi->Write("p_th_x_vs_xi");
0307 auto g_th_x_RMS_vs_xi = std::make_unique<TGraphErrors>();
0308 profileToRMSGraph(p_th_x_vs_xi.get(), g_th_x_RMS_vs_xi.get());
0309 g_th_x_RMS_vs_xi->Write("g_th_x_RMS_vs_xi");
0310
0311 h_th_y->Write("h_th_y");
0312 h2_th_y_vs_xi->Write("h2_th_y_vs_xi");
0313 p_th_y_vs_xi->Write("p_th_y_vs_xi");
0314 auto g_th_y_RMS_vs_xi = std::make_unique<TGraphErrors>();
0315 profileToRMSGraph(p_th_y_vs_xi.get(), g_th_y_RMS_vs_xi.get());
0316 g_th_y_RMS_vs_xi->Write("g_th_y_RMS_vs_xi");
0317
0318 h_vtx_y->Write("h_vtx_y");
0319 h2_vtx_y_vs_xi->Write("h2_vtx_y_vs_xi");
0320 p_vtx_y_vs_xi->Write("p_vtx_y_vs_xi");
0321 auto g_vtx_y_RMS_vs_xi = std::make_unique<TGraphErrors>();
0322 profileToRMSGraph(p_vtx_y_vs_xi.get(), g_vtx_y_RMS_vs_xi.get());
0323 g_vtx_y_RMS_vs_xi->Write("g_vtx_y_RMS_vs_xi");
0324
0325 h_t->Scale(1., "width");
0326
0327 h_t_unif->Write("h_t_unif");
0328 h_t->Write("h_t");
0329 h_t_xi_range1->Write("h_t_xi_range1");
0330 h_t_xi_range2->Write("h_t_xi_range2");
0331 h_t_xi_range3->Write("h_t_xi_range3");
0332
0333 h2_t_vs_xi->Write("h2_t_vs_xi");
0334
0335 h_time->Write("h_time");
0336 h_time_unc->Write("h_time_unc");
0337 p_time_unc_vs_x_ClCo->Write("p_time_unc_vs_x_ClCo");
0338 p_time_unc_vs_xi->Write("p_time_unc_vs_xi");
0339
0340 TDirectory *d_top = gDirectory;
0341
0342 gDirectory = d_top->mkdir("h_xi_nTracks");
0343 for (const auto p : m_h_xi_nTracks) {
0344 char buf[100];
0345 sprintf(buf, "h_xi_nTracks_%u", p.first);
0346 p.second->Write(buf);
0347 }
0348
0349 gDirectory = d_top;
0350
0351 h_xi_n1f1->Write("h_xi_n1f1");
0352
0353 h2_x_timing_vs_x_tracking_ClCo->Write("h2_x_timing_vs_x_tracking_ClCo");
0354
0355 h_de_x_timing_vs_tracking->Write("h_de_x_timing_vs_tracking");
0356 h_de_x_rel_timing_vs_tracking->Write("h_de_x_rel_timing_vs_tracking");
0357 h_de_x_match_timing_vs_tracking->Write("h_de_x_match_timing_vs_tracking");
0358
0359 h_de_x_timing_vs_tracking_ClCo->Write("h_de_x_timing_vs_tracking_ClCo");
0360 h_de_x_rel_timing_vs_tracking_ClCo->Write("h_de_x_rel_timing_vs_tracking_ClCo");
0361 h_de_x_match_timing_vs_tracking_ClCo->Write("h_de_x_match_timing_vs_tracking_ClCo");
0362
0363 h2_y_vs_x_tt0_ClCo->Write("h2_y_vs_x_tt0_ClCo");
0364 h2_y_vs_x_tt1_ClCo->Write("h2_y_vs_x_tt1_ClCo");
0365 h2_y_vs_x_ttm_ClCo->Write("h2_y_vs_x_ttm_ClCo");
0366 }
0367 };
0368
0369 std::map<unsigned int, MultiRPPlots> multiRPPlots_;
0370
0371 struct SingleMultiCorrelationPlots {
0372 std::unique_ptr<TH2D> h2_xi_mu_vs_xi_si;
0373 std::unique_ptr<TH1D> h_xi_diff_mu_si, h_xi_diff_si_mu;
0374
0375 std::unique_ptr<TH2D> h2_xi_diff_si_mu_vs_xi_mu;
0376 std::unique_ptr<TProfile> p_xi_diff_si_mu_vs_xi_mu;
0377
0378 std::unique_ptr<TH2D> h2_th_y_mu_vs_th_y_si;
0379
0380 SingleMultiCorrelationPlots()
0381 : h2_xi_mu_vs_xi_si(new TH2D("", ";#xi_{single};#xi_{multi}", 100, 0., 0.3, 100, 0., 0.3)),
0382 h_xi_diff_mu_si(new TH1D("", ";#xi_{multi} - #xi_{single}", 100, -0.1, +0.1)),
0383 h_xi_diff_si_mu(new TH1D("", ";#xi_{single} - #xi_{multi}", 100, -0.1, +0.1)),
0384 h2_xi_diff_si_mu_vs_xi_mu(
0385 new TH2D("", ";#xi_{multi};#xi_{single} - #xi_{multi}", 100, 0., 0.3, 100, -0.10, 0.10)),
0386 p_xi_diff_si_mu_vs_xi_mu(new TProfile("", ";#xi_{multi};#xi_{single} - #xi_{multi}", 100, 0., 0.3)),
0387 h2_th_y_mu_vs_th_y_si(
0388 new TH2D("", ";#theta^{*}_{y,si};#theta^{*}_{y,mu}", 100, -500E-6, +500E-6, 100, -500E-6, +500E-6)) {}
0389
0390 void fill(const reco::ForwardProton &p_single, const reco::ForwardProton &p_multi) {
0391 if (p_single.validFit() && p_multi.validFit()) {
0392 h2_xi_mu_vs_xi_si->Fill(p_single.xi(), p_multi.xi());
0393 h_xi_diff_mu_si->Fill(p_multi.xi() - p_single.xi());
0394 h_xi_diff_si_mu->Fill(p_single.xi() - p_multi.xi());
0395
0396 h2_xi_diff_si_mu_vs_xi_mu->Fill(p_multi.xi(), p_single.xi() - p_multi.xi());
0397 p_xi_diff_si_mu_vs_xi_mu->Fill(p_multi.xi(), p_single.xi() - p_multi.xi());
0398
0399 const double th_y_si = p_single.thetaY();
0400 const double th_y_mu = p_multi.thetaY();
0401
0402 h2_th_y_mu_vs_th_y_si->Fill(th_y_si, th_y_mu);
0403 }
0404 }
0405
0406 void write() const {
0407 h2_xi_mu_vs_xi_si->Write("h2_xi_mu_vs_xi_si");
0408 h_xi_diff_mu_si->Write("h_xi_diff_mu_si");
0409 h_xi_diff_si_mu->Write("h_xi_diff_si_mu");
0410
0411 h2_xi_diff_si_mu_vs_xi_mu->Write("h2_xi_diff_si_mu_vs_xi_mu");
0412 p_xi_diff_si_mu_vs_xi_mu->Write("p_xi_diff_si_mu_vs_xi_mu");
0413
0414 h2_th_y_mu_vs_th_y_si->Write("h2_th_y_mu_vs_th_y_si");
0415 }
0416 };
0417
0418 std::map<unsigned int, SingleMultiCorrelationPlots> singleMultiCorrelationPlots_;
0419
0420 struct ArmCorrelationPlots {
0421 std::unique_ptr<TH1D> h_xi_si_diffNF;
0422 std::unique_ptr<TH2D> h2_xi_si_diffNF_vs_xi_mu;
0423 std::unique_ptr<TProfile> p_xi_si_diffNF_vs_xi_mu;
0424
0425 std::unique_ptr<TH1D> h_th_y_si_diffNF;
0426 std::unique_ptr<TProfile> p_th_y_si_diffNF_vs_xi_mu;
0427
0428 ArmCorrelationPlots()
0429 : h_xi_si_diffNF(new TH1D("", ";#xi_{sF} - #xi_{sN}", 100, -0.02, +0.02)),
0430 h2_xi_si_diffNF_vs_xi_mu(new TH2D("", ";#xi_{m};#xi_{sF} - #xi_{sN}", 100, 0., 0.3, 100, -0.02, +0.02)),
0431 p_xi_si_diffNF_vs_xi_mu(new TProfile("", ";#xi_{m};#xi_{sF} - #xi_{sN}", 100, 0., 0.3)),
0432 h_th_y_si_diffNF(new TH1D("", ";#theta_{y,sF} - #theta_{y,sN}", 100, -100E-6, +100E-6)),
0433 p_th_y_si_diffNF_vs_xi_mu(new TProfile("", ";#xi_{m};#theta_{y,sF} - #theta_{y,sN}", 100, 0., 0.3)) {}
0434
0435 void fill(const reco::ForwardProton &p_s_N, const reco::ForwardProton &p_s_F, const reco::ForwardProton &p_m) {
0436 if (!p_s_N.validFit() || !p_s_F.validFit() || !p_m.validFit())
0437 return;
0438
0439 const double th_y_s_N = p_s_N.thetaY();
0440 const double th_y_s_F = p_s_F.thetaY();
0441
0442 h_xi_si_diffNF->Fill(p_s_F.xi() - p_s_N.xi());
0443 h2_xi_si_diffNF_vs_xi_mu->Fill(p_m.xi(), p_s_F.xi() - p_s_N.xi());
0444 p_xi_si_diffNF_vs_xi_mu->Fill(p_m.xi(), p_s_F.xi() - p_s_N.xi());
0445
0446 h_th_y_si_diffNF->Fill(th_y_s_F - th_y_s_N);
0447 p_th_y_si_diffNF_vs_xi_mu->Fill(p_m.xi(), th_y_s_F - th_y_s_N);
0448 }
0449
0450 void write() const {
0451 h_xi_si_diffNF->Write("h_xi_si_diffNF");
0452 h2_xi_si_diffNF_vs_xi_mu->Write("h2_xi_si_diffNF_vs_xi_mu");
0453 p_xi_si_diffNF_vs_xi_mu->Write("p_xi_si_diffNF_vs_xi_mu");
0454
0455 h_th_y_si_diffNF->Write("h_th_y_si_diffNF");
0456 p_th_y_si_diffNF_vs_xi_mu->Write("p_th_y_si_diffNF_vs_xi_mu");
0457 }
0458 };
0459
0460 std::map<unsigned int, ArmCorrelationPlots> armCorrelationPlots_;
0461
0462 std::unique_ptr<TProfile> p_x_L_diffNF_vs_x_L_N_, p_x_R_diffNF_vs_x_R_N_;
0463 std::unique_ptr<TProfile> p_y_L_diffNF_vs_y_L_N_, p_y_R_diffNF_vs_y_R_N_;
0464
0465 signed int n_non_empty_events_;
0466 };
0467
0468
0469
0470 using namespace std;
0471 using namespace edm;
0472
0473
0474
0475 CTPPSProtonReconstructionPlotter::CTPPSProtonReconstructionPlotter(const edm::ParameterSet &ps)
0476 : tokenTracks_(consumes<CTPPSLocalTrackLiteCollection>(ps.getParameter<edm::InputTag>("tagTracks"))),
0477 tokenRecoProtonsSingleRP_(
0478 consumes<reco::ForwardProtonCollection>(ps.getParameter<InputTag>("tagRecoProtonsSingleRP"))),
0479 tokenRecoProtonsMultiRP_(
0480 consumes<reco::ForwardProtonCollection>(ps.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
0481 geometryESToken_(esConsumes()),
0482 ppsAssociationCutsToken_(esConsumes<PPSAssociationCuts, PPSAssociationCutsRcd>()),
0483
0484 rpId_45_N_(ps.getParameter<unsigned int>("rpId_45_N")),
0485 rpId_45_F_(ps.getParameter<unsigned int>("rpId_45_F")),
0486 rpId_56_N_(ps.getParameter<unsigned int>("rpId_56_N")),
0487 rpId_56_F_(ps.getParameter<unsigned int>("rpId_56_F")),
0488
0489 outputFile_(ps.getParameter<string>("outputFile")),
0490 maxNonEmptyEvents_(ps.getUntrackedParameter<signed int>("maxNonEmptyEvents", -1)),
0491
0492 p_x_L_diffNF_vs_x_L_N_(new TProfile("p_x_L_diffNF_vs_x_L_N", ";x_{LN};x_{LF} - x_{LN}", 100, 0., +20.)),
0493 p_x_R_diffNF_vs_x_R_N_(new TProfile("p_x_R_diffNF_vs_x_R_N", ";x_{RN};x_{RF} - x_{RN}", 100, 0., +20.)),
0494
0495 p_y_L_diffNF_vs_y_L_N_(new TProfile("p_y_L_diffNF_vs_y_L_N", ";y_{LN};y_{LF} - y_{LN}", 100, -20., +20.)),
0496 p_y_R_diffNF_vs_y_R_N_(new TProfile("p_y_R_diffNF_vs_y_R_N", ";y_{RN};y_{RF} - y_{RN}", 100, -20., +20.)),
0497
0498 n_non_empty_events_(0) {}
0499
0500
0501
0502 void CTPPSProtonReconstructionPlotter::CalculateTimingTrackingDistance(const reco::ForwardProton &proton,
0503 const CTPPSLocalTrackLite &tr_ti,
0504 const CTPPSGeometry &geometry,
0505 double &x_tr,
0506 double &x_ti,
0507 double &de_x,
0508 double &de_x_unc) {
0509
0510 const auto &tr_i = *proton.contributingLocalTracks().at(0);
0511 const auto &tr_j = *proton.contributingLocalTracks().at(1);
0512
0513 const double z_i = geometry.rpTranslation(tr_i.rpId()).z();
0514 const double z_j = geometry.rpTranslation(tr_j.rpId()).z();
0515
0516
0517 const double z_ti = -geometry.rpTranslation(tr_ti.rpId()).z();
0518 const double f_i = (z_ti - z_j) / (z_i - z_j), f_j = (z_i - z_ti) / (z_i - z_j);
0519 const double x_inter = f_i * tr_i.x() + f_j * tr_j.x();
0520 const double x_inter_unc_sq = f_i * f_i * tr_i.xUnc() * tr_i.xUnc() + f_j * f_j * tr_j.xUnc() * tr_j.xUnc();
0521
0522
0523 x_tr = x_inter;
0524 x_ti = tr_ti.x();
0525
0526 de_x = tr_ti.x() - x_inter;
0527 de_x_unc = sqrt(tr_ti.xUnc() * tr_ti.xUnc() + x_inter_unc_sq);
0528 }
0529
0530
0531
0532 void CTPPSProtonReconstructionPlotter::analyze(const edm::Event &event, const edm::EventSetup &iSetup) {
0533
0534 edm::Handle<CTPPSLocalTrackLiteCollection> hTracks;
0535 event.getByToken(tokenTracks_, hTracks);
0536
0537 Handle<reco::ForwardProtonCollection> hRecoProtonsSingleRP;
0538 event.getByToken(tokenRecoProtonsSingleRP_, hRecoProtonsSingleRP);
0539
0540 Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
0541 event.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
0542
0543 if (!hRecoProtonsSingleRP->empty())
0544 n_non_empty_events_++;
0545
0546 if (maxNonEmptyEvents_ > 0 && n_non_empty_events_ > maxNonEmptyEvents_)
0547 throw cms::Exception("CTPPSProtonReconstructionPlotter") << "Number of non empty events reached maximum.";
0548
0549
0550 const auto &geometry = iSetup.getData(geometryESToken_);
0551 edm::ESHandle<PPSAssociationCuts> ppsAssociationCuts = iSetup.getHandle(ppsAssociationCutsToken_);
0552
0553
0554 const CTPPSLocalTrackLite *tr_L_N = nullptr;
0555 const CTPPSLocalTrackLite *tr_L_F = nullptr;
0556 const CTPPSLocalTrackLite *tr_R_N = nullptr;
0557 const CTPPSLocalTrackLite *tr_R_F = nullptr;
0558 std::map<unsigned int, unsigned int> armTrackCounter, armTimingTrackCounter;
0559 std::map<unsigned int, unsigned int> armTrackCounter_N, armTrackCounter_F;
0560
0561 for (const auto &tr : *hTracks) {
0562 CTPPSDetId rpId(tr.rpId());
0563 unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0564
0565 if (decRPId == rpId_45_N_) {
0566 tr_L_N = &tr;
0567 armTrackCounter_N[0]++;
0568 }
0569 if (decRPId == rpId_45_F_) {
0570 tr_L_F = &tr;
0571 armTrackCounter_F[0]++;
0572 }
0573 if (decRPId == rpId_56_N_) {
0574 tr_R_N = &tr;
0575 armTrackCounter_N[1]++;
0576 }
0577 if (decRPId == rpId_56_F_) {
0578 tr_R_F = &tr;
0579 armTrackCounter_F[1]++;
0580 }
0581
0582 armTrackCounter[rpId.arm()]++;
0583
0584 const bool trackerRP =
0585 (rpId.subdetId() == CTPPSDetId::sdTrackingStrip || rpId.subdetId() == CTPPSDetId::sdTrackingPixel);
0586 if (!trackerRP)
0587 armTimingTrackCounter[rpId.arm()]++;
0588 }
0589
0590 if (tr_L_N && tr_L_F) {
0591 p_x_L_diffNF_vs_x_L_N_->Fill(tr_L_N->x(), tr_L_F->x() - tr_L_N->x());
0592 p_y_L_diffNF_vs_y_L_N_->Fill(tr_L_N->y(), tr_L_F->y() - tr_L_N->y());
0593 }
0594
0595 if (tr_R_N && tr_R_F) {
0596 p_x_R_diffNF_vs_x_R_N_->Fill(tr_R_N->x(), tr_R_F->x() - tr_R_N->x());
0597 p_y_R_diffNF_vs_y_R_N_->Fill(tr_R_N->y(), tr_R_F->y() - tr_R_N->y());
0598 }
0599
0600
0601 std::map<unsigned int, unsigned int> singleRPMultiplicity, multiRPMultiplicity;
0602 singleRPMultiplicity[rpId_45_N_] = singleRPMultiplicity[rpId_45_F_] = singleRPMultiplicity[rpId_56_N_] =
0603 singleRPMultiplicity[rpId_56_F_] = 0;
0604 multiRPMultiplicity[0] = multiRPMultiplicity[1] = 0;
0605
0606
0607 for (const auto &proton : *hRecoProtonsSingleRP) {
0608
0609 const auto iter = proton.contributingLocalTracks().begin();
0610 if (iter == proton.contributingLocalTracks().end()) {
0611 continue;
0612 }
0613 const auto pcLTiter = *iter;
0614 assert(pcLTiter.isNonnull());
0615 CTPPSDetId rpId(pcLTiter->rpId());
0616 unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0617
0618 const bool n1f1 = (armTrackCounter_N[rpId.arm()] == 1 && armTrackCounter_F[rpId.arm()] == 1);
0619
0620 singleRPPlots_[decRPId].fill(proton, armTrackCounter[rpId.arm()], n1f1);
0621
0622 if (proton.validFit())
0623 singleRPMultiplicity[decRPId]++;
0624 }
0625
0626 for (const auto it : singleRPMultiplicity)
0627 singleRPPlots_[it.first].h_multiplicity->Fill(it.second);
0628
0629
0630 for (const auto &proton : *hRecoProtonsMultiRP) {
0631
0632 const auto &pcLTiter = *proton.contributingLocalTracks().begin();
0633 assert(pcLTiter.isNonnull());
0634 CTPPSDetId rpId(pcLTiter->rpId());
0635 unsigned int armId = rpId.arm();
0636
0637 const bool n1f1 = (armTrackCounter_N[armId] == 1 && armTrackCounter_F[armId] == 1);
0638
0639 multiRPPlots_[armId].fill(proton, armTrackCounter[armId], n1f1);
0640
0641 if (proton.validFit())
0642 multiRPMultiplicity[armId]++;
0643 }
0644
0645 for (const auto it : multiRPMultiplicity) {
0646 const auto &pl = multiRPPlots_[it.first];
0647 pl.h_multiplicity->Fill(it.second);
0648 pl.h2_timing_tracks_vs_prot_mult->Fill(it.second, armTimingTrackCounter[it.first]);
0649 }
0650
0651
0652 map<unsigned int, bool> clCo;
0653 clCo[0] = (singleRPMultiplicity[rpId_45_N_] && singleRPMultiplicity[rpId_45_F_] && multiRPMultiplicity[0] == 1);
0654 clCo[1] = (singleRPMultiplicity[rpId_56_N_] && singleRPMultiplicity[rpId_56_F_] && multiRPMultiplicity[1] == 1);
0655
0656
0657 for (const auto &proton : *hRecoProtonsMultiRP) {
0658 if (!proton.validFit())
0659 continue;
0660
0661 CTPPSDetId rpId_proton((*proton.contributingLocalTracks().begin())->rpId());
0662 unsigned int armId = rpId_proton.arm();
0663 const auto &pl = multiRPPlots_[armId];
0664
0665 for (const auto &tr : *hTracks) {
0666 CTPPSDetId rpId_tr(tr.rpId());
0667 if (rpId_tr.arm() != armId)
0668 continue;
0669
0670 const bool trackerRP =
0671 (rpId_tr.subdetId() == CTPPSDetId::sdTrackingStrip || rpId_tr.subdetId() == CTPPSDetId::sdTrackingPixel);
0672 if (trackerRP)
0673 continue;
0674
0675 double x_tr = -1., x_ti = -1.;
0676 double de_x = 0., de_x_unc = 0.;
0677 CalculateTimingTrackingDistance(proton, tr, geometry, x_tr, x_ti, de_x, de_x_unc);
0678
0679 const double rd = (de_x_unc > 0.) ? de_x / de_x_unc : -1E10;
0680 const auto &ac = ppsAssociationCuts->getAssociationCuts(armId);
0681 const bool match = (ac.getTiTrMin() <= fabs(rd) && fabs(rd) <= ac.getTiTrMax());
0682
0683 pl.h_de_x_timing_vs_tracking->Fill(de_x);
0684 pl.h_de_x_rel_timing_vs_tracking->Fill(rd);
0685 pl.h_de_x_match_timing_vs_tracking->Fill(match ? 1. : 0.);
0686
0687 if (clCo[armId] && armTimingTrackCounter[armId] == 1) {
0688 pl.h2_x_timing_vs_x_tracking_ClCo->Fill(x_tr, x_ti);
0689
0690 pl.h_de_x_timing_vs_tracking_ClCo->Fill(de_x);
0691 pl.h_de_x_rel_timing_vs_tracking_ClCo->Fill(rd);
0692 pl.h_de_x_match_timing_vs_tracking_ClCo->Fill(match ? 1. : 0.);
0693
0694 pl.p_time_unc_vs_x_ClCo->Fill(x_tr, proton.timeError());
0695 }
0696 }
0697 }
0698
0699
0700 for (const auto &proton : *hRecoProtonsMultiRP) {
0701 if (!proton.validFit())
0702 continue;
0703
0704 CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId());
0705 unsigned int armId = rpId.arm();
0706 const auto &pl = multiRPPlots_[armId];
0707 const auto &nTimingTracks = armTimingTrackCounter[armId];
0708
0709 if (!clCo[armId])
0710 continue;
0711
0712 double x_ref = 0., y_ref = 0.;
0713 if (armId == 0) {
0714 x_ref = tr_L_N->x();
0715 y_ref = tr_L_N->y();
0716 }
0717 if (armId == 1) {
0718 x_ref = tr_R_N->x();
0719 y_ref = tr_R_N->y();
0720 }
0721
0722 if (nTimingTracks == 0)
0723 pl.h2_y_vs_x_tt0_ClCo->Fill(x_ref, y_ref);
0724 if (nTimingTracks == 1)
0725 pl.h2_y_vs_x_tt1_ClCo->Fill(x_ref, y_ref);
0726 if (nTimingTracks > 1)
0727 pl.h2_y_vs_x_ttm_ClCo->Fill(x_ref, y_ref);
0728 }
0729
0730
0731 for (const auto &proton_m : *hRecoProtonsMultiRP) {
0732 if (proton_m.contributingLocalTracks().begin() == proton_m.contributingLocalTracks().end()) {
0733 continue;
0734 }
0735 CTPPSDetId rpId_m((*proton_m.contributingLocalTracks().begin())->rpId());
0736 unsigned int arm = rpId_m.arm();
0737
0738 const reco::ForwardProton *p_s_N = nullptr, *p_s_F = nullptr;
0739
0740 for (const auto &proton_s : *hRecoProtonsSingleRP) {
0741
0742 if (proton_s.contributingLocalTracks().empty()) {
0743 continue;
0744 }
0745 const auto key_s = proton_s.contributingLocalTracks()[0].key();
0746 bool compatible = false;
0747 for (const auto &tr_m : proton_m.contributingLocalTracks()) {
0748 if (tr_m.key() == key_s) {
0749 compatible = true;
0750 break;
0751 }
0752 }
0753
0754 if (!compatible)
0755 continue;
0756
0757
0758 CTPPSDetId rpId_s((*proton_s.contributingLocalTracks().begin())->rpId());
0759 const unsigned int idx = rpId_s.arm() * 1000 + rpId_s.station() * 100 + rpId_s.rp() * 10 + rpId_s.arm();
0760 singleMultiCorrelationPlots_[idx].fill(proton_s, proton_m);
0761
0762
0763 const unsigned int rpDecId_s = rpId_s.arm() * 100 + rpId_s.station() * 10 + rpId_s.rp();
0764 if (rpDecId_s == rpId_45_N_ || rpDecId_s == rpId_56_N_)
0765 p_s_N = &proton_s;
0766 if (rpDecId_s == rpId_45_F_ || rpDecId_s == rpId_56_F_)
0767 p_s_F = &proton_s;
0768 }
0769
0770
0771 if (p_s_N && p_s_F)
0772 armCorrelationPlots_[arm].fill(*p_s_N, *p_s_F, proton_m);
0773 }
0774 }
0775
0776
0777
0778 void CTPPSProtonReconstructionPlotter::endJob() {
0779 LogInfo("CTPPSProtonReconstructionPlotter") << "n_non_empty_events = " << n_non_empty_events_;
0780
0781 auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0782
0783 p_x_L_diffNF_vs_x_L_N_->Write();
0784 p_x_R_diffNF_vs_x_R_N_->Write();
0785
0786 p_y_L_diffNF_vs_y_L_N_->Write();
0787 p_y_R_diffNF_vs_y_R_N_->Write();
0788
0789 TDirectory *d_singleRPPlots = f_out->mkdir("singleRPPlots");
0790 for (const auto &it : singleRPPlots_) {
0791 gDirectory = d_singleRPPlots->mkdir(Form("rp%u", it.first));
0792 it.second.write();
0793 }
0794
0795 TDirectory *d_multiRPPlots = f_out->mkdir("multiRPPlots");
0796 for (const auto &it : multiRPPlots_) {
0797 gDirectory = d_multiRPPlots->mkdir(Form("arm%u", it.first));
0798 it.second.write();
0799 }
0800
0801 TDirectory *d_singleMultiCorrelationPlots = f_out->mkdir("singleMultiCorrelationPlots");
0802 for (const auto &it : singleMultiCorrelationPlots_) {
0803 unsigned int si_rp = it.first / 10;
0804 unsigned int mu_arm = it.first % 10;
0805
0806 gDirectory = d_singleMultiCorrelationPlots->mkdir(Form("si_rp%u_mu_arm%u", si_rp, mu_arm));
0807 it.second.write();
0808 }
0809
0810 TDirectory *d_armCorrelationPlots = f_out->mkdir("armCorrelationPlots");
0811 for (const auto &it : armCorrelationPlots_) {
0812 gDirectory = d_armCorrelationPlots->mkdir(Form("arm%u", it.first));
0813 it.second.write();
0814 }
0815 }
0816
0817
0818
0819 DEFINE_FWK_MODULE(CTPPSProtonReconstructionPlotter);
