Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-10-27 04:18:22

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005 
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012 
0013 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0014 
0015 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0016 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0017 
0018 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0019 
0020 #include "DataFormats/ProtonReco/interface/ForwardProton.h"
0021 #include "DataFormats/ProtonReco/interface/ForwardProtonFwd.h"
0022 #include "DataFormats/CTPPSReco/interface/CTPPSLocalTrackLite.h"
0023 
0024 #include "TFile.h"
0025 #include "TH1D.h"
0026 #include "TH2D.h"
0027 #include "TProfile.h"
0028 #include "TGraphErrors.h"
0029 
0030 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0031 
0032 #include <map>
0033 #include <string>
0034 
0035 //----------------------------------------------------------------------------------------------------
0036 
0037 class CTPPSProtonReconstructionSimulationValidator : public edm::one::EDAnalyzer<> {
0038 public:
0039   explicit CTPPSProtonReconstructionSimulationValidator(const edm::ParameterSet &);
0040 
0041 private:
0042   void analyze(const edm::Event &, const edm::EventSetup &) override;
0043   void endJob() override;
0044 
0045   void fillPlots(unsigned int meth_idx,
0046                  unsigned int idx,
0047                  const reco::ForwardProton &rec_pr,
0048                  const HepMC::FourVector &vtx,
0049                  const HepMC::FourVector &mom,
0050                  const LHCInfo &lhcInfo);
0051 
0052   edm::EDGetTokenT<edm::HepMCProduct> tokenHepMCBeforeSmearing_;
0053   edm::EDGetTokenT<edm::HepMCProduct> tokenHepMCAfterSmearing_;
0054 
0055   edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsSingleRP_;
0056   edm::EDGetTokenT<reco::ForwardProtonCollection> tokenRecoProtonsMultiRP_;
0057 
0058   edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoESToken_;
0059 
0060   std::string outputFile_;
0061 
0062   struct PlotGroup {
0063     std::unique_ptr<TH1D> h_de_xi;
0064     std::unique_ptr<TProfile> p_de_xi_vs_xi_simu;
0065     std::unique_ptr<TH2D> h_xi_reco_vs_xi_simu;
0066 
0067     std::unique_ptr<TH1D> h_de_th_x;
0068     std::unique_ptr<TProfile> p_de_th_x_vs_xi_simu;
0069 
0070     std::unique_ptr<TH1D> h_de_th_y;
0071     std::unique_ptr<TProfile> p_de_th_y_vs_xi_simu;
0072 
0073     std::unique_ptr<TH1D> h_de_vtx_y;
0074     std::unique_ptr<TProfile> p_de_vtx_y_vs_xi_simu;
0075 
0076     std::unique_ptr<TH1D> h_de_t;
0077     std::unique_ptr<TProfile> p_de_t_vs_xi_simu, p_de_t_vs_t_simu;
0078     std::unique_ptr<TH2D> h2_de_t_vs_t_simu;
0079 
0080     PlotGroup()
0081         : h_de_xi(new TH1D("", ";#xi_{reco} - #xi_{simu}", 100, 0., 0.)),
0082           p_de_xi_vs_xi_simu(new TProfile("", ";#xi_{simu};#xi_{reco} - #xi_{simu}", 50, 0., 0.25)),
0083           h_xi_reco_vs_xi_simu(new TH2D("", ";#xi_{simu};#xi_{reco}", 100, 0., 0.30, 100, 0., 0.30)),
0084 
0085           h_de_th_x(new TH1D("", ";#theta_{x,reco} - #theta_{x,simu}", 100, 0., 0.)),
0086           p_de_th_x_vs_xi_simu(new TProfile("", ";#xi_{simu};#theta_{x,reco} - #theta_{x,simu}", 50, 0., 0.25)),
0087 
0088           h_de_th_y(new TH1D("", ";#theta_{y,reco} - #theta_{y,simu}", 100, 0., 0.)),
0089           p_de_th_y_vs_xi_simu(new TProfile("", ";#xi_{simu};#theta_{y,reco} - #theta_{y,simu}", 50, 0., 0.25)),
0090 
0091           h_de_vtx_y(new TH1D("", ";vtx_{y,reco} - vtx_{y,simu}   (mm)", 100, 0., 0.)),
0092           p_de_vtx_y_vs_xi_simu(new TProfile("", ";#xi_{simu};vtx_{y,reco} - vtx_{y,simu} (mm)", 50, 0., 0.25)),
0093 
0094           h_de_t(new TH1D("", ";t_{reco} - t_{simu}", 100, -1., +1.)),
0095           p_de_t_vs_xi_simu(new TProfile("", ";xi_{simu};t_{reco} - t_{simu}", 50, 0., 0.25)),
0096           p_de_t_vs_t_simu(new TProfile("", ";t_{simu};t_{reco} - t_{simu}", 20, 0., 5.)),
0097           h2_de_t_vs_t_simu(new TH2D("", ";t_{simu};t_{reco} - t_{simu}", 150, 0., 5., 300, -2., +2.)) {}
0098 
0099     static TGraphErrors profileToRMSGraph(TProfile *p, const char *name = "") {
0100       TGraphErrors gr_err;
0101       gr_err.SetName(name);
0102 
0103       for (int bi = 1; bi <= p->GetNbinsX(); ++bi) {
0104         double c = p->GetBinCenter(bi);
0105         double w = p->GetBinWidth(bi);
0106 
0107         double N = p->GetBinEntries(bi);
0108         double Sy = p->GetBinContent(bi) * N;
0109         double Syy = p->GetSumw2()->At(bi);
0110 
0111         double si_sq = Syy / N - Sy * Sy / N / N;
0112         double si = (si_sq >= 0.) ? sqrt(si_sq) : 0.;
0113         double si_unc_sq = si_sq / 2. / N;  // Gaussian approximation
0114         double si_unc = (si_unc_sq >= 0.) ? sqrt(si_unc_sq) : 0.;
0115 
0116         int idx = gr_err.GetN();
0117         gr_err.SetPoint(idx, c, si);
0118         gr_err.SetPointError(idx, w / 2., si_unc);
0119       }
0120 
0121       return gr_err;
0122     }
0123 
0124     void write() const {
0125       h_xi_reco_vs_xi_simu->Write("h_xi_reco_vs_xi_simu");
0126       h_de_xi->Write("h_de_xi");
0127       p_de_xi_vs_xi_simu->Write("p_de_xi_vs_xi_simu");
0128       profileToRMSGraph(p_de_xi_vs_xi_simu.get(), "g_rms_de_xi_vs_xi_simu").Write();
0129 
0130       h_de_th_x->Write("h_de_th_x");
0131       p_de_th_x_vs_xi_simu->Write("p_de_th_x_vs_xi_simu");
0132       profileToRMSGraph(p_de_th_x_vs_xi_simu.get(), "g_rms_de_th_x_vs_xi_simu").Write();
0133 
0134       h_de_th_y->Write("h_de_th_y");
0135       p_de_th_y_vs_xi_simu->Write("p_de_th_y_vs_xi_simu");
0136       profileToRMSGraph(p_de_th_y_vs_xi_simu.get(), "g_rms_de_th_y_vs_xi_simu").Write();
0137 
0138       h_de_vtx_y->Write("h_de_vtx_y");
0139       p_de_vtx_y_vs_xi_simu->Write("p_de_vtx_y_vs_xi_simu");
0140       profileToRMSGraph(p_de_vtx_y_vs_xi_simu.get(), "g_rms_de_vtx_y_vs_xi_simu").Write();
0141 
0142       h_de_t->Write("h_de_t");
0143       p_de_t_vs_xi_simu->Write("p_de_t_vs_xi_simu");
0144       profileToRMSGraph(p_de_t_vs_xi_simu.get(), "g_rms_de_t_vs_xi_simu").Write();
0145       p_de_t_vs_t_simu->Write("p_de_t_vs_t_simu");
0146       profileToRMSGraph(p_de_t_vs_t_simu.get(), "g_rms_de_t_vs_t_simu").Write();
0147       h2_de_t_vs_t_simu->Write("h2_de_t_vs_t_simu");
0148     }
0149   };
0150 
0151   std::map<unsigned int, std::map<unsigned int, PlotGroup>> plots_;
0152 
0153   struct DoubleArmPlotGroup {
0154     std::unique_ptr<TH2D> h2_t_sh_vs_vtx_t, h2_t_dh_vs_vtx_z;
0155     std::unique_ptr<TH1D> h_t_sh_minus_vtx_t, h_t_dh_minus_vtx_z;
0156 
0157     DoubleArmPlotGroup()
0158         : h2_t_sh_vs_vtx_t(new TH2D("", ";vtx_t   (mm);(t_56 + t_45)/2   (mm)", 100, -250., +250., 100, -250., +250.)),
0159           h2_t_dh_vs_vtx_z(new TH2D("", ";vtx_z   (mm);(t_56 - t_45)/2   (mm)", 100, -250., +250., 100, -250., +250.)),
0160           h_t_sh_minus_vtx_t(new TH1D("", ";(t_56 + t_45)/2 - vtx_t   (mm)", 100, -100., +100.)),
0161           h_t_dh_minus_vtx_z(new TH1D("", ";(t_56 - t_45)/2 - vtx_z   (mm)", 100, -100., +100.)) {}
0162 
0163     void fill(double time_45, double time_56, double vtx_z, double vtx_t) {
0164       const double t_sum_half = (time_56 + time_45) / 2. * CLHEP::c_light;
0165       const double t_dif_half = (time_56 - time_45) / 2. * CLHEP::c_light;
0166 
0167       h2_t_sh_vs_vtx_t->Fill(vtx_t, t_sum_half);
0168       h_t_sh_minus_vtx_t->Fill(t_sum_half - vtx_t);
0169 
0170       h2_t_dh_vs_vtx_z->Fill(vtx_z, t_dif_half);
0171       h_t_dh_minus_vtx_z->Fill(t_dif_half - vtx_z);
0172     }
0173 
0174     void write() const {
0175       h2_t_sh_vs_vtx_t->Write("h2_t_sh_vs_vtx_t");
0176       h_t_sh_minus_vtx_t->Write("h_t_sh_minus_vtx_t");
0177 
0178       h2_t_dh_vs_vtx_z->Write("h2_t_dh_vs_vtx_z");
0179       h_t_dh_minus_vtx_z->Write("h_t_dh_minus_vtx_z");
0180     }
0181   };
0182 
0183   DoubleArmPlotGroup double_arm_plots_;
0184 };
0185 
0186 //----------------------------------------------------------------------------------------------------
0187 
0188 using namespace std;
0189 using namespace edm;
0190 using namespace HepMC;
0191 
0192 //----------------------------------------------------------------------------------------------------
0193 
0194 CTPPSProtonReconstructionSimulationValidator::CTPPSProtonReconstructionSimulationValidator(
0195     const edm::ParameterSet &iConfig)
0196     : tokenHepMCBeforeSmearing_(
0197           consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMCBeforeSmearing"))),
0198       tokenHepMCAfterSmearing_(
0199           consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMCAfterSmearing"))),
0200       tokenRecoProtonsSingleRP_(
0201           consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsSingleRP"))),
0202       tokenRecoProtonsMultiRP_(
0203           consumes<reco::ForwardProtonCollection>(iConfig.getParameter<InputTag>("tagRecoProtonsMultiRP"))),
0204       lhcInfoESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0205       outputFile_(iConfig.getParameter<string>("outputFile")) {}
0206 
0207 //----------------------------------------------------------------------------------------------------
0208 
0209 void CTPPSProtonReconstructionSimulationValidator::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0210   // get conditions
0211   const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0212 
0213   // get input
0214   edm::Handle<edm::HepMCProduct> hHepMCBeforeSmearing;
0215   iEvent.getByToken(tokenHepMCBeforeSmearing_, hHepMCBeforeSmearing);
0216   HepMC::GenEvent *hepMCEventBeforeSmearing = (HepMC::GenEvent *)hHepMCBeforeSmearing->GetEvent();
0217 
0218   edm::Handle<edm::HepMCProduct> hHepMCAfterSmearing;
0219   iEvent.getByToken(tokenHepMCAfterSmearing_, hHepMCAfterSmearing);
0220   HepMC::GenEvent *hepMCEventAfterSmearing = (HepMC::GenEvent *)hHepMCAfterSmearing->GetEvent();
0221 
0222   Handle<reco::ForwardProtonCollection> hRecoProtonsSingleRP;
0223   iEvent.getByToken(tokenRecoProtonsSingleRP_, hRecoProtonsSingleRP);
0224 
0225   Handle<reco::ForwardProtonCollection> hRecoProtonsMultiRP;
0226   iEvent.getByToken(tokenRecoProtonsMultiRP_, hRecoProtonsMultiRP);
0227 
0228   // extract vertex position
0229   bool vertex_set = false;
0230   FourVector vtx;
0231   for (auto it = hepMCEventAfterSmearing->vertices_begin(); it != hepMCEventAfterSmearing->vertices_end(); ++it) {
0232     if (vertex_set) {
0233       LogError("CTPPSProtonReconstructionSimulationValidator") << "Multiple vertices found.";
0234       return;
0235     }
0236 
0237     vertex_set = true;
0238     vtx = (*it)->position();
0239   }
0240 
0241   // extract forward protons
0242   bool proton_45_set = false;
0243   bool proton_56_set = false;
0244   FourVector mom_45, mom_56;
0245 
0246   for (auto it = hepMCEventBeforeSmearing->particles_begin(); it != hepMCEventBeforeSmearing->particles_end(); ++it) {
0247     const auto &part = *it;
0248 
0249     // accept only stable non-beam protons
0250     if (part->pdg_id() != 2212)
0251       continue;
0252 
0253     if (part->status() != 1)
0254       continue;
0255 
0256     if (part->is_beam())
0257       continue;
0258 
0259     const auto &mom = part->momentum();
0260 
0261     if (mom.e() < 4500.)
0262       continue;
0263 
0264     if (mom.z() > 0) {
0265       // 45
0266       if (proton_45_set) {
0267         LogError("CTPPSProtonReconstructionSimulationValidator") << "Found multiple protons in sector 45.";
0268         return;
0269       }
0270 
0271       proton_45_set = true;
0272       mom_45 = mom;
0273     } else {
0274       // 56
0275       if (proton_56_set) {
0276         LogError("CTPPSProtonReconstructionSimulationValidator") << "Found multiple protons in sector 56.";
0277         return;
0278       }
0279 
0280       proton_56_set = true;
0281       mom_56 = mom;
0282     }
0283   }
0284 
0285   // single-arm comparison
0286   for (const auto &handle : {hRecoProtonsSingleRP, hRecoProtonsMultiRP}) {
0287     for (const auto &rec_pr : *handle) {
0288       if (!rec_pr.validFit())
0289         continue;
0290 
0291       unsigned int idx = 12345;
0292 
0293       bool mom_set = false;
0294       FourVector mom;
0295 
0296       if (rec_pr.lhcSector() == reco::ForwardProton::LHCSector::sector45) {
0297         idx = 0;
0298         mom_set = proton_45_set;
0299         mom = mom_45;
0300       }
0301       if (rec_pr.lhcSector() == reco::ForwardProton::LHCSector::sector56) {
0302         idx = 1;
0303         mom_set = proton_56_set;
0304         mom = mom_56;
0305       }
0306 
0307       if (!mom_set)
0308         continue;
0309 
0310       unsigned int meth_idx = 1234;
0311 
0312       if (rec_pr.method() == reco::ForwardProton::ReconstructionMethod::singleRP) {
0313         meth_idx = 0;
0314 
0315         CTPPSDetId rpId((*rec_pr.contributingLocalTracks().begin())->rpId());
0316         idx = 100 * rpId.arm() + 10 * rpId.station() + rpId.rp();
0317       }
0318 
0319       if (rec_pr.method() == reco::ForwardProton::ReconstructionMethod::multiRP)
0320         meth_idx = 1;
0321 
0322       fillPlots(meth_idx, idx, rec_pr, vtx, mom, lhcInfo);
0323     }
0324   }
0325 
0326   // double-arm comparison
0327   bool time_set_45 = false, time_set_56 = false;
0328   double time_45 = 0., time_56 = 0.;
0329   for (const auto &rec_pr : *hRecoProtonsMultiRP) {
0330     if (!rec_pr.validFit())
0331       continue;
0332 
0333     // skip protons with no timing information
0334     if (rec_pr.timeError() < 1E-3)
0335       continue;
0336 
0337     if (rec_pr.lhcSector() == reco::ForwardProton::LHCSector::sector45) {
0338       time_set_45 = true;
0339       time_45 = rec_pr.time();
0340     }
0341     if (rec_pr.lhcSector() == reco::ForwardProton::LHCSector::sector56) {
0342       time_set_56 = true;
0343       time_56 = rec_pr.time();
0344     }
0345   }
0346 
0347   if (time_set_45 && time_set_56)
0348     double_arm_plots_.fill(time_45, time_56, vtx.z(), vtx.t());
0349 }
0350 
0351 //----------------------------------------------------------------------------------------------------
0352 
0353 void CTPPSProtonReconstructionSimulationValidator::fillPlots(unsigned int meth_idx,
0354                                                              unsigned int idx,
0355                                                              const reco::ForwardProton &rec_pr,
0356                                                              const HepMC::FourVector &vtx,
0357                                                              const HepMC::FourVector &mom,
0358                                                              const LHCInfo &lhcInfo) {
0359   const double p_nom = lhcInfo.energy();
0360   const double xi_simu = (p_nom - mom.rho()) / p_nom;
0361   const double th_x_simu = mom.x() / mom.rho();
0362   const double th_y_simu = mom.y() / mom.rho();
0363   const double vtx_y_simu = vtx.y();
0364   const double th_simu = sqrt(th_x_simu * th_x_simu + th_y_simu * th_y_simu);
0365   const double t_simu = -reco::ForwardProton::calculateT(p_nom, mom.rho(), th_simu);
0366 
0367   const double xi_reco = rec_pr.xi();
0368   const double th_x_reco = rec_pr.thetaX();
0369   const double th_y_reco = rec_pr.thetaY();
0370   const double vtx_y_reco = rec_pr.vertex().y() * 10.;  // conversion: cm --> mm
0371   const double t_reco = -rec_pr.t();
0372 
0373   auto &plt = plots_[meth_idx][idx];
0374 
0375   plt.h_xi_reco_vs_xi_simu->Fill(xi_simu, xi_reco);
0376   plt.h_de_xi->Fill(xi_reco - xi_simu);
0377   plt.p_de_xi_vs_xi_simu->Fill(xi_simu, xi_reco - xi_simu);
0378 
0379   plt.h_de_th_x->Fill(th_x_reco - th_x_simu);
0380   plt.p_de_th_x_vs_xi_simu->Fill(xi_simu, th_x_reco - th_x_simu);
0381 
0382   plt.h_de_th_y->Fill(th_y_reco - th_y_simu);
0383   plt.p_de_th_y_vs_xi_simu->Fill(xi_simu, th_y_reco - th_y_simu);
0384 
0385   plt.h_de_vtx_y->Fill(vtx_y_reco - vtx_y_simu);
0386   plt.p_de_vtx_y_vs_xi_simu->Fill(xi_simu, vtx_y_reco - vtx_y_simu);
0387 
0388   plt.h_de_t->Fill(t_reco - t_simu);
0389   plt.p_de_t_vs_xi_simu->Fill(xi_simu, t_reco - t_simu);
0390   plt.p_de_t_vs_t_simu->Fill(t_simu, t_reco - t_simu);
0391   plt.h2_de_t_vs_t_simu->Fill(t_simu, t_reco - t_simu);
0392 }
0393 
0394 //----------------------------------------------------------------------------------------------------
0395 
0396 void CTPPSProtonReconstructionSimulationValidator::endJob() {
0397   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0398 
0399   for (const auto &mit : plots_) {
0400     const char *method = (mit.first == 0) ? "single rp" : "multi rp";
0401     TDirectory *d_method = f_out->mkdir(method);
0402 
0403     for (const auto &eit : mit.second) {
0404       gDirectory = d_method->mkdir(Form("%i", eit.first));
0405       eit.second.write();
0406     }
0407   }
0408 
0409   gDirectory = f_out->mkdir("double arm");
0410   double_arm_plots_.write();
0411 }
0412 
0413 //----------------------------------------------------------------------------------------------------
0414 
0415 DEFINE_FWK_MODULE(CTPPSProtonReconstructionSimulationValidator);