CTPPS/plugins/CTPPSHepMCDistributionPlotter.cc

0001 /****************************************************************************
0002  * Authors:
0003  *   Jan Kašpar
0004  ****************************************************************************/
0005
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/EventSetup.h"
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013
0014 #include "DataFormats/CTPPSDetId/interface/CTPPSDetId.h"
0015
0016 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0017
0018 #include "CondFormats/RunInfo/interface/LHCInfo.h"
0019 #include "CondFormats/DataRecord/interface/LHCInfoRcd.h"
0020
0021 #include "TFile.h"
0022 #include "TH1D.h"
0023
0024 #include <string>
0025
0026 //----------------------------------------------------------------------------------------------------
0027
0028 class CTPPSHepMCDistributionPlotter : public edm::one::EDAnalyzer<> {
0029 public:
0030   explicit CTPPSHepMCDistributionPlotter(const edm::ParameterSet &);
0031
0032 private:
0033   void analyze(const edm::Event &, const edm::EventSetup &) override;
0034   void endJob() override;
0035
0036   edm::EDGetTokenT<edm::HepMCProduct> tokenHepMC_;
0037   edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoESToken_;
0038   std::string outputFile_;
0039
0040   std::unique_ptr<TH1D> h_vtx_x_, h_vtx_y_, h_vtx_z_, h_vtx_t_;
0041   std::unique_ptr<TH1D> h_xi_, h_th_x_, h_th_y_;
0042 };
0043
0044 //----------------------------------------------------------------------------------------------------
0045
0046 using namespace std;
0047 using namespace edm;
0048 using namespace HepMC;
0049
0050 //----------------------------------------------------------------------------------------------------
0051
0052 CTPPSHepMCDistributionPlotter::CTPPSHepMCDistributionPlotter(const edm::ParameterSet &iConfig)
0053     : tokenHepMC_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMC"))),
0054       lhcInfoESToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0055       outputFile_(iConfig.getParameter<string>("outputFile")),
0056
0057       h_vtx_x_(new TH1D("h_vtx_x", ";vtx_x   (mm)", 100, 0., 0.)),
0058       h_vtx_y_(new TH1D("h_vtx_y", ";vtx_y   (mm)", 100, 0., 0.)),
0059       h_vtx_z_(new TH1D("h_vtx_z", ";vtx_z   (mm)", 100, 0., 0.)),
0060       h_vtx_t_(new TH1D("h_vtx_t", ";vtx_t   (mm)", 100, 0., 0.)),
0061
0062       h_xi_(new TH1D("h_xi", ";#xi", 100, 0., 0.30)),
0063       h_th_x_(new TH1D("h_th_x", ";#theta^{*}_{x}", 100, -300E-6, +300E-6)),
0064       h_th_y_(new TH1D("h_th_y", ";#theta^{*}_{y}", 100, -300E-6, +300E-6)) {}
0065
0066 //----------------------------------------------------------------------------------------------------
0067
0068 void CTPPSHepMCDistributionPlotter::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0069   // get conditions
0070   const auto &lhcInfo = iSetup.getData(lhcInfoESToken_);
0071
0072   // get input
0073   edm::Handle<edm::HepMCProduct> hHepMC;
0074   iEvent.getByToken(tokenHepMC_, hHepMC);
0075   HepMC::GenEvent *hepMCEvent = (HepMC::GenEvent *)hHepMC->GetEvent();
0076
0077   // plot vertices
0078   for (HepMC::GenEvent::vertex_iterator vit = hepMCEvent->vertices_begin(); vit != hepMCEvent->vertices_end(); ++vit) {
0079     const auto pos = (*vit)->position();
0080     h_vtx_x_->Fill(pos.x());
0081     h_vtx_y_->Fill(pos.y());
0082     h_vtx_z_->Fill(pos.z());
0083     h_vtx_t_->Fill(pos.t());
0084   }
0085
0086   // extract protons
0087   for (auto it = hepMCEvent->particles_begin(); it != hepMCEvent->particles_end(); ++it) {
0088     const auto &part = *it;
0089
0090     // accept only stable non-beam protons
0091     if (part->pdg_id() != 2212)
0092       continue;
0093
0094     if (part->status() != 1)
0095       continue;
0096
0097     if (part->is_beam())
0098       continue;
0099
0100     const auto &mom = part->momentum();
0101     const double p_nom = lhcInfo.energy();
0102
0103     if (mom.rho() / p_nom < 0.7)
0104       continue;
0105
0106     const double xi_simu = (p_nom - mom.e()) / p_nom;
0107     const double th_x_simu = mom.x() / mom.rho();
0108     const double th_y_simu = mom.y() / mom.rho();
0109
0110     h_xi_->Fill(xi_simu);
0111     h_th_x_->Fill(th_x_simu);
0112     h_th_y_->Fill(th_y_simu);
0113   }
0114 }
0115
0116 //----------------------------------------------------------------------------------------------------
0117
0118 void CTPPSHepMCDistributionPlotter::endJob() {
0119   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0120
0121   h_vtx_x_->Write();
0122   h_vtx_y_->Write();
0123   h_vtx_z_->Write();
0124   h_vtx_t_->Write();
0125
0126   h_xi_->Write();
0127   h_th_x_->Write();
0128   h_th_y_->Write();
0129 }
0130
0131 //----------------------------------------------------------------------------------------------------
0132
0133 DEFINE_FWK_MODULE(CTPPSHepMCDistributionPlotter);