File indexing completed on 2024-04-06 12:31:58
0001
0002
0003
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 "CondTools/RunInfo/interface/LHCInfoCombined.h"
0019
0020 #include "TFile.h"
0021 #include "TH1D.h"
0022
0023 #include <string>
0024
0025
0026
0027 class CTPPSHepMCDistributionPlotter : public edm::one::EDAnalyzer<> {
0028 public:
0029 explicit CTPPSHepMCDistributionPlotter(const edm::ParameterSet &);
0030
0031 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0032
0033 private:
0034 void analyze(const edm::Event &, const edm::EventSetup &) override;
0035 void endJob() override;
0036
0037 const edm::EDGetTokenT<edm::HepMCProduct> tokenHepMC_;
0038
0039 const edm::ESGetToken<LHCInfo, LHCInfoRcd> lhcInfoToken_;
0040 const edm::ESGetToken<LHCInfoPerLS, LHCInfoPerLSRcd> lhcInfoPerLSToken_;
0041 const edm::ESGetToken<LHCInfoPerFill, LHCInfoPerFillRcd> lhcInfoPerFillToken_;
0042 const bool useNewLHCInfo_;
0043
0044 const std::string outputFile_;
0045
0046 std::unique_ptr<TH1D> h_vtx_x_, h_vtx_y_, h_vtx_z_, h_vtx_t_;
0047 std::unique_ptr<TH1D> h_xi_, h_th_x_, h_th_y_;
0048 };
0049
0050
0051
0052 using namespace std;
0053 using namespace edm;
0054 using namespace HepMC;
0055
0056
0057
0058 CTPPSHepMCDistributionPlotter::CTPPSHepMCDistributionPlotter(const edm::ParameterSet &iConfig)
0059 : tokenHepMC_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("tagHepMC"))),
0060
0061 lhcInfoToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoLabel")))),
0062 lhcInfoPerLSToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoPerLSLabel")))),
0063 lhcInfoPerFillToken_(esConsumes(ESInputTag("", iConfig.getParameter<std::string>("lhcInfoPerFillLabel")))),
0064 useNewLHCInfo_(iConfig.getParameter<bool>("useNewLHCInfo")),
0065
0066 outputFile_(iConfig.getParameter<string>("outputFile")),
0067
0068 h_vtx_x_(new TH1D("h_vtx_x", ";vtx_x (mm)", 100, 0., 0.)),
0069 h_vtx_y_(new TH1D("h_vtx_y", ";vtx_y (mm)", 100, 0., 0.)),
0070 h_vtx_z_(new TH1D("h_vtx_z", ";vtx_z (mm)", 100, 0., 0.)),
0071 h_vtx_t_(new TH1D("h_vtx_t", ";vtx_t (mm)", 100, 0., 0.)),
0072
0073 h_xi_(new TH1D("h_xi", ";#xi", 100, 0., 0.30)),
0074 h_th_x_(new TH1D("h_th_x", ";#theta^{*}_{x}", 100, -300E-6, +300E-6)),
0075 h_th_y_(new TH1D("h_th_y", ";#theta^{*}_{y}", 100, -300E-6, +300E-6)) {}
0076
0077
0078
0079 void CTPPSHepMCDistributionPlotter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0080 edm::ParameterSetDescription desc;
0081
0082 desc.add<std::string>("lhcInfoLabel", "")->setComment("label of the LHCInfo record");
0083 desc.add<std::string>("lhcInfoPerLSLabel", "")->setComment("label of the LHCInfoPerLS record");
0084 desc.add<std::string>("lhcInfoPerFillLabel", "")->setComment("label of the LHCInfoPerFill record");
0085 desc.add<bool>("useNewLHCInfo", false)->setComment("flag whether to use new LHCInfoPer* records or old LHCInfo");
0086
0087 desc.add<std::string>("outputFile", "")->setComment("output file");
0088
0089 descriptions.add("ctppsHepMCDistributionPlotterDefault", desc);
0090 }
0091
0092 void CTPPSHepMCDistributionPlotter::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0093
0094 LHCInfoCombined lhcInfoCombined(iSetup, lhcInfoPerLSToken_, lhcInfoPerFillToken_, lhcInfoToken_, useNewLHCInfo_);
0095
0096
0097 edm::Handle<edm::HepMCProduct> hHepMC;
0098 iEvent.getByToken(tokenHepMC_, hHepMC);
0099 HepMC::GenEvent *hepMCEvent = (HepMC::GenEvent *)hHepMC->GetEvent();
0100
0101
0102 for (HepMC::GenEvent::vertex_iterator vit = hepMCEvent->vertices_begin(); vit != hepMCEvent->vertices_end(); ++vit) {
0103 const auto pos = (*vit)->position();
0104 h_vtx_x_->Fill(pos.x());
0105 h_vtx_y_->Fill(pos.y());
0106 h_vtx_z_->Fill(pos.z());
0107 h_vtx_t_->Fill(pos.t());
0108 }
0109
0110
0111 for (auto it = hepMCEvent->particles_begin(); it != hepMCEvent->particles_end(); ++it) {
0112 const auto &part = *it;
0113
0114
0115 if (part->pdg_id() != 2212)
0116 continue;
0117
0118 if (part->status() != 1)
0119 continue;
0120
0121 if (part->is_beam())
0122 continue;
0123
0124 const auto &mom = part->momentum();
0125 const double p_nom = lhcInfoCombined.energy;
0126
0127 if (mom.rho() / p_nom < 0.7)
0128 continue;
0129
0130 const double xi_simu = (p_nom - mom.e()) / p_nom;
0131 const double th_x_simu = mom.x() / mom.rho();
0132 const double th_y_simu = mom.y() / mom.rho();
0133
0134 h_xi_->Fill(xi_simu);
0135 h_th_x_->Fill(th_x_simu);
0136 h_th_y_->Fill(th_y_simu);
0137 }
0138 }
0139
0140
0141
0142 void CTPPSHepMCDistributionPlotter::endJob() {
0143 auto f_out = std::make_unique<TFile>(outputFile_.c_str(), "recreate");
0144
0145 h_vtx_x_->Write();
0146 h_vtx_y_->Write();
0147 h_vtx_z_->Write();
0148 h_vtx_t_->Write();
0149
0150 h_xi_->Write();
0151 h_th_x_->Write();
0152 h_th_y_->Write();
0153 }
0154
0155
0156
0157 DEFINE_FWK_MODULE(CTPPSHepMCDistributionPlotter);
