CTPPS/plugins/CTPPSBeamSmearingValidator.cc

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
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011
0012 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0013
0014 #include "TFile.h"
0015 #include "TH1D.h"
0016
0017 #include <map>
0018
0019 //----------------------------------------------------------------------------------------------------
0020
0021 class CTPPSBeamSmearingValidator : public edm::one::EDAnalyzer<> {
0022 public:
0023   explicit CTPPSBeamSmearingValidator(const edm::ParameterSet &);
0024
0025   ~CTPPSBeamSmearingValidator() override {}
0026
0027 private:
0028   void analyze(const edm::Event &, const edm::EventSetup &) override;
0029
0030   void endJob() override;
0031
0032   edm::EDGetTokenT<edm::HepMCProduct> tokenBeforeSmearing_;
0033   edm::EDGetTokenT<edm::HepMCProduct> tokenAfterSmearing_;
0034
0035   std::string outputFile_;
0036
0037   std::unique_ptr<TH1D> h_de_vtx_x_, h_de_vtx_y_, h_de_vtx_z_, h_de_vtx_t_;
0038
0039   struct SectorPlots {
0040     std::unique_ptr<TH1D> h_de_th_x, h_de_th_y, h_de_p;
0041
0042     SectorPlots()
0043         : h_de_th_x(new TH1D("", ";#Delta#theta_{x}   (rad)", 100, 0., 0.)),
0044           h_de_th_y(new TH1D("", ";#Delta#theta_{y}   (rad)", 100, 0., 0.)),
0045           h_de_p(new TH1D("", ";#Deltap   (GeV)", 100, 0., 0.)) {}
0046
0047     void write() const {
0048       h_de_th_x->Write("h_de_th_x");
0049       h_de_th_y->Write("h_de_th_y");
0050       h_de_p->Write("h_de_p");
0051     }
0052   };
0053
0054   std::map<unsigned int, SectorPlots> sectorPlots_;
0055 };
0056
0057 //----------------------------------------------------------------------------------------------------
0058
0059 using namespace std;
0060 using namespace edm;
0061 using namespace HepMC;
0062
0063 //----------------------------------------------------------------------------------------------------
0064
0065 CTPPSBeamSmearingValidator::CTPPSBeamSmearingValidator(const edm::ParameterSet &iConfig)
0066     : tokenBeforeSmearing_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagBeforeSmearing"))),
0067       tokenAfterSmearing_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagAfterSmearing"))),
0068       outputFile_(iConfig.getParameter<string>("outputFile")),
0069       h_de_vtx_x_(new TH1D("h_de_vtx_x", ";#Delta vtx_{x}   (mm)", 100, 0., 0.)),
0070       h_de_vtx_y_(new TH1D("h_de_vtx_y", ";#Delta vtx_{y}   (mm)", 100, 0., 0.)),
0071       h_de_vtx_z_(new TH1D("h_de_vtx_z", ";#Delta vtx_{z}   (mm)", 100, 0., 0.)),
0072       h_de_vtx_t_(new TH1D("h_de_vtx_t", ";#Delta vtx_{t}   (mm)", 100, 0., 0.)) {}
0073
0074 //----------------------------------------------------------------------------------------------------
0075
0076 void CTPPSBeamSmearingValidator::analyze(const edm::Event &iEvent, const edm::EventSetup &) {
0077   // get input
0078   edm::Handle<edm::HepMCProduct> hBeforeSmearing;
0079   iEvent.getByToken(tokenBeforeSmearing_, hBeforeSmearing);
0080   HepMC::GenEvent *orig = (HepMC::GenEvent *)hBeforeSmearing->GetEvent();
0081
0082   edm::Handle<edm::HepMCProduct> hAfterSmearing;
0083   iEvent.getByToken(tokenAfterSmearing_, hAfterSmearing);
0084   HepMC::GenEvent *smear = (HepMC::GenEvent *)hAfterSmearing->GetEvent();
0085
0086   // vertices
0087   GenEvent::vertex_const_iterator vold, vnew;
0088   for (vold = orig->vertices_begin(), vnew = smear->vertices_begin();
0089        vold != orig->vertices_end() && vnew != smear->vertices_end();
0090        ++vold, ++vnew) {
0091     const FourVector &vo = (*vold)->position();
0092     const FourVector &vn = (*vnew)->position();
0093
0094     // HepMC gives vertex in mm
0095     h_de_vtx_x_->Fill(vn.x() - vo.x());
0096     h_de_vtx_y_->Fill(vn.y() - vo.y());
0097     h_de_vtx_z_->Fill(vn.z() - vo.z());
0098     h_de_vtx_t_->Fill(vn.t() - vo.t());
0099   }
0100
0101   // particles
0102   GenEvent::particle_const_iterator pold, pnew;
0103   for (pold = orig->particles_begin(), pnew = smear->particles_begin();
0104        pold != orig->particles_end() && pnew != smear->particles_end();
0105        ++pold, ++pnew) {
0106     FourVector o = (*pold)->momentum(), n = (*pnew)->momentum();
0107
0108     // determine direction region
0109     signed int idx = -1;
0110     const double thetaLim = 0.01;  // rad
0111     double th = o.theta();
0112
0113     if (th < thetaLim)
0114       idx = 0;
0115     if (th > (M_PI - thetaLim))
0116       idx = 1;
0117
0118     if (idx < 0)
0119       continue;
0120
0121     /*
0122         cout << "particle\n\told: [" << o.x() << ", " << o.y() << ", " << o.z() << ", " << o.t()
0123         << "]\n\tnew: [" << n.x() << ", " << n.y() << ", " << n.z() << ", " << n.t()
0124         << "]\n\tregion: " << idx << endl;
0125     */
0126
0127     // fill histograms
0128     auto &sp = sectorPlots_[idx];
0129
0130     double othx = o.x() / o.rho(), othy = o.y() / o.rho();
0131     double nthx = n.x() / n.rho(), nthy = n.y() / n.rho();
0132
0133     sp.h_de_p->Fill(n.rho() - o.rho());
0134
0135     sp.h_de_th_x->Fill(nthx - othx);
0136     sp.h_de_th_y->Fill(nthy - othy);
0137   }
0138 }
0139
0140 //----------------------------------------------------------------------------------------------------
0141
0142 void CTPPSBeamSmearingValidator::endJob() {
0143   auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0144
0145   h_de_vtx_x_->Write();
0146   h_de_vtx_y_->Write();
0147   h_de_vtx_z_->Write();
0148   h_de_vtx_t_->Write();
0149
0150   gDirectory = f_out->mkdir("sector 45");
0151   sectorPlots_[0].write();
0152
0153   gDirectory = f_out->mkdir("sector 56");
0154   sectorPlots_[1].write();
0155 }
0156
0157 //----------------------------------------------------------------------------------------------------
0158
0159 DEFINE_FWK_MODULE(CTPPSBeamSmearingValidator);