Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-24 02:22:26

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   // Book histograms
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   // total X0
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   // Fill from the MB collection
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 //define this as a plug-in
0340 DEFINE_FWK_MODULE(MaterialBudgetHcalAnalysis);