File indexing completed on 2023-03-17 11:27:28
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
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
0160 #include "FWCore/Framework/interface/MakerMacros.h"
0161 DEFINE_FWK_MODULE(MaterialBudgetVolumeAnalysis);
