Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-18 09:15:34

0001 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0002 #include "DataFormats/Math/interface/GeantUnits.h"
0003 
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/Framework/interface/Frameworkfwd.h"
0007 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ServiceRegistry/interface/Service.h"
0012 #include "FWCore/Utilities/interface/InputTag.h"
0013 
0014 #include "SimDataFormats/CaloHit/interface/MaterialInformation.h"
0015 
0016 #include "TProfile.h"
0017 #include "TProfile2D.h"
0018 
0019 #include <iostream>
0020 #include <string>
0021 #include <vector>
0022 
0023 //#define EDM_ML_DEBUG
0024 
0025 using namespace geant_units::operators;
0026 
0027 class MaterialBudgetVolumeAnalysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0028 public:
0029   explicit MaterialBudgetVolumeAnalysis(edm::ParameterSet const&);
0030 
0031   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0032 
0033 private:
0034   void analyze(edm::Event const&, edm::EventSetup const&) override;
0035   void bookHisto();
0036 
0037   const std::vector<std::string> names_;
0038   const edm::InputTag tag_;
0039   const int binEta_, binPhi_;
0040   const double etaLow_, etaHigh_, phiLow_, phiHigh_;
0041   edm::EDGetTokenT<edm::MaterialInformationContainer> tok_info_;
0042   std::vector<TProfile*> meStepEta_, meStepPhi_;
0043   std::vector<TProfile*> meRadLEta_, meRadLPhi_;
0044   std::vector<TProfile*> meIntLEta_, meIntLPhi_;
0045   std::vector<TProfile2D*> meStepEtaPhi_, meRadLEtaPhi_, meIntLEtaPhi_;
0046 };
0047 
0048 MaterialBudgetVolumeAnalysis::MaterialBudgetVolumeAnalysis(const edm::ParameterSet& p)
0049     : names_(p.getParameter<std::vector<std::string> >("names")),
0050       tag_(p.getParameter<edm::InputTag>("inputTag")),
0051       binEta_(p.getParameter<int>("nBinEta")),
0052       binPhi_(p.getParameter<int>("nBinPhi")),
0053       etaLow_(p.getParameter<double>("etaLow")),
0054       etaHigh_(p.getParameter<double>("etaHigh")),
0055       phiLow_(-1._pi),
0056       phiHigh_(1._pi) {
0057   usesResource(TFileService::kSharedResource);
0058   tok_info_ = consumes<edm::MaterialInformationContainer>(tag_);
0059 
0060   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolumeAnalysis: Eta plot: NX " << binEta_ << " Range "
0061                                      << -etaLow_ << ":" << etaHigh_ << " Phi plot: NX " << binPhi_ << " Range "
0062                                      << -1._pi << ":" << 1._pi << " for " << names_.size() << " detectors from "
0063                                      << tag_;
0064   std::ostringstream st1;
0065   for (unsigned int k = 0; k < names_.size(); ++k)
0066     st1 << " [" << k << "] " << names_[k];
0067   edm::LogVerbatim("MaterialBudget") << "MaterialBudgetVolume: " << st1.str();
0068   bookHisto();
0069 }
0070 
0071 void MaterialBudgetVolumeAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0072   edm::ParameterSetDescription desc;
0073   std::vector<std::string> names = {
0074       "BEAM", "BEAM1", "BEAM2", "BEAM3", "BEAM4", "Tracker", "ECAL", "HCal", "MUON", "VCAL", "MGNT", "OQUA", "CALOEC"};
0075   desc.add<std::vector<std::string> >("names", names);
0076   desc.add<edm::InputTag>("inputTag", edm::InputTag("g4SimHits", "MaterialInformation"));
0077   desc.add<int>("nBinEta", 300);
0078   desc.add<int>("nBinPhi", 180);
0079   desc.add<double>("etaLow", -6.0);
0080   desc.add<double>("etaHigh", 6.0);
0081   descriptions.add("materialBudgetVolumeAnalysis", desc);
0082 }
0083 
0084 void MaterialBudgetVolumeAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0085   edm::Handle<edm::MaterialInformationContainer> materialInformationContainer;
0086   iEvent.getByToken(tok_info_, materialInformationContainer);
0087 #ifdef EDM_ML_DEBUG
0088   unsigned int nsize(0), ntot(0), nused(0);
0089 #endif
0090   if (materialInformationContainer.isValid()) {
0091 #ifdef EDM_ML_DEBUG
0092     nsize = materialInformationContainer->size();
0093 #endif
0094     for (const auto& it : *(materialInformationContainer.product())) {
0095 #ifdef EDM_ML_DEBUG
0096       ntot++;
0097 #endif
0098       if (std::find(names_.begin(), names_.end(), it.vname()) != names_.end()) {
0099 #ifdef EDM_ML_DEBUG
0100         nused++;
0101 #endif
0102         unsigned int k =
0103             static_cast<unsigned int>(std::find(names_.begin(), names_.end(), it.vname()) - names_.begin());
0104         meStepEta_[k]->Fill(it.trackEta(), it.stepLength());
0105         meRadLEta_[k]->Fill(it.trackEta(), it.radiationLength());
0106         meIntLEta_[k]->Fill(it.trackEta(), it.interactionLength());
0107         meStepPhi_[k]->Fill(it.trackPhi(), it.stepLength());
0108         meRadLPhi_[k]->Fill(it.trackPhi(), it.radiationLength());
0109         meIntLPhi_[k]->Fill(it.trackPhi(), it.interactionLength());
0110         meStepEtaPhi_[k]->Fill(it.trackEta(), it.trackPhi(), it.stepLength());
0111         meRadLEtaPhi_[k]->Fill(it.trackEta(), it.trackPhi(), it.radiationLength());
0112         meIntLEtaPhi_[k]->Fill(it.trackEta(), it.trackPhi(), it.interactionLength());
0113       }
0114     }
0115   }
0116 #ifdef EDM_ML_DEBUG
0117   edm::LogVerbatim("MaterialBudget") << "MaterialInformation with " << nsize << ":" << ntot << " elements of which "
0118                                      << nused << " are used";
0119 #endif
0120 }
0121 
0122 void MaterialBudgetVolumeAnalysis::bookHisto() {
0123   edm::Service<TFileService> fs;
0124   char name[40], title[100];
0125   for (unsigned int k = 0; k < names_.size(); ++k) {
0126     sprintf(name, "stepEta%s", names_[k].c_str());
0127     sprintf(title, "MB(Step) vs #eta for %s", names_[k].c_str());
0128     meStepEta_.emplace_back(fs->make<TProfile>(name, title, binEta_, etaLow_, etaHigh_));
0129     sprintf(name, "radlEta%s", names_[k].c_str());
0130     sprintf(title, "MB(X0) vs #eta for %s", names_[k].c_str());
0131     meRadLEta_.emplace_back(fs->make<TProfile>(name, title, binEta_, etaLow_, etaHigh_));
0132     sprintf(name, "intlEta%s", names_[k].c_str());
0133     sprintf(title, "MB(L0) vs #eta for %s", names_[k].c_str());
0134     meIntLEta_.emplace_back(fs->make<TProfile>(name, title, binEta_, etaLow_, etaHigh_));
0135     sprintf(name, "stepPhi%s", names_[k].c_str());
0136     sprintf(title, "MB(Step) vs #phi for %s", names_[k].c_str());
0137     meStepPhi_.emplace_back(fs->make<TProfile>(name, title, binPhi_, phiLow_, phiHigh_));
0138     sprintf(name, "radlPhi%s", names_[k].c_str());
0139     sprintf(title, "MB(X0) vs #phi for %s", names_[k].c_str());
0140     meRadLPhi_.emplace_back(fs->make<TProfile>(name, title, binPhi_, phiLow_, phiHigh_));
0141     sprintf(name, "intlPhi%s", names_[k].c_str());
0142     sprintf(title, "MB(L0) vs #phi for %s", names_[k].c_str());
0143     meIntLPhi_.emplace_back(fs->make<TProfile>(name, title, binPhi_, phiLow_, phiHigh_));
0144     sprintf(name, "stepEtaPhi%s", names_[k].c_str());
0145     sprintf(title, "MB(Step) vs #eta and #phi for %s", names_[k].c_str());
0146     meStepEtaPhi_.emplace_back(
0147         fs->make<TProfile2D>(name, title, binEta_ / 2, etaLow_, etaHigh_, binPhi_ / 2, phiLow_, phiHigh_));
0148     sprintf(name, "radlEtaPhi%s", names_[k].c_str());
0149     sprintf(title, "MB(X0) vs #eta and #phi for %s", names_[k].c_str());
0150     meRadLEtaPhi_.emplace_back(
0151         fs->make<TProfile2D>(name, title, binEta_ / 2, etaLow_, etaHigh_, binPhi_ / 2, phiLow_, phiHigh_));
0152     sprintf(name, "intlEtaPhi%s", names_[k].c_str());
0153     sprintf(title, "MB(L0) vs #eta and #phi for %s", names_[k].c_str());
0154     meIntLEtaPhi_.emplace_back(
0155         fs->make<TProfile2D>(name, title, binEta_ / 2, etaLow_, etaHigh_, binPhi_ / 2, phiLow_, phiHigh_));
0156   }
0157 }
0158 
0159 // define this as a plug-in
0160 #include "FWCore/Framework/interface/MakerMacros.h"
0161 DEFINE_FWK_MODULE(MaterialBudgetVolumeAnalysis);