Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-09-23 03:14:24

0001 /****************************************************************************
0002 * Authors:
0003 *  Jan Kašpar (jan.kaspar@gmail.com)
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;  // Gaussian approximation
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         //p_time_unc_vs_x_ClCo filled in ClCo code below
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   // identify tracking-RP tracks
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   // interpolation from tracking RPs
0517   const double z_ti = -geometry.rpTranslation(tr_ti.rpId()).z();  // the minus sign fixes a bug in the diamond geometry
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   // save distance
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   // get input
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   // get conditions
0550   const auto &geometry = iSetup.getData(geometryESToken_);
0551   edm::ESHandle<PPSAssociationCuts> ppsAssociationCuts = iSetup.getHandle(ppsAssociationCutsToken_);
0552 
0553   // track plots
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   // initialise multiplicity counters
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   // make single-RP-reco plots
0607   for (const auto &proton : *hRecoProtonsSingleRP) {
0608     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId());
0609     unsigned int decRPId = rpId.arm() * 100 + rpId.station() * 10 + rpId.rp();
0610 
0611     const bool n1f1 = (armTrackCounter_N[rpId.arm()] == 1 && armTrackCounter_F[rpId.arm()] == 1);
0612 
0613     singleRPPlots_[decRPId].fill(proton, armTrackCounter[rpId.arm()], n1f1);
0614 
0615     if (proton.validFit())
0616       singleRPMultiplicity[decRPId]++;
0617   }
0618 
0619   for (const auto it : singleRPMultiplicity)
0620     singleRPPlots_[it.first].h_multiplicity->Fill(it.second);
0621 
0622   // make multi-RP-reco plots
0623   for (const auto &proton : *hRecoProtonsMultiRP) {
0624     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId());
0625     unsigned int armId = rpId.arm();
0626 
0627     const bool n1f1 = (armTrackCounter_N[armId] == 1 && armTrackCounter_F[armId] == 1);
0628 
0629     multiRPPlots_[armId].fill(proton, armTrackCounter[armId], n1f1);
0630 
0631     if (proton.validFit())
0632       multiRPMultiplicity[armId]++;
0633   }
0634 
0635   for (const auto it : multiRPMultiplicity) {
0636     const auto &pl = multiRPPlots_[it.first];
0637     pl.h_multiplicity->Fill(it.second);
0638     pl.h2_timing_tracks_vs_prot_mult->Fill(it.second, armTimingTrackCounter[it.first]);
0639   }
0640 
0641   // define "clean condition" for each arm
0642   map<unsigned int, bool> clCo;
0643   clCo[0] = (singleRPMultiplicity[rpId_45_N_] && singleRPMultiplicity[rpId_45_F_] && multiRPMultiplicity[0] == 1);
0644   clCo[1] = (singleRPMultiplicity[rpId_56_N_] && singleRPMultiplicity[rpId_56_F_] && multiRPMultiplicity[1] == 1);
0645 
0646   // plot distances between multi-RP protons and timing tracks in the same arm
0647   for (const auto &proton : *hRecoProtonsMultiRP) {
0648     if (!proton.validFit())
0649       continue;
0650 
0651     CTPPSDetId rpId_proton((*proton.contributingLocalTracks().begin())->rpId());
0652     unsigned int armId = rpId_proton.arm();
0653     const auto &pl = multiRPPlots_[armId];
0654 
0655     for (const auto &tr : *hTracks) {
0656       CTPPSDetId rpId_tr(tr.rpId());
0657       if (rpId_tr.arm() != armId)
0658         continue;
0659 
0660       const bool trackerRP =
0661           (rpId_tr.subdetId() == CTPPSDetId::sdTrackingStrip || rpId_tr.subdetId() == CTPPSDetId::sdTrackingPixel);
0662       if (trackerRP)
0663         continue;
0664 
0665       double x_tr = -1., x_ti = -1.;
0666       double de_x = 0., de_x_unc = 0.;
0667       CalculateTimingTrackingDistance(proton, tr, geometry, x_tr, x_ti, de_x, de_x_unc);
0668 
0669       const double rd = (de_x_unc > 0.) ? de_x / de_x_unc : -1E10;
0670       const auto &ac = ppsAssociationCuts->getAssociationCuts(armId);
0671       const bool match = (ac.getTiTrMin() <= fabs(rd) && fabs(rd) <= ac.getTiTrMax());
0672 
0673       pl.h_de_x_timing_vs_tracking->Fill(de_x);
0674       pl.h_de_x_rel_timing_vs_tracking->Fill(rd);
0675       pl.h_de_x_match_timing_vs_tracking->Fill(match ? 1. : 0.);
0676 
0677       if (clCo[armId] && armTimingTrackCounter[armId] == 1) {
0678         pl.h2_x_timing_vs_x_tracking_ClCo->Fill(x_tr, x_ti);
0679 
0680         pl.h_de_x_timing_vs_tracking_ClCo->Fill(de_x);
0681         pl.h_de_x_rel_timing_vs_tracking_ClCo->Fill(rd);
0682         pl.h_de_x_match_timing_vs_tracking_ClCo->Fill(match ? 1. : 0.);
0683 
0684         pl.p_time_unc_vs_x_ClCo->Fill(x_tr, proton.timeError());
0685       }
0686     }
0687   }
0688 
0689   // plot xy maps
0690   for (const auto &proton : *hRecoProtonsMultiRP) {
0691     if (!proton.validFit())
0692       continue;
0693 
0694     CTPPSDetId rpId((*proton.contributingLocalTracks().begin())->rpId());
0695     unsigned int armId = rpId.arm();
0696     const auto &pl = multiRPPlots_[armId];
0697     const auto &nTimingTracks = armTimingTrackCounter[armId];
0698 
0699     if (!clCo[armId])
0700       continue;
0701 
0702     double x_ref = 0., y_ref = 0.;
0703     if (armId == 0) {
0704       x_ref = tr_L_N->x();
0705       y_ref = tr_L_N->y();
0706     }
0707     if (armId == 1) {
0708       x_ref = tr_R_N->x();
0709       y_ref = tr_R_N->y();
0710     }
0711 
0712     if (nTimingTracks == 0)
0713       pl.h2_y_vs_x_tt0_ClCo->Fill(x_ref, y_ref);
0714     if (nTimingTracks == 1)
0715       pl.h2_y_vs_x_tt1_ClCo->Fill(x_ref, y_ref);
0716     if (nTimingTracks > 1)
0717       pl.h2_y_vs_x_ttm_ClCo->Fill(x_ref, y_ref);
0718   }
0719 
0720   // make correlation plots
0721   for (const auto &proton_m : *hRecoProtonsMultiRP) {
0722     CTPPSDetId rpId_m((*proton_m.contributingLocalTracks().begin())->rpId());
0723     unsigned int arm = rpId_m.arm();
0724 
0725     const reco::ForwardProton *p_s_N = nullptr, *p_s_F = nullptr;
0726 
0727     for (const auto &proton_s : *hRecoProtonsSingleRP) {
0728       // skip if source of single-RP reco not included in multi-RP reco
0729       const auto key_s = proton_s.contributingLocalTracks()[0].key();
0730       bool compatible = false;
0731       for (const auto &tr_m : proton_m.contributingLocalTracks()) {
0732         if (tr_m.key() == key_s) {
0733           compatible = true;
0734           break;
0735         }
0736       }
0737 
0738       if (!compatible)
0739         continue;
0740 
0741       // fill single-multi plots
0742       CTPPSDetId rpId_s((*proton_s.contributingLocalTracks().begin())->rpId());
0743       const unsigned int idx = rpId_s.arm() * 1000 + rpId_s.station() * 100 + rpId_s.rp() * 10 + rpId_s.arm();
0744       singleMultiCorrelationPlots_[idx].fill(proton_s, proton_m);
0745 
0746       // select protons for arm-correlation plots
0747       const unsigned int rpDecId_s = rpId_s.arm() * 100 + rpId_s.station() * 10 + rpId_s.rp();
0748       if (rpDecId_s == rpId_45_N_ || rpDecId_s == rpId_56_N_)
0749         p_s_N = &proton_s;
0750       if (rpDecId_s == rpId_45_F_ || rpDecId_s == rpId_56_F_)
0751         p_s_F = &proton_s;
0752     }
0753 
0754     // fill arm-correlation plots
0755     if (p_s_N && p_s_F)
0756       armCorrelationPlots_[arm].fill(*p_s_N, *p_s_F, proton_m);
0757   }
0758 }
0759 
0760 //----------------------------------------------------------------------------------------------------
0761 
0762 void CTPPSProtonReconstructionPlotter::endJob() {
0763   LogInfo("CTPPSProtonReconstructionPlotter") << "n_non_empty_events = " << n_non_empty_events_;
0764 
0765   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0766 
0767   p_x_L_diffNF_vs_x_L_N_->Write();
0768   p_x_R_diffNF_vs_x_R_N_->Write();
0769 
0770   p_y_L_diffNF_vs_y_L_N_->Write();
0771   p_y_R_diffNF_vs_y_R_N_->Write();
0772 
0773   TDirectory *d_singleRPPlots = f_out->mkdir("singleRPPlots");
0774   for (const auto &it : singleRPPlots_) {
0775     gDirectory = d_singleRPPlots->mkdir(Form("rp%u", it.first));
0776     it.second.write();
0777   }
0778 
0779   TDirectory *d_multiRPPlots = f_out->mkdir("multiRPPlots");
0780   for (const auto &it : multiRPPlots_) {
0781     gDirectory = d_multiRPPlots->mkdir(Form("arm%u", it.first));
0782     it.second.write();
0783   }
0784 
0785   TDirectory *d_singleMultiCorrelationPlots = f_out->mkdir("singleMultiCorrelationPlots");
0786   for (const auto &it : singleMultiCorrelationPlots_) {
0787     unsigned int si_rp = it.first / 10;
0788     unsigned int mu_arm = it.first % 10;
0789 
0790     gDirectory = d_singleMultiCorrelationPlots->mkdir(Form("si_rp%u_mu_arm%u", si_rp, mu_arm));
0791     it.second.write();
0792   }
0793 
0794   TDirectory *d_armCorrelationPlots = f_out->mkdir("armCorrelationPlots");
0795   for (const auto &it : armCorrelationPlots_) {
0796     gDirectory = d_armCorrelationPlots->mkdir(Form("arm%u", it.first));
0797     it.second.write();
0798   }
0799 }
0800 
0801 //----------------------------------------------------------------------------------------------------
0802 
0803 DEFINE_FWK_MODULE(CTPPSProtonReconstructionPlotter);