Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:06:33 

0001 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0002 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0003 
0004 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0012 #include "FWCore/Utilities/interface/Exception.h"
0013 #include "FWCore/Utilities/interface/InputTag.h"
0014 
0015 #include "FWCore/ServiceRegistry/interface/Service.h"
0016 
0017 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0018 #include "Geometry/HcalCommonData/interface/HcalHitRelabeller.h"
0019 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0020 
0021 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0022 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0023 
0024 #include "TH1D.h"
0025 #include "TH2D.h"
0026 #include "TProfile.h"
0027 
0028 #include <fstream>
0029 #include <iostream>
0030 #include <map>
0031 #include <string>
0032 #include <vector>
0033 
0034 class HcalSimHitCheck : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0035 public:
0036   HcalSimHitCheck(const edm::ParameterSet &ps);
0037 
0038   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0039 
0040 protected:
0041   void beginRun(edm::Run const &, edm::EventSetup const &) override;
0042   void analyze(edm::Event const &, edm::EventSetup const &) override;
0043   void endRun(edm::Run const &, edm::EventSetup const &) override {}
0044 
0045   void analyzeHits(std::vector<PCaloHit> &);
0046 
0047 private:
0048   const HcalDDDRecConstants *hcons_;
0049   int maxDepthHB_, maxDepthHE_;
0050   int maxDepthHO_, maxDepthHF_;
0051   int maxDepth_;
0052 
0053   int iphi_bins;
0054   float iphi_min, iphi_max;
0055   int ieta_bins_HB;
0056   float ieta_min_HB, ieta_max_HB;
0057   int ieta_bins_HE;
0058   float ieta_min_HE, ieta_max_HE;
0059   int ieta_bins_HO;
0060   float ieta_min_HO, ieta_max_HO;
0061   int ieta_bins_HF;
0062   float ieta_min_HF, ieta_max_HF;
0063 
0064   const std::string g4Label, hcalHits, outFile_;
0065   const int verbose_;
0066   const bool checkHit_, testNumber_, hep17_;
0067 
0068   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hits_;
0069   const edm::ESGetToken<HcalDDDRecConstants, HcalRecNumberingRecord> tok_HRNDC_;
0070 
0071   TH1D *meAllNHit_, *meBadDetHit_, *meBadSubHit_, *meBadIdHit_;
0072   TH1D *meHBNHit_, *meHENHit_, *meHONHit_, *meHFNHit_;
0073   TH1D *meDetectHit_, *meSubdetHit_, *meDepthHit_, *meEtaHit_;
0074   TH2D *meEtaPhiHit_;
0075   std::vector<TH2D *> meEtaPhiHitDepth_;
0076   TH1D *mePhiHit_, *mePhiHitb_, *meEnergyHit_, *meTimeHit_, *meTimeWHit_;
0077   TH1D *meHBDepHit_, *meHEDepHit_, *meHODepHit_, *meHFDepHit_, *meHFDepHitw_;
0078   TH1D *meHBEtaHit_, *meHEEtaHit_, *meHOEtaHit_, *meHFEtaHit_;
0079   TH1D *meHBPhiHit_, *meHEPhiHit_, *meHOPhiHit_, *meHFPhiHit_;
0080   TH1D *meHBEneHit_, *meHEEneHit_, *meHOEneHit_, *meHFEneHit_;
0081   TH2D *meHBEneMap_, *meHEEneMap_, *meHOEneMap_, *meHFEneMap_;
0082   TH1D *meHBEneSum_, *meHEEneSum_, *meHOEneSum_, *meHFEneSum_;
0083   TProfile *meHBEneSum_vs_ieta_, *meHEEneSum_vs_ieta_, *meHOEneSum_vs_ieta_, *meHFEneSum_vs_ieta_;
0084   TH1D *meHBTimHit_, *meHETimHit_, *meHOTimHit_, *meHFTimHit_;
0085   TH1D *meHBEneHit2_, *meHEEneHit2_, *meHOEneHit2_, *meHFEneHit2_;
0086   TH1D *meHBL10Ene_, *meHEL10Ene_, *meHOL10Ene_, *meHFL10Ene_;
0087   TProfile *meHBL10EneP_, *meHEL10EneP_, *meHOL10EneP_, *meHFL10EneP_;
0088   TH1D *meHEP17EneHit_, *meHEP17EneHit2_;
0089 };
0090 
0091 HcalSimHitCheck::HcalSimHitCheck(const edm::ParameterSet &ps)
0092     : g4Label(ps.getParameter<std::string>("moduleLabel")),
0093       hcalHits(ps.getParameter<std::string>("HitCollection")),
0094       outFile_(ps.getParameter<std::string>("outputFile")),
0095       verbose_(ps.getParameter<int>("Verbose")),
0096       checkHit_(true),
0097       testNumber_(ps.getParameter<bool>("TestNumber")),
0098       hep17_(ps.getParameter<bool>("hep17")),
0099       tok_hits_(consumes<edm::PCaloHitContainer>(edm::InputTag(g4Label, hcalHits))),
0100       tok_HRNDC_(esConsumes<HcalDDDRecConstants, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {
0101   edm::LogVerbatim("HcalSim") << "Module Label: " << g4Label << "   Hits: " << hcalHits << " / " << checkHit_
0102                               << "   Output: " << outFile_;
0103 }
0104 
0105 void HcalSimHitCheck::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0106   edm::ParameterSetDescription desc;
0107   desc.add<std::string>("moduleLabel", "g4SimHits");
0108   desc.add<std::string>("HitCollection", "HcalHits");
0109   desc.add<std::string>("outputFile", "hcHit.root");
0110   desc.add<int>("Verbose", 0);
0111   desc.add<bool>("TestNumber", true);
0112   desc.add<bool>("hep17", false);
0113   descriptions.add("hcalSimHitCheck", desc);
0114 }
0115 
0116 void HcalSimHitCheck::beginRun(const edm::Run &, const edm::EventSetup &es) {
0117   hcons_ = &es.getData(tok_HRNDC_);
0118   maxDepthHB_ = hcons_->getMaxDepth(0);
0119   maxDepthHE_ = hcons_->getMaxDepth(1);
0120   maxDepthHF_ = hcons_->getMaxDepth(2);
0121   maxDepthHO_ = hcons_->getMaxDepth(3);
0122   maxDepth_ = (maxDepthHB_ > maxDepthHE_ ? maxDepthHB_ : maxDepthHE_);
0123   maxDepth_ = (maxDepth_ > maxDepthHF_ ? maxDepth_ : maxDepthHF_);
0124   maxDepth_ = (maxDepth_ > maxDepthHO_ ? maxDepth_ : maxDepthHO_);
0125 
0126   // Get Phi segmentation from geometry, use the max phi number so that all iphi
0127   // values are included.
0128 
0129   int NphiMax = hcons_->getNPhi(0);
0130 
0131   NphiMax = (hcons_->getNPhi(1) > NphiMax ? hcons_->getNPhi(1) : NphiMax);
0132   NphiMax = (hcons_->getNPhi(2) > NphiMax ? hcons_->getNPhi(2) : NphiMax);
0133   NphiMax = (hcons_->getNPhi(3) > NphiMax ? hcons_->getNPhi(3) : NphiMax);
0134 
0135   // Center the iphi bins on the integers
0136   iphi_min = 0.5;
0137   iphi_max = NphiMax + 0.5;
0138   iphi_bins = (int)(iphi_max - iphi_min);
0139 
0140   int iEtaHBMax = hcons_->getEtaRange(0).second;
0141   int iEtaHEMax = std::max(hcons_->getEtaRange(1).second, 1);
0142   int iEtaHFMax = hcons_->getEtaRange(2).second;
0143   int iEtaHOMax = hcons_->getEtaRange(3).second;
0144 
0145   // Retain classic behavior, all plots have same ieta range.
0146   // Comment out code to allow each subdetector to have its on range
0147 
0148   int iEtaMax = (iEtaHBMax > iEtaHEMax ? iEtaHBMax : iEtaHEMax);
0149   iEtaMax = (iEtaMax > iEtaHFMax ? iEtaMax : iEtaHFMax);
0150   iEtaMax = (iEtaMax > iEtaHOMax ? iEtaMax : iEtaHOMax);
0151 
0152   iEtaHBMax = iEtaMax;
0153   iEtaHEMax = iEtaMax;
0154   iEtaHFMax = iEtaMax;
0155   iEtaHOMax = iEtaMax;
0156 
0157   // Give an empty bin around the subdet ieta range to make it clear that all
0158   // ieta rings have been included
0159   ieta_min_HB = -iEtaHBMax - 1.5;
0160   ieta_max_HB = iEtaHBMax + 1.5;
0161   ieta_bins_HB = (int)(ieta_max_HB - ieta_min_HB);
0162 
0163   ieta_min_HE = -iEtaHEMax - 1.5;
0164   ieta_max_HE = iEtaHEMax + 1.5;
0165   ieta_bins_HE = (int)(ieta_max_HE - ieta_min_HE);
0166 
0167   ieta_min_HF = -iEtaHFMax - 1.5;
0168   ieta_max_HF = iEtaHFMax + 1.5;
0169   ieta_bins_HF = (int)(ieta_max_HF - ieta_min_HF);
0170 
0171   ieta_min_HO = -iEtaHOMax - 1.5;
0172   ieta_max_HO = iEtaHOMax + 1.5;
0173   ieta_bins_HO = (int)(ieta_max_HO - ieta_min_HO);
0174 
0175   Char_t hname[100];
0176   Char_t htitle[100];
0177 
0178   edm::Service<TFileService> fs;
0179 
0180   // Histograms for Hits
0181   if (checkHit_) {
0182     meAllNHit_ = fs->make<TH1D>("Hit01", "Number of Hits in HCal", 20000, 0., 20000.);
0183     meBadDetHit_ = fs->make<TH1D>("Hit02", "Hits with wrong Det", 100, 0., 100.);
0184     meBadSubHit_ = fs->make<TH1D>("Hit03", "Hits with wrong Subdet", 100, 0., 100.);
0185     meBadIdHit_ = fs->make<TH1D>("Hit04", "Hits with wrong ID", 100, 0., 100.);
0186     meHBNHit_ = fs->make<TH1D>("Hit05", "Number of Hits in HB", 20000, 0., 20000.);
0187     meHENHit_ = fs->make<TH1D>("Hit06", "Number of Hits in HE", 10000, 0., 10000.);
0188     meHONHit_ = fs->make<TH1D>("Hit07", "Number of Hits in HO", 10000, 0., 10000.);
0189     meHFNHit_ = fs->make<TH1D>("Hit08", "Number of Hits in HF", 10000, 0., 10000.);
0190     meHBNHit_->Sumw2();
0191     meHENHit_->Sumw2();
0192     meHONHit_->Sumw2();
0193     meHFNHit_->Sumw2();
0194     meDetectHit_ = fs->make<TH1D>("Hit09", "Detector ID", 50, 0., 50.);
0195     meSubdetHit_ = fs->make<TH1D>("Hit10", "Subdetectors in HCal", 50, 0., 50.);
0196     meDepthHit_ = fs->make<TH1D>("Hit11", "Depths in HCal", 20, 0., 20.);
0197     meEtaHit_ = fs->make<TH1D>("Hit12", "Eta in HCal", ieta_bins_HF, ieta_min_HF, ieta_max_HF);
0198     meEtaPhiHit_ = fs->make<TH2D>(
0199         "Hit12b", "Eta-phi in HCal", ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max);
0200     for (int depth = 1; depth <= maxDepth_; depth++) {
0201       sprintf(hname, "Hit12bd%d", depth);
0202       sprintf(htitle, "Eta-phi in HCal d%d", depth);
0203       meEtaPhiHitDepth_.emplace_back(
0204           fs->make<TH2D>(hname, htitle, ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max));
0205     }
0206     // KC: There are different phi segmentation schemes, this plot uses wider
0207     // bins to represent the most sparse segmentation
0208     mePhiHit_ = fs->make<TH1D>("Hit13", "Phi in HCal (HB,HO)", iphi_bins, iphi_min, iphi_max);
0209     mePhiHitb_ = fs->make<TH1D>("Hit13b", "Phi in HCal (HE,HF)", iphi_bins, iphi_min, iphi_max);
0210     meEnergyHit_ = fs->make<TH1D>("Hit14", "Energy in HCal", 2000, 0., 20.);
0211     meTimeHit_ = fs->make<TH1D>("Hit15", "Time in HCal", 528, 0., 528.);
0212     meTimeWHit_ = fs->make<TH1D>("Hit16", "Time in HCal (E wtd)", 528, 0., 528.);
0213     meHBDepHit_ = fs->make<TH1D>("Hit17", "Depths in HB", 20, 0., 20.);
0214     meHEDepHit_ = fs->make<TH1D>("Hit18", "Depths in HE", 20, 0., 20.);
0215     meHODepHit_ = fs->make<TH1D>("Hit19", "Depths in HO", 20, 0., 20.);
0216     meHFDepHit_ = fs->make<TH1D>("Hit20", "Depths in HF", 20, 0., 20.);
0217     meHBDepHit_->Sumw2();
0218     meHEDepHit_->Sumw2();
0219     meHODepHit_->Sumw2();
0220     meHFDepHit_->Sumw2();
0221     meHFDepHitw_ = fs->make<TH1D>("Hit20b", "Depths in HF (p.e. weighted)", 20, 0., 20.);
0222     meHBEtaHit_ = fs->make<TH1D>("Hit21", "Eta in HB", ieta_bins_HB, ieta_min_HB, ieta_max_HB);
0223     meHEEtaHit_ = fs->make<TH1D>("Hit22", "Eta in HE", ieta_bins_HE, ieta_min_HE, ieta_max_HE);
0224     meHOEtaHit_ = fs->make<TH1D>("Hit23", "Eta in HO", ieta_bins_HO, ieta_min_HO, ieta_max_HO);
0225     meHFEtaHit_ = fs->make<TH1D>("Hit24", "Eta in HF", ieta_bins_HF, ieta_min_HF, ieta_max_HF);
0226     meHBEtaHit_->Sumw2();
0227     meHEEtaHit_->Sumw2();
0228     meHOEtaHit_->Sumw2();
0229     meHFEtaHit_->Sumw2();
0230     meHBPhiHit_ = fs->make<TH1D>("Hit25", "Phi in HB", iphi_bins, iphi_min, iphi_max);
0231     meHEPhiHit_ = fs->make<TH1D>("Hit26", "Phi in HE", iphi_bins, iphi_min, iphi_max);
0232     meHOPhiHit_ = fs->make<TH1D>("Hit27", "Phi in HO", iphi_bins, iphi_min, iphi_max);
0233     meHFPhiHit_ = fs->make<TH1D>("Hit28", "Phi in HF", iphi_bins, iphi_min, iphi_max);
0234     meHBPhiHit_->Sumw2();
0235     meHEPhiHit_->Sumw2();
0236     meHOPhiHit_->Sumw2();
0237     meHFPhiHit_->Sumw2();
0238     meHBEneHit_ = fs->make<TH1D>("Hit29", "Energy in HB", 2000, 0., 20.);
0239     meHEEneHit_ = fs->make<TH1D>("Hit30", "Energy in HE", 500, 0., 5.);
0240     meHEP17EneHit_ = fs->make<TH1D>("Hit30b", "Energy in HEP17", 500, 0., 5.);
0241     meHOEneHit_ = fs->make<TH1D>("Hit31", "Energy in HO", 500, 0., 5.);
0242     meHFEneHit_ = fs->make<TH1D>("Hit32", "Energy in HF", 1001, -0.5, 1000.5);
0243     meHBEneHit_->Sumw2();
0244     meHEEneHit_->Sumw2();
0245     meHOEneHit_->Sumw2();
0246     meHFEneHit_->Sumw2();
0247 
0248     // HxEneMap, HxEneSum, HxEneSum_vs_ieta plot the sum of the simhits energy
0249     // within a single ieta-iphi tower.
0250 
0251     meHBEneMap_ =
0252         fs->make<TH2D>("HBEneMap", "HBEneMap", ieta_bins_HB, ieta_min_HB, ieta_max_HB, iphi_bins, iphi_min, iphi_max);
0253     meHEEneMap_ =
0254         fs->make<TH2D>("HEEneMap", "HEEneMap", ieta_bins_HE, ieta_min_HE, ieta_max_HE, iphi_bins, iphi_min, iphi_max);
0255     meHOEneMap_ =
0256         fs->make<TH2D>("HOEneMap", "HOEneMap", ieta_bins_HO, ieta_min_HO, ieta_max_HO, iphi_bins, iphi_min, iphi_max);
0257     meHFEneMap_ =
0258         fs->make<TH2D>("HFEneMap", "HFEneMap", ieta_bins_HF, ieta_min_HF, ieta_max_HF, iphi_bins, iphi_min, iphi_max);
0259 
0260     meHBEneSum_ = fs->make<TH1D>("HBEneSum", "HBEneSum", 2000, 0., 20.);
0261     meHEEneSum_ = fs->make<TH1D>("HEEneSum", "HEEneSum", 500, 0., 5.);
0262     meHOEneSum_ = fs->make<TH1D>("HOEneSum", "HOEneSum", 500, 0., 5.);
0263     meHFEneSum_ = fs->make<TH1D>("HFEneSum", "HFEneSum", 1001, -0.5, 1000.5);
0264     meHBEneSum_->Sumw2();
0265     meHEEneSum_->Sumw2();
0266     meHOEneSum_->Sumw2();
0267     meHFEneSum_->Sumw2();
0268 
0269     meHBEneSum_vs_ieta_ = fs->make<TProfile>(
0270         "HBEneSum_vs_ieta", "HBEneSum_vs_ieta", ieta_bins_HB, ieta_min_HB, ieta_max_HB, -10.5, 2000.5);
0271     meHEEneSum_vs_ieta_ = fs->make<TProfile>(
0272         "HEEneSum_vs_ieta", "HEEneSum_vs_ieta", ieta_bins_HE, ieta_min_HE, ieta_max_HE, -10.5, 2000.5);
0273     meHOEneSum_vs_ieta_ = fs->make<TProfile>(
0274         "HOEneSum_vs_ieta", "HOEneSum_vs_ieta", ieta_bins_HO, ieta_min_HO, ieta_max_HO, -10.5, 2000.5);
0275     meHFEneSum_vs_ieta_ = fs->make<TProfile>(
0276         "HFEneSum_vs_ieta", "HFEneSum_vs_ieta", ieta_bins_HF, ieta_min_HF, ieta_max_HF, -10.5, 2000.5);
0277 
0278     meHBTimHit_ = fs->make<TH1D>("Hit33", "Time in HB", 528, 0., 528.);
0279     meHETimHit_ = fs->make<TH1D>("Hit34", "Time in HE", 528, 0., 528.);
0280     meHOTimHit_ = fs->make<TH1D>("Hit35", "Time in HO", 528, 0., 528.);
0281     meHFTimHit_ = fs->make<TH1D>("Hit36", "Time in HF", 528, 0., 528.);
0282     meHBTimHit_->Sumw2();
0283     meHETimHit_->Sumw2();
0284     meHOTimHit_->Sumw2();
0285     meHFTimHit_->Sumw2();
0286     // These are the zoomed in energy ranges
0287     meHBEneHit2_ = fs->make<TH1D>("Hit37", "Energy in HB 2", 100, 0., 0.0001);
0288     meHEEneHit2_ = fs->make<TH1D>("Hit38", "Energy in HE 2", 100, 0., 0.0001);
0289     meHEP17EneHit2_ = fs->make<TH1D>("Hit38b", "Energy in HEP17 2", 100, 0., 0.0001);
0290     meHOEneHit2_ = fs->make<TH1D>("Hit39", "Energy in HO 2", 100, 0., 0.0001);
0291     meHFEneHit2_ = fs->make<TH1D>("Hit40", "Energy in HF 2", 100, 0.5, 100.5);
0292     meHBL10Ene_ = fs->make<TH1D>("Hit41", "Log10Energy in HB", 140, -10., 4.);
0293     meHEL10Ene_ = fs->make<TH1D>("Hit42", "Log10Energy in HE", 140, -10., 4.);
0294     meHFL10Ene_ = fs->make<TH1D>("Hit43", "Log10Energy in HF", 50, -1., 4.);
0295     meHOL10Ene_ = fs->make<TH1D>("Hit44", "Log10Energy in HO", 140, -10., 4.);
0296     meHBL10EneP_ = fs->make<TProfile>("Hit45", "Log10Energy in HB vs Hit contribution", 140, -10., 4., 0., 1.);
0297     meHEL10EneP_ = fs->make<TProfile>("Hit46", "Log10Energy in HE vs Hit contribution", 140, -10., 4., 0., 1.);
0298     meHFL10EneP_ = fs->make<TProfile>("Hit47", "Log10Energy in HF vs Hit contribution", 140, -10., 4., 0., 1.);
0299     meHOL10EneP_ = fs->make<TProfile>("Hit48", "Log10Energy in HO vs Hit contribution", 140, -10., 4., 0., 1.);
0300   }
0301 }
0302 
0303 void HcalSimHitCheck::analyze(const edm::Event &e, const edm::EventSetup &) {
0304   if (verbose_ > 0)
0305     edm::LogVerbatim("HcalSim") << "Run = " << e.id().run() << " Event = " << e.id().event();
0306 
0307   bool getHits = false;
0308   if (checkHit_) {
0309     const edm::Handle<edm::PCaloHitContainer> &hitsHcal = e.getHandle(tok_hits_);
0310     if (hitsHcal.isValid()) {
0311       getHits = true;
0312       if (verbose_ > 0)
0313         edm::LogVerbatim("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
0314       std::vector<PCaloHit> caloHits;
0315       caloHits.insert(caloHits.end(), hitsHcal->begin(), hitsHcal->end());
0316       if (verbose_ > 0)
0317         edm::LogVerbatim("HcalSim") << "HcalValidation: Hit buffer " << caloHits.size();
0318       analyzeHits(caloHits);
0319     } else {
0320       if (verbose_ > 0)
0321         edm::LogVerbatim("HcalSim") << "HcalValidation: Input flags Hits " << getHits;
0322     }
0323   }
0324 }
0325 
0326 void HcalSimHitCheck::analyzeHits(std::vector<PCaloHit> &hits) {
0327   int nHit = hits.size();
0328   int nHB = 0, nHE = 0, nHO = 0, nHF = 0, nBad1 = 0, nBad2 = 0, nBad = 0;
0329   std::vector<double> encontHB(140, 0.);
0330   std::vector<double> encontHE(140, 0.);
0331   std::vector<double> encontHF(140, 0.);
0332   std::vector<double> encontHO(140, 0.);
0333   double entotHB = 0, entotHE = 0, entotHF = 0, entotHO = 0;
0334 
0335   double HBEneMap[ieta_bins_HB][iphi_bins];
0336   double HEEneMap[ieta_bins_HE][iphi_bins];
0337   double HOEneMap[ieta_bins_HO][iphi_bins];
0338   double HFEneMap[ieta_bins_HF][iphi_bins];
0339 
0340   // Works in ieta_min_Hx is < 0
0341   int eta_offset_HB = -(int)ieta_min_HB;
0342   int eta_offset_HE = -(int)ieta_min_HE;
0343   int eta_offset_HO = -(int)ieta_min_HO;
0344   int eta_offset_HF = -(int)ieta_min_HF;
0345 
0346   for (int i = 0; i < ieta_bins_HB; i++) {
0347     for (int j = 0; j < iphi_bins; j++) {
0348       HBEneMap[i][j] = 0.;
0349     }
0350   }
0351 
0352   for (int i = 0; i < ieta_bins_HE; i++) {
0353     for (int j = 0; j < iphi_bins; j++) {
0354       HEEneMap[i][j] = 0.;
0355     }
0356   }
0357 
0358   for (int i = 0; i < ieta_bins_HO; i++) {
0359     for (int j = 0; j < iphi_bins; j++) {
0360       HOEneMap[i][j] = 0.;
0361     }
0362   }
0363 
0364   for (int i = 0; i < ieta_bins_HF; i++) {
0365     for (int j = 0; j < iphi_bins; j++) {
0366       HFEneMap[i][j] = 0.;
0367     }
0368   }
0369 
0370   for (int i = 0; i < nHit; i++) {
0371     double energy = hits[i].energy();
0372     double log10en = log10(energy);
0373     int log10i = int((log10en + 10.) * 10.);
0374     double time = hits[i].time();
0375     unsigned int id = hits[i].id();
0376     int det, subdet, depth, eta, phi;
0377     HcalDetId hid;
0378     if (testNumber_)
0379       hid = HcalHitRelabeller::relabel(id, hcons_);
0380     else
0381       hid = HcalDetId(id);
0382     det = hid.det();
0383     subdet = hid.subdet();
0384     depth = hid.depth();
0385     eta = hid.ieta();
0386     phi = hid.iphi();
0387 
0388     if (verbose_ > 1)
0389       edm::LogVerbatim("HcalSim") << "Hit[" << i << "] ID " << std::hex << id << std::dec << " Det " << det << " Sub "
0390                                   << subdet << " depth " << depth << " Eta " << eta << " Phi " << phi << " E " << energy
0391                                   << " time " << time;
0392     if (det == 4) {  // Check DetId.h
0393       if (subdet == static_cast<int>(HcalBarrel))
0394         nHB++;
0395       else if (subdet == static_cast<int>(HcalEndcap))
0396         nHE++;
0397       else if (subdet == static_cast<int>(HcalOuter))
0398         nHO++;
0399       else if (subdet == static_cast<int>(HcalForward))
0400         nHF++;
0401       else {
0402         nBad++;
0403         nBad2++;
0404       }
0405     } else {
0406       nBad++;
0407       nBad1++;
0408     }
0409 
0410     meDetectHit_->Fill(double(det));
0411     if (det == 4) {
0412       meSubdetHit_->Fill(double(subdet));
0413       meDepthHit_->Fill(double(depth));
0414       meEtaHit_->Fill(double(eta));
0415       meEtaPhiHit_->Fill(double(eta), double(phi));
0416       meEtaPhiHitDepth_[depth - 1]->Fill(double(eta), double(phi));
0417 
0418       // We will group the phi plots by HB,HO and HE,HF since these groups share
0419       // similar segmentation schemes
0420       if (subdet == static_cast<int>(HcalBarrel))
0421         mePhiHit_->Fill(double(phi));
0422       else if (subdet == static_cast<int>(HcalEndcap))
0423         mePhiHitb_->Fill(double(phi));
0424       else if (subdet == static_cast<int>(HcalOuter))
0425         mePhiHit_->Fill(double(phi));
0426       else if (subdet == static_cast<int>(HcalForward))
0427         mePhiHitb_->Fill(double(phi));
0428 
0429       // KC: HF energy is in photoelectrons rather than eV, so it will not be
0430       // included in total HCal energy
0431       if (subdet != static_cast<int>(HcalForward)) {
0432         meEnergyHit_->Fill(energy);
0433 
0434         // Since the HF energy is a different scale it does not make sense to
0435         // include it in the Energy Weighted Plot
0436         meTimeWHit_->Fill(double(time), energy);
0437       }
0438       meTimeHit_->Fill(time);
0439 
0440       if (subdet == static_cast<int>(HcalBarrel)) {
0441         meHBDepHit_->Fill(double(depth));
0442         meHBEtaHit_->Fill(double(eta));
0443         meHBPhiHit_->Fill(double(phi));
0444         meHBEneHit_->Fill(energy);
0445         meHBEneHit2_->Fill(energy);
0446         meHBTimHit_->Fill(time);
0447         meHBL10Ene_->Fill(log10en);
0448         if (log10i >= 0 && log10i < 140)
0449           encontHB[log10i] += energy;
0450         entotHB += energy;
0451 
0452         HBEneMap[eta + eta_offset_HB][phi - 1] += energy;
0453 
0454       } else if (subdet == static_cast<int>(HcalEndcap)) {
0455         meHEDepHit_->Fill(double(depth));
0456         meHEEtaHit_->Fill(double(eta));
0457         meHEPhiHit_->Fill(double(phi));
0458 
0459         bool isHEP17 = (phi >= 63) && (phi <= 66) && (eta > 0);
0460         if (hep17_) {
0461           if (!isHEP17) {
0462             meHEEneHit_->Fill(energy);
0463             meHEEneHit2_->Fill(energy);
0464           } else {
0465             meHEP17EneHit_->Fill(energy);
0466             meHEP17EneHit2_->Fill(energy);
0467           }
0468         } else {
0469           meHEEneHit_->Fill(energy);
0470           meHEEneHit2_->Fill(energy);
0471         }
0472 
0473         meHETimHit_->Fill(time);
0474         meHEL10Ene_->Fill(log10en);
0475         if (log10i >= 0 && log10i < 140)
0476           encontHE[log10i] += energy;
0477         entotHE += energy;
0478 
0479         HEEneMap[eta + eta_offset_HE][phi - 1] += energy;
0480 
0481       } else if (subdet == static_cast<int>(HcalOuter)) {
0482         meHODepHit_->Fill(double(depth));
0483         meHOEtaHit_->Fill(double(eta));
0484         meHOPhiHit_->Fill(double(phi));
0485         meHOEneHit_->Fill(energy);
0486         meHOEneHit2_->Fill(energy);
0487         meHOTimHit_->Fill(time);
0488         meHOL10Ene_->Fill(log10en);
0489         if (log10i >= 0 && log10i < 140)
0490           encontHO[log10i] += energy;
0491         entotHO += energy;
0492 
0493         HOEneMap[eta + eta_offset_HO][phi - 1] += energy;
0494 
0495       } else if (subdet == static_cast<int>(HcalForward)) {
0496         meHFDepHit_->Fill(double(depth));
0497         meHFDepHitw_->Fill(double(depth), energy);
0498         meHFEtaHit_->Fill(double(eta));
0499         meHFPhiHit_->Fill(double(phi));
0500         meHFEneHit_->Fill(energy);
0501         meHFEneHit2_->Fill(energy);
0502         meHFTimHit_->Fill(time);
0503         meHFL10Ene_->Fill(log10en);
0504         if (log10i >= 0 && log10i < 140)
0505           encontHF[log10i] += energy;
0506         entotHF += energy;
0507 
0508         HFEneMap[eta + eta_offset_HF][phi - 1] += energy;
0509       }
0510     }
0511   }
0512   if (entotHB != 0)
0513     for (int i = 0; i < 140; i++)
0514       meHBL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHB[i] / entotHB);
0515   if (entotHE != 0)
0516     for (int i = 0; i < 140; i++)
0517       meHEL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHE[i] / entotHE);
0518   if (entotHF != 0)
0519     for (int i = 0; i < 140; i++)
0520       meHFL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHF[i] / entotHF);
0521   if (entotHO != 0)
0522     for (int i = 0; i < 140; i++)
0523       meHOL10EneP_->Fill(-10. + (float(i) + 0.5) / 10., encontHO[i] / entotHO);
0524 
0525   meAllNHit_->Fill(double(nHit));
0526   meBadDetHit_->Fill(double(nBad1));
0527   meBadSubHit_->Fill(double(nBad2));
0528   meBadIdHit_->Fill(double(nBad));
0529   meHBNHit_->Fill(double(nHB));
0530   meHENHit_->Fill(double(nHE));
0531   meHONHit_->Fill(double(nHO));
0532   meHFNHit_->Fill(double(nHF));
0533 
0534   for (int i = 0; i < ieta_bins_HB; i++) {
0535     for (int j = 0; j < iphi_bins; j++) {
0536       if (HBEneMap[i][j] != 0) {
0537         meHBEneSum_->Fill(HBEneMap[i][j]);
0538         meHBEneSum_vs_ieta_->Fill((i - eta_offset_HB), HBEneMap[i][j]);
0539         meHBEneMap_->Fill((i - eta_offset_HB), j + 1, HBEneMap[i][j]);
0540       }
0541     }
0542   }
0543 
0544   for (int i = 0; i < ieta_bins_HE; i++) {
0545     for (int j = 0; j < iphi_bins; j++) {
0546       if (HEEneMap[i][j] != 0) {
0547         meHEEneSum_->Fill(HEEneMap[i][j]);
0548         meHEEneSum_vs_ieta_->Fill((i - eta_offset_HE), HEEneMap[i][j]);
0549         meHEEneMap_->Fill((i - eta_offset_HE), j + 1, HEEneMap[i][j]);
0550       }
0551     }
0552   }
0553 
0554   for (int i = 0; i < ieta_bins_HO; i++) {
0555     for (int j = 0; j < iphi_bins; j++) {
0556       if (HOEneMap[i][j] != 0) {
0557         meHOEneSum_->Fill(HOEneMap[i][j]);
0558         meHOEneSum_vs_ieta_->Fill((i - eta_offset_HO), HOEneMap[i][j]);
0559         meHOEneMap_->Fill((i - eta_offset_HO), j + 1, HOEneMap[i][j]);
0560       }
0561     }
0562   }
0563 
0564   for (int i = 0; i < ieta_bins_HF; i++) {
0565     for (int j = 0; j < iphi_bins; j++) {
0566       if (HFEneMap[i][j] != 0) {
0567         meHFEneSum_->Fill(HFEneMap[i][j]);
0568         meHFEneSum_vs_ieta_->Fill((i - eta_offset_HF), HFEneMap[i][j]);
0569         meHFEneMap_->Fill((i - eta_offset_HF), j + 1, HFEneMap[i][j]);
0570       }
0571     }
0572   }
0573 
0574   if (verbose_ > 0)
0575     edm::LogVerbatim("HcalSim") << "HcalSimHitCheck::analyzeHits: HB " << nHB << " HE " << nHE << " HO " << nHO
0576                                 << " HF " << nHF << " Bad " << nBad << " All " << nHit;
0577 }
0578 
0579 #include "FWCore/Framework/interface/MakerMacros.h"
0580 DEFINE_FWK_MODULE(HcalSimHitCheck);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!