File indexing completed on 2024-04-06 12:32:14
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/ServiceRegistry/interface/Service.h"
0009 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0010
0011 #include "DataFormats/Math/interface/GeantUnits.h"
0012 #include "SimDataFormats/ValidationFormats/interface/MaterialAccountingCalo.h"
0013
0014 #include <TH1F.h>
0015 #include <TH2F.h>
0016 #include <TProfile.h>
0017 #include <TProfile2D.h>
0018
0019 #include <string>
0020 #include <vector>
0021
0022 using namespace geant_units::operators;
0023
0024 class MaterialBudgetHcalAnalysis : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0025 public:
0026 MaterialBudgetHcalAnalysis(const edm::ParameterSet &p);
0027 ~MaterialBudgetHcalAnalysis() override = default;
0028
0029 static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0030
0031 private:
0032 void analyze(edm::Event const &, edm::EventSetup const &) override;
0033 void beginJob() override;
0034
0035 static const uint32_t maxSet_ = 25, maxSet2_ = 9;
0036 const int binEta_, binPhi_;
0037 const double maxEta_, etaLow_, etaHigh_, etaLowMin_, etaLowMax_, etaMidMin_;
0038 const double etaMidMax_, etaHighMin_, etaHighMax_, etaMinP_, etaMaxP_;
0039 const edm::InputTag labelMBCalo_;
0040 const edm::EDGetTokenT<MaterialAccountingCaloCollection> tokMBCalo_;
0041 TH1F *me400_[maxSet_], *me800_[maxSet_], *me1300_[maxSet2_];
0042 TH2F *me1200_[maxSet_], *me1400_[maxSet2_];
0043 TProfile *me100_[maxSet_], *me200_[maxSet_], *me300_[maxSet_];
0044 TProfile *me500_[maxSet_], *me600_[maxSet_], *me700_[maxSet_];
0045 TProfile *me1500_[maxSet2_];
0046 TProfile *me1600_[maxSet_], *me1700_[maxSet_], *me1800_[maxSet_];
0047 TProfile *me1900_[maxSet_], *me2000_[maxSet_], *me2100_[maxSet_];
0048 TProfile *me2200_[maxSet_], *me2300_[maxSet_], *me2400_[maxSet_];
0049 TProfile2D *me900_[maxSet_], *me1000_[maxSet_], *me1100_[maxSet_];
0050 };
0051
0052 MaterialBudgetHcalAnalysis::MaterialBudgetHcalAnalysis(const edm::ParameterSet &p)
0053 : binEta_(p.getParameter<int>("nBinEta")),
0054 binPhi_(p.getParameter<int>("nBinPhi")),
0055 maxEta_(p.getParameter<double>("maxEta")),
0056 etaLow_(p.getParameter<double>("etaLow")),
0057 etaHigh_(p.getParameter<double>("etaHigh")),
0058 etaLowMin_(p.getParameter<double>("etaLowMin")),
0059 etaLowMax_(p.getParameter<double>("etaLowMax")),
0060 etaMidMin_(p.getParameter<double>("etaMidMin")),
0061 etaMidMax_(p.getParameter<double>("etaMidMax")),
0062 etaHighMin_(p.getParameter<double>("etaHighMin")),
0063 etaHighMax_(p.getParameter<double>("etaHighMax")),
0064 etaMinP_(p.getParameter<double>("etaMinP")),
0065 etaMaxP_(p.getParameter<double>("etaMaxP")),
0066 labelMBCalo_(p.getParameter<edm::InputTag>("labelMBCaloLabel")),
0067 tokMBCalo_(consumes<MaterialAccountingCaloCollection>(labelMBCalo_)) {
0068 usesResource(TFileService::kSharedResource);
0069 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalAnalysis: == Eta plot: NX " << binEta_ << " Range "
0070 << -maxEta_ << ":" << maxEta_ << " Phi plot: NX " << binPhi_ << " Range " << -1._pi
0071 << ":" << 1._pi << " (Eta limit " << etaLow_ << ":" << etaHigh_ << ")"
0072 << " Eta range (" << etaLowMin_ << ":" << etaLowMax_ << "), (" << etaMidMin_ << ":"
0073 << etaMidMax_ << "), (" << etaHighMin_ << ":" << etaHighMax_
0074 << ") Debug for eta range " << etaMinP_ << ":" << etaMaxP_;
0075 }
0076
0077 void MaterialBudgetHcalAnalysis::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0078 edm::ParameterSetDescription desc;
0079 desc.add<int>("nBinEta", 260);
0080 desc.add<int>("nBinPhi", 180);
0081 desc.add<double>("maxEta", 5.2);
0082 desc.add<double>("etaLow", -5.2);
0083 desc.add<double>("etaHigh", 5.2);
0084 desc.add<double>("etaMinP", 5.2);
0085 desc.add<double>("etaMaxP", 0.0);
0086 desc.add<double>("etaLowMin", 0.783);
0087 desc.add<double>("etaLowMax", 0.870);
0088 desc.add<double>("etaMidMin", 2.650);
0089 desc.add<double>("etaMidMax", 2.868);
0090 desc.add<double>("etaHighMin", 2.868);
0091 desc.add<double>("etaHighMax", 3.000);
0092 desc.add<edm::InputTag>("labelMBCaloLabel", edm::InputTag("g4SimHits", "HcalMatBCalo"));
0093 descriptions.add("materialBudgetHcalAnalysis", desc);
0094 }
0095
0096 void MaterialBudgetHcalAnalysis::beginJob() {
0097
0098 edm::Service<TFileService> tfile;
0099
0100 if (!tfile.isAvailable())
0101 throw cms::Exception("BadConfig") << "TFileService unavailable: "
0102 << "please add it to config file";
0103
0104 double maxPhi = 1._pi;
0105 edm::LogVerbatim("MaterialBudgetFull") << "MaterialBudgetHcalAnalysis: Booking user histos === with " << binEta_
0106 << " bins in eta from " << -maxEta_ << " to " << maxEta_ << " and " << binPhi_
0107 << " bins in phi from " << -maxPhi << " to " << maxPhi;
0108
0109 std::string iter;
0110 std::string range0 = "(" + std::to_string(etaMidMin_) + ":" + std::to_string(etaMidMax_) + ") ";
0111 std::string range1 = "(" + std::to_string(etaHighMin_) + ":" + std::to_string(etaHighMax_) + ") ";
0112 std::string range2 = "(" + std::to_string(etaLowMin_) + ":" + std::to_string(etaLowMax_) + ") ";
0113
0114 for (uint32_t i = 0; i < maxSet_; i++) {
0115 iter = std::to_string(i);
0116 me100_[i] = tfile->make<TProfile>(
0117 std::to_string(i + 100).c_str(), ("MB(X0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0118 me200_[i] = tfile->make<TProfile>(
0119 std::to_string(i + 200).c_str(), ("MB(L0) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0120 me300_[i] = tfile->make<TProfile>(
0121 std::to_string(i + 300).c_str(), ("MB(Step) prof Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0122 me400_[i] = tfile->make<TH1F>(
0123 std::to_string(i + 400).c_str(), ("Eta in region " + iter).c_str(), binEta_, -maxEta_, maxEta_);
0124 me500_[i] = tfile->make<TProfile>(
0125 std::to_string(i + 500).c_str(), ("MB(X0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0126 me600_[i] = tfile->make<TProfile>(
0127 std::to_string(i + 600).c_str(), ("MB(L0) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0128 me700_[i] = tfile->make<TProfile>(
0129 std::to_string(i + 700).c_str(), ("MB(Step) prof Ph in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0130 me800_[i] =
0131 tfile->make<TH1F>(std::to_string(i + 800).c_str(), ("Phi in region " + iter).c_str(), binPhi_, -maxPhi, maxPhi);
0132 me900_[i] = tfile->make<TProfile2D>(std::to_string(i + 900).c_str(),
0133 ("MB(X0) prof Eta Phi in region " + iter).c_str(),
0134 binEta_ / 2,
0135 -maxEta_,
0136 maxEta_,
0137 binPhi_ / 2,
0138 -maxPhi,
0139 maxPhi);
0140 me1000_[i] = tfile->make<TProfile2D>(std::to_string(i + 1000).c_str(),
0141 ("MB(L0) prof Eta Phi in region " + iter).c_str(),
0142 binEta_ / 2,
0143 -maxEta_,
0144 maxEta_,
0145 binPhi_ / 2,
0146 -maxPhi,
0147 maxPhi);
0148 me1100_[i] = tfile->make<TProfile2D>(std::to_string(i + 1100).c_str(),
0149 ("MB(Step) prof Eta Phi in region " + iter).c_str(),
0150 binEta_ / 2,
0151 -maxEta_,
0152 maxEta_,
0153 binPhi_ / 2,
0154 -maxPhi,
0155 maxPhi);
0156 me1200_[i] = tfile->make<TH2F>(std::to_string(i + 1200).c_str(),
0157 ("Eta vs Phi in region " + iter).c_str(),
0158 binEta_ / 2,
0159 -maxEta_,
0160 maxEta_,
0161 binPhi_ / 2,
0162 -maxPhi,
0163 maxPhi);
0164 me1600_[i] = tfile->make<TProfile>(std::to_string(i + 1600).c_str(),
0165 ("MB(X0) prof Ph in region " + range0 + iter).c_str(),
0166 binPhi_,
0167 -maxPhi,
0168 maxPhi);
0169 me1700_[i] = tfile->make<TProfile>(std::to_string(i + 1700).c_str(),
0170 ("MB(L0) prof Ph in region " + range0 + iter).c_str(),
0171 binPhi_,
0172 -maxPhi,
0173 maxPhi);
0174 me1800_[i] = tfile->make<TProfile>(std::to_string(i + 1800).c_str(),
0175 ("MB(Step) prof Ph in region " + range0 + iter).c_str(),
0176 binPhi_,
0177 -maxPhi,
0178 maxPhi);
0179 me1900_[i] = tfile->make<TProfile>(std::to_string(i + 1900).c_str(),
0180 ("MB(X0) prof Ph in region " + range1 + iter).c_str(),
0181 binPhi_,
0182 -maxPhi,
0183 maxPhi);
0184 me2000_[i] = tfile->make<TProfile>(std::to_string(i + 2000).c_str(),
0185 ("MB(L0) prof Ph in region " + range1 + iter).c_str(),
0186 binPhi_,
0187 -maxPhi,
0188 maxPhi);
0189 me2100_[i] = tfile->make<TProfile>(std::to_string(i + 2100).c_str(),
0190 ("MB(Step) prof Ph in region " + range1 + iter).c_str(),
0191 binPhi_,
0192 -maxPhi,
0193 maxPhi);
0194 me2200_[i] = tfile->make<TProfile>(std::to_string(i + 2200).c_str(),
0195 ("MB(X0) prof Ph in region " + range2 + iter).c_str(),
0196 binPhi_,
0197 -maxPhi,
0198 maxPhi);
0199 me2300_[i] = tfile->make<TProfile>(std::to_string(i + 2300).c_str(),
0200 ("MB(L0) prof Ph in region " + range2 + iter).c_str(),
0201 binPhi_,
0202 -maxPhi,
0203 maxPhi);
0204 me2400_[i] = tfile->make<TProfile>(std::to_string(i + 2400).c_str(),
0205 ("MB(Step) prof Ph in region " + range2 + iter).c_str(),
0206 binPhi_,
0207 -maxPhi,
0208 maxPhi);
0209 }
0210 for (uint32_t i = 0; i < maxSet2_; i++) {
0211 iter = std::to_string(i);
0212 me1300_[i] = tfile->make<TH1F>(std::to_string(i + 1300).c_str(),
0213 ("Events with layers Hit (0 all, 1 HB, ..) for " + iter).c_str(),
0214 binEta_,
0215 -maxEta_,
0216 maxEta_);
0217 me1400_[i] = tfile->make<TH2F>(std::to_string(i + 1400).c_str(),
0218 ("Eta vs Phi for layers hit in " + iter).c_str(),
0219 binEta_ / 2,
0220 -maxEta_,
0221 maxEta_,
0222 binPhi_ / 2,
0223 -maxPhi,
0224 maxPhi);
0225 me1500_[i] = tfile->make<TProfile>(std::to_string(i + 1500).c_str(),
0226 ("Number of layers crossed (0 all, 1 HB, ..) for " + iter).c_str(),
0227 binEta_,
0228 -maxEta_,
0229 maxEta_);
0230 }
0231
0232 edm::LogVerbatim("MaterialBudget") << "MaterialBudgetHcalAnalysis: Booking user histos done ===";
0233 }
0234
0235 void MaterialBudgetHcalAnalysis::analyze(edm::Event const &iEvent, edm::EventSetup const &iSetup) {
0236 #ifdef EDM_ML_DEBUG
0237 edm::LogVerbatim("MaterialBudgetFull") << "Run " << iEvent.id().run() << " Event " << iEvent.id().event()
0238 << " Luminosity " << iEvent.luminosityBlock() << " Bunch "
0239 << iEvent.bunchCrossing();
0240 #endif
0241
0242
0243 auto const &hcalMBColl = iEvent.getHandle(tokMBCalo_);
0244 if (hcalMBColl.isValid()) {
0245 auto hcalMB = hcalMBColl.product();
0246 #ifdef EDM_ML_DEBUG
0247 edm::LogVerbatim("MaterialBudgetFull")
0248 << "Finds HcalMaterialBudgetCollection with " << hcalMB->size() << " entries";
0249 #endif
0250
0251 for (auto itr = hcalMB->begin(); itr != hcalMB->end(); ++itr) {
0252 for (uint32_t ii = 0; ii < itr->m_stepLen.size(); ++ii) {
0253 #ifdef EDM_ML_DEBUG
0254 if ((std::abs(itr->m_eta) >= etaMinP_) && (std::abs(itr->m_eta) <= etaMaxP_))
0255 edm::LogVerbatim("MaterialBudget")
0256 << "MaterialBudgetHcalAnalysis:FillHisto called with index " << ii << " integrated step "
0257 << itr->m_stepLen[ii] << " X0 " << itr->m_radLen[ii] << " Lamda " << itr->m_intLen[ii];
0258 #endif
0259 if (ii < maxSet_) {
0260 me100_[ii]->Fill(itr->m_eta, itr->m_radLen[ii]);
0261 me200_[ii]->Fill(itr->m_eta, itr->m_intLen[ii]);
0262 me300_[ii]->Fill(itr->m_eta, itr->m_stepLen[ii]);
0263 me400_[ii]->Fill(itr->m_eta);
0264
0265 if (itr->m_eta >= etaLow_ && itr->m_eta <= etaHigh_) {
0266 me500_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
0267 me600_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
0268 me700_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
0269 me800_[ii]->Fill(itr->m_phi);
0270 }
0271
0272 me900_[ii]->Fill(itr->m_eta, itr->m_phi, itr->m_radLen[ii]);
0273 me1000_[ii]->Fill(itr->m_eta, itr->m_phi, itr->m_intLen[ii]);
0274 me1100_[ii]->Fill(itr->m_eta, itr->m_phi, itr->m_stepLen[ii]);
0275 me1200_[ii]->Fill(itr->m_eta, itr->m_phi);
0276
0277 if ((std::abs(itr->m_eta) >= etaMidMin_) && (std::abs(itr->m_eta) <= etaMidMax_)) {
0278 me1600_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
0279 me1700_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
0280 me1800_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
0281 }
0282
0283 if ((std::abs(itr->m_eta) >= etaHighMin_) && (std::abs(itr->m_eta) <= etaHighMax_)) {
0284 me1900_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
0285 me2000_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
0286 me2100_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
0287 }
0288
0289 if ((std::abs(itr->m_eta) >= etaLowMin_) && (std::abs(itr->m_eta) <= etaLowMax_)) {
0290 me2200_[ii]->Fill(itr->m_phi, itr->m_radLen[ii]);
0291 me2300_[ii]->Fill(itr->m_phi, itr->m_intLen[ii]);
0292 me2400_[ii]->Fill(itr->m_phi, itr->m_stepLen[ii]);
0293 }
0294 }
0295 }
0296
0297 me1300_[0]->Fill(itr->m_eta);
0298 me1400_[0]->Fill(itr->m_eta, itr->m_phi);
0299 if (itr->m_layers[0] > 0) {
0300 me1300_[1]->Fill(itr->m_eta);
0301 me1400_[1]->Fill(itr->m_eta, itr->m_phi);
0302 }
0303 if (itr->m_layers[0] >= 16) {
0304 me1300_[2]->Fill(itr->m_eta);
0305 me1400_[2]->Fill(itr->m_eta, itr->m_phi);
0306 }
0307 if (itr->m_layers[1] > 0) {
0308 me1300_[3]->Fill(itr->m_eta);
0309 me1400_[3]->Fill(itr->m_eta, itr->m_phi);
0310 }
0311 if (itr->m_layers[1] >= 16) {
0312 me1300_[4]->Fill(itr->m_eta);
0313 me1400_[4]->Fill(itr->m_eta, itr->m_phi);
0314 }
0315 if (itr->m_layers[2] > 0) {
0316 me1300_[5]->Fill(itr->m_eta);
0317 me1400_[5]->Fill(itr->m_eta, itr->m_phi);
0318 }
0319 if (itr->m_layers[2] >= 2) {
0320 me1300_[6]->Fill(itr->m_eta);
0321 me1400_[6]->Fill(itr->m_eta, itr->m_phi);
0322 }
0323 if (itr->m_layers[3] > 0) {
0324 me1300_[7]->Fill(itr->m_eta);
0325 me1400_[7]->Fill(itr->m_eta, itr->m_phi);
0326 }
0327 if (itr->m_layers[0] > 0 || itr->m_layers[1] > 0 || (itr->m_layers[3] > 0 && std::abs(itr->m_eta) > 3.0)) {
0328 me1300_[8]->Fill(itr->m_eta);
0329 me1400_[8]->Fill(itr->m_eta, itr->m_phi);
0330 }
0331 me1500_[0]->Fill(itr->m_eta, (double)(itr->m_layers[0] + itr->m_layers[1] + itr->m_layers[2] + itr->m_layers[3]));
0332 me1500_[1]->Fill(itr->m_eta, (double)(itr->m_layers[0]));
0333 me1500_[2]->Fill(itr->m_eta, (double)(itr->m_layers[1]));
0334 me1500_[4]->Fill(itr->m_eta, (double)(itr->m_layers[3]));
0335 }
0336 }
0337 }
0338
0339
0340 DEFINE_FWK_MODULE(MaterialBudgetHcalAnalysis);
