Forward/test/SimG4FluxAnalyzer.cc

0001 // system include files
0002 #include <iostream>
0003 #include <fstream>
0004 #include <map>
0005 #include <sstream>
0006 #include <string>
0007 #include <vector>
0008
0009 // user include files
0010 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0011
0012 #include "FWCore/Framework/interface/Frameworkfwd.h"
0013 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022
0023 #include "SimDataFormats/CaloTest/interface/ParticleFlux.h"
0024
0025 // Root objects
0026 #include "TROOT.h"
0027 #include "TSystem.h"
0028 #include "TFile.h"
0029 #include "TTree.h"
0030
0031 //#define EDM_ML_DEBUG
0032
0033 class SimG4FluxAnalyzer : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0034 public:
0035   explicit SimG4FluxAnalyzer(edm::ParameterSet const&);
0036   ~SimG4FluxAnalyzer() override;
0037
0038   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0039
0040 private:
0041   void beginJob() override;
0042   void beginRun(edm::Run const&, edm::EventSetup const&) override {}
0043   void endRun(edm::Run const&, edm::EventSetup const&) override {}
0044   void analyze(edm::Event const&, edm::EventSetup const&) override;
0045
0046   edm::Service<TFileService> fs_;
0047   std::vector<std::string> lvNames_;
0048   std::vector<edm::EDGetTokenT<ParticleFlux>> tok_PF_;
0049   TTree* tree_;
0050
0051   std::vector<std::string> detName_;
0052   std::vector<int> detId_, pdgId_, vxType_;
0053   std::vector<float> tof_, vtxX_, vtxY_, vtxZ_;
0054   std::vector<float> hitPtX_, hitPtY_, hitPtZ_;
0055   std::vector<float> momX_, momY_, momZ_;
0056 };
0057
0058 SimG4FluxAnalyzer::SimG4FluxAnalyzer(const edm::ParameterSet& iConfig) {
0059   usesResource("TFileService");
0060
0061   //now do whatever initialization is needed
0062   lvNames_ = iConfig.getParameter<std::vector<std::string>>("LVNames");
0063 #ifdef EDM_ML_DEBUG
0064   std::ostringstream st1;
0065   for (const auto& name : lvNames_)
0066     st1 << " " << name;
0067   edm::LogVerbatim("SimG4FluxProducer") << "SimG4FluxAnalyzer:: for " << lvNames_.size() << " names:" << st1.str();
0068 #endif
0069
0070   for (const auto& name : lvNames_) {
0071     std::string tagn = name + "ParticleFlux";
0072     tok_PF_.push_back(consumes<ParticleFlux>(edm::InputTag("g4SimHits", tagn)));
0073 #ifdef EDM_ML_DEBUG
0074     edm::LogVerbatim("SimG4FluxProducer") << "Flux source " << edm::InputTag("g4SimHits", tagn);
0075 #endif
0076   }
0077 }
0078
0079 SimG4FluxAnalyzer::~SimG4FluxAnalyzer() {}
0080
0081 void SimG4FluxAnalyzer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0082   edm::ParameterSetDescription desc;
0083   desc.setUnknown();
0084   std::vector<std::string> lvnames = {
0085       "TotemT1Part1", "TotemT1Part2", "TotemT1Part3", "TotemT2Part1", "TotemT2Part2", "TotemT2Part3"};
0086   desc.add<std::vector<std::string>>("LVNames", lvnames);
0087   descriptions.add("SimG4FluxAnalyzer", desc);
0088 }
0089
0090 void SimG4FluxAnalyzer::beginJob() {
0091   tree_ = fs_->make<TTree>("Flux", "ParticleFlux");
0092   tree_->Branch("DetectorName", &detName_);
0093   tree_->Branch("DetectorID", &detId_);
0094   tree_->Branch("ParticleCode", &pdgId_);
0095   tree_->Branch("VertexType", &vxType_);
0096   tree_->Branch("ParticleTOF", &tof_);
0097   tree_->Branch("VertexX", &vtxX_);
0098   tree_->Branch("VertexY", &vtxY_);
0099   tree_->Branch("VertexZ", &vtxZ_);
0100   tree_->Branch("HitPointX", &hitPtX_);
0101   tree_->Branch("HitPointY", &hitPtY_);
0102   tree_->Branch("HitPointZ", &hitPtZ_);
0103   tree_->Branch("MomentumX", &momX_);
0104   tree_->Branch("MomentumY", &momY_);
0105   tree_->Branch("MomentumZ", &momZ_);
0106 }
0107
0108 void SimG4FluxAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0109   //Loop over all flux containers
0110   detName_.clear();
0111   detId_.clear();
0112   pdgId_.clear();
0113   vxType_.clear();
0114   tof_.clear();
0115   vtxX_.clear();
0116   vtxY_.clear();
0117   vtxZ_.clear();
0118   hitPtX_.clear();
0119   hitPtY_.clear();
0120   hitPtZ_.clear();
0121   momX_.clear();
0122   momY_.clear();
0123   momZ_.clear();
0124 #ifdef EDM_ML_DEBUG
0125   unsigned int k(0);
0126 #endif
0127   for (const auto& token : tok_PF_) {
0128     edm::Handle<ParticleFlux> pFlux;
0129     iEvent.getByToken(token, pFlux);
0130     if (pFlux.isValid()) {
0131       const ParticleFlux* pflux = pFlux.product();
0132       const std::string& name = pflux->getName();
0133       int id = pflux->getId();
0134       std::vector<ParticleFlux::flux> flux = pflux->getFlux();
0135 #ifdef EDM_ML_DEBUG
0136       edm::LogVerbatim("SimG4FluxProducer")
0137           << "SimG4FluxAnalyzer:: ParticleFlux for " << lvNames_[k] << " has " << pflux->getComponents() << " entries";
0138       ++k;
0139       unsigned k1(0);
0140 #endif
0141       for (const auto& element : flux) {
0142         detName_.push_back(name);
0143         detId_.push_back(id);
0144         pdgId_.push_back(element.pdgId);
0145         vxType_.push_back(element.vxType);
0146         tof_.push_back(element.tof);
0147         vtxX_.push_back(element.vertex.X());
0148         vtxY_.push_back(element.vertex.Y());
0149         vtxZ_.push_back(element.vertex.Z());
0150         hitPtX_.push_back(element.hitPoint.X());
0151         hitPtY_.push_back(element.hitPoint.Y());
0152         hitPtZ_.push_back(element.hitPoint.Z());
0153         momX_.push_back(element.momentum.X());
0154         momY_.push_back(element.momentum.Y());
0155         momZ_.push_back(element.momentum.Z());
0156 #ifdef EDM_ML_DEBUG
0157         edm::LogVerbatim("SimG4FluxProducer")
0158             << "Flux[" << k1 << "] PDGId " << element.pdgId << " VT " << element.vxType << " ToF " << element.tof
0159             << " Vertex " << element.vertex << " Hit " << element.hitPoint << " p " << element.momentum;
0160         ++k1;
0161 #endif
0162       }
0163     }
0164   }
0165 #ifdef EDM_ML_DEBUG
0166   edm::LogVerbatim("SimG4FluxProducer") << "All flux compnents have " << detName_.size() << " entries";
0167 #endif
0168   if (!detName_.empty())
0169     tree_->Fill();
0170 }
0171
0172 //define this as a plug-in
0173 DEFINE_FWK_MODULE(SimG4FluxAnalyzer);