File indexing completed on 2021-10-27 04:18:22
0001
0002
0003
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;
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
0211 const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0212
0213
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
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
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
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
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
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
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
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
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.;
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);
