Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-12-21 23:30:15

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