Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:28 

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