Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-04-29 23:10:25

0001 #include "DQMOffline/JetMET/interface/HCALRecHitAnalyzer.h"
0002 // author: Bobby Scurlock, University of Florida
0003 // first version 12/7/2006
0004 // modified: Mike Schmitt
0005 // date: 03.05.2007
0006 // note: 1) code rewrite. 2.) changed to loop over all hcal detids;
0007 //       not only those within calotowers.
0008 
0009 #include "FWCore/Framework/interface/MakerMacros.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 
0012 #include "FWCore/PluginManager/interface/ModuleDef.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/ESHandle.h"
0015 #include "DataFormats/Common/interface/Handle.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 
0019 //#include "Geometry/Vector/interface/GlobalPoint.h"
0020 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0021 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0022 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0023 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0024 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0025 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0026 
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0029 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0030 
0031 #include <memory>
0032 #include <vector>
0033 #include <utility>
0034 #include <ostream>
0035 #include <fstream>
0036 #include <string>
0037 #include <algorithm>
0038 #include <cmath>
0039 #include <TLorentzVector.h>
0040 #include "DQMServices/Core/interface/DQMStore.h"
0041 
0042 #define DEBUG(X)                   \
0043   {                                \
0044     if (debug_) {                  \
0045       std::cout << X << std::endl; \
0046     }                              \
0047   }
0048 
0049 HCALRecHitAnalyzer::HCALRecHitAnalyzer(const edm::ParameterSet& iConfig) {
0050   // Retrieve Information from the Configuration File
0051   hBHERecHitsLabel_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("HBHERecHitsLabel"));
0052   hORecHitsLabel_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("HORecHitsLabel"));
0053   hFRecHitsLabel_ = consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("HFRecHitsLabel"));
0054   caloGeomToken_ = esConsumes<edm::Transition::BeginRun>();
0055   debug_ = iConfig.getParameter<bool>("Debug");
0056   finebinning_ = iConfig.getUntrackedParameter<bool>("FineBinning");
0057   FolderName_ = iConfig.getUntrackedParameter<std::string>("FolderName");
0058 }
0059 
0060 void HCALRecHitAnalyzer::dqmBeginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0061   Nevents = 0;
0062   FillGeometry(iSetup);
0063 }
0064 
0065 void HCALRecHitAnalyzer::bookHistograms(DQMStore::IBooker& ibooker, edm::Run const& iRun, edm::EventSetup const&) {
0066   ibooker.setCurrentFolder(FolderName_ + "/geometry");
0067   hHCAL_ieta_iphi_HBMap = ibooker.book2D("METTask_HCAL_ieta_iphi_HBMap", "", 83, -41, 42, 72, 1, 73);
0068   hHCAL_ieta_iphi_HEMap = ibooker.book2D("METTask_HCAL_ieta_iphi_HEMap", "", 83, -41, 42, 72, 1, 73);
0069   hHCAL_ieta_iphi_HFMap = ibooker.book2D("METTask_HCAL_ieta_iphi_HFMap", "", 83, -41, 42, 72, 1, 73);
0070   hHCAL_ieta_iphi_HOMap = ibooker.book2D("METTask_HCAL_ieta_iphi_HOMap", "", 83, -41, 42, 72, 1, 73);
0071   hHCAL_ieta_iphi_etaMap = ibooker.book2D("METTask_HCAL_ieta_iphi_etaMap", "", 83, -41, 42, 72, 1, 73);
0072   hHCAL_ieta_iphi_phiMap = ibooker.book2D("METTask_HCAL_ieta_iphi_phiMap", "", 83, -41, 42, 72, 1, 73);
0073   hHCAL_ieta_detaMap = ibooker.book1D("METTask_HCAL_ieta_detaMap", "", 83, -41, 42);
0074   hHCAL_ieta_dphiMap = ibooker.book1D("METTask_HCAL_ieta_dphiMap", "", 83, -41, 42);
0075 
0076   // Initialize bins for geometry to -999 because z = 0 is a valid entry
0077   for (int i = 1; i <= 83; i++) {
0078     hHCAL_ieta_detaMap->setBinContent(i, -999);
0079     hHCAL_ieta_dphiMap->setBinContent(i, -999);
0080 
0081     for (int j = 1; j <= 72; j++) {
0082       hHCAL_ieta_iphi_HBMap->setBinContent(i, j, 0);
0083       hHCAL_ieta_iphi_HEMap->setBinContent(i, j, 0);
0084       hHCAL_ieta_iphi_HFMap->setBinContent(i, j, 0);
0085       hHCAL_ieta_iphi_HOMap->setBinContent(i, j, 0);
0086       hHCAL_ieta_iphi_etaMap->setBinContent(i, j, -999);
0087       hHCAL_ieta_iphi_phiMap->setBinContent(i, j, -999);
0088     }
0089   }
0090 
0091   //    ibooker.setCurrentFolder("RecoMETV/MET_HCAL/data");
0092   ibooker.setCurrentFolder(FolderName_);
0093   //--Store number of events used
0094   hHCAL_Nevents = ibooker.book1D("METTask_HCAL_Nevents", "", 1, 0, 1);
0095   //--Data integrated over all events and stored by HCAL(ieta,iphi)
0096 
0097   hHCAL_D1_energy_ieta_iphi = ibooker.book2D("METTask_HCAL_D1_energy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0098   hHCAL_D2_energy_ieta_iphi = ibooker.book2D("METTask_HCAL_D2_energy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0099   hHCAL_D3_energy_ieta_iphi = ibooker.book2D("METTask_HCAL_D3_energy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0100   hHCAL_D4_energy_ieta_iphi = ibooker.book2D("METTask_HCAL_D4_energy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0101 
0102   hHCAL_D1_Minenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D1_Minenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0103   hHCAL_D2_Minenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D2_Minenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0104   hHCAL_D3_Minenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D3_Minenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0105   hHCAL_D4_Minenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D4_Minenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0106 
0107   hHCAL_D1_Maxenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D1_Maxenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0108   hHCAL_D2_Maxenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D2_Maxenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0109   hHCAL_D3_Maxenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D3_Maxenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0110   hHCAL_D4_Maxenergy_ieta_iphi = ibooker.book2D("METTask_HCAL_D4_Maxenergy_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0111 
0112   // need to initialize those
0113   for (int i = 1; i <= 83; i++)
0114     for (int j = 1; j <= 73; j++) {
0115       hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(i, j, -999);
0116       hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(i, j, -999);
0117       hHCAL_D3_Maxenergy_ieta_iphi->setBinContent(i, j, -999);
0118       hHCAL_D4_Maxenergy_ieta_iphi->setBinContent(i, j, -999);
0119 
0120       hHCAL_D1_Minenergy_ieta_iphi->setBinContent(i, j, 14000);
0121       hHCAL_D2_Minenergy_ieta_iphi->setBinContent(i, j, 14000);
0122       hHCAL_D3_Minenergy_ieta_iphi->setBinContent(i, j, 14000);
0123       hHCAL_D4_Minenergy_ieta_iphi->setBinContent(i, j, 14000);
0124     }
0125 
0126   hHCAL_D1_Occ_ieta_iphi = ibooker.book2D("METTask_HCAL_D1_Occ_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0127   hHCAL_D2_Occ_ieta_iphi = ibooker.book2D("METTask_HCAL_D2_Occ_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0128   hHCAL_D3_Occ_ieta_iphi = ibooker.book2D("METTask_HCAL_D3_Occ_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0129   hHCAL_D4_Occ_ieta_iphi = ibooker.book2D("METTask_HCAL_D4_Occ_ieta_iphi", "", 83, -41, 42, 72, 1, 73);
0130   //--Data over eta-rings
0131 
0132   // CaloTower values
0133 
0134   if (finebinning_) {
0135     hHCAL_D1_energyvsieta = ibooker.book2D("METTask_HCAL_D1_energyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0136     hHCAL_D2_energyvsieta = ibooker.book2D("METTask_HCAL_D2_energyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0137     hHCAL_D3_energyvsieta = ibooker.book2D("METTask_HCAL_D3_energyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0138     hHCAL_D4_energyvsieta = ibooker.book2D("METTask_HCAL_D4_energyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0139 
0140     hHCAL_D1_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D1_Minenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0141     hHCAL_D2_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D2_Minenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0142     hHCAL_D3_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D3_Minenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0143     hHCAL_D4_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D4_Minenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0144 
0145     hHCAL_D1_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D1_Maxenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0146     hHCAL_D2_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D2_Maxenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0147     hHCAL_D3_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D3_Maxenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0148     hHCAL_D4_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D4_Maxenergyvsieta", "", 83, -41, 42, 20101, -100, 2001);
0149 
0150     // Integrated over phi
0151     hHCAL_D1_Occvsieta = ibooker.book2D("METTask_HCAL_D1_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0152     hHCAL_D2_Occvsieta = ibooker.book2D("METTask_HCAL_D2_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0153     hHCAL_D3_Occvsieta = ibooker.book2D("METTask_HCAL_D3_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0154     hHCAL_D4_Occvsieta = ibooker.book2D("METTask_HCAL_D4_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0155 
0156     hHCAL_D1_SETvsieta = ibooker.book2D("METTask_HCAL_D1_SETvsieta", "", 83, -41, 42, 20001, 0, 2001);
0157     hHCAL_D2_SETvsieta = ibooker.book2D("METTask_HCAL_D2_SETvsieta", "", 83, -41, 42, 20001, 0, 2001);
0158     hHCAL_D3_SETvsieta = ibooker.book2D("METTask_HCAL_D3_SETvsieta", "", 83, -41, 42, 20001, 0, 2001);
0159     hHCAL_D4_SETvsieta = ibooker.book2D("METTask_HCAL_D4_SETvsieta", "", 83, -41, 42, 20001, 0, 2001);
0160 
0161     hHCAL_D1_METvsieta = ibooker.book2D("METTask_HCAL_D1_METvsieta", "", 83, -41, 42, 20001, 0, 2001);
0162     hHCAL_D2_METvsieta = ibooker.book2D("METTask_HCAL_D2_METvsieta", "", 83, -41, 42, 20001, 0, 2001);
0163     hHCAL_D3_METvsieta = ibooker.book2D("METTask_HCAL_D3_METvsieta", "", 83, -41, 42, 20001, 0, 2001);
0164     hHCAL_D4_METvsieta = ibooker.book2D("METTask_HCAL_D4_METvsieta", "", 83, -41, 42, 20001, 0, 2001);
0165 
0166     hHCAL_D1_METPhivsieta = ibooker.book2D("METTask_HCAL_D1_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0167     hHCAL_D2_METPhivsieta = ibooker.book2D("METTask_HCAL_D2_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0168     hHCAL_D3_METPhivsieta = ibooker.book2D("METTask_HCAL_D3_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0169     hHCAL_D4_METPhivsieta = ibooker.book2D("METTask_HCAL_D4_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0170 
0171     hHCAL_D1_MExvsieta = ibooker.book2D("METTask_HCAL_D1_MExvsieta", "", 83, -41, 42, 10001, -500, 501);
0172     hHCAL_D2_MExvsieta = ibooker.book2D("METTask_HCAL_D2_MExvsieta", "", 83, -41, 42, 10001, -500, 501);
0173     hHCAL_D3_MExvsieta = ibooker.book2D("METTask_HCAL_D3_MExvsieta", "", 83, -41, 42, 10001, -500, 501);
0174     hHCAL_D4_MExvsieta = ibooker.book2D("METTask_HCAL_D4_MExvsieta", "", 83, -41, 42, 10001, -500, 501);
0175 
0176     hHCAL_D1_MEyvsieta = ibooker.book2D("METTask_HCAL_D1_MEyvsieta", "", 83, -41, 42, 10001, -500, 501);
0177     hHCAL_D2_MEyvsieta = ibooker.book2D("METTask_HCAL_D2_MEyvsieta", "", 83, -41, 42, 10001, -500, 501);
0178     hHCAL_D3_MEyvsieta = ibooker.book2D("METTask_HCAL_D3_MEyvsieta", "", 83, -41, 42, 10001, -500, 501);
0179     hHCAL_D4_MEyvsieta = ibooker.book2D("METTask_HCAL_D4_MEyvsieta", "", 83, -41, 42, 10001, -500, 501);
0180   } else {
0181     hHCAL_D1_energyvsieta = ibooker.book2D("METTask_HCAL_D1_energyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0182     hHCAL_D2_energyvsieta = ibooker.book2D("METTask_HCAL_D2_energyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0183     hHCAL_D3_energyvsieta = ibooker.book2D("METTask_HCAL_D3_energyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0184     hHCAL_D4_energyvsieta = ibooker.book2D("METTask_HCAL_D4_energyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0185 
0186     hHCAL_D1_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D1_Minenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0187     hHCAL_D2_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D2_Minenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0188     hHCAL_D3_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D3_Minenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0189     hHCAL_D4_Minenergyvsieta = ibooker.book2D("METTask_HCAL_D4_Minenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0190 
0191     hHCAL_D1_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D1_Maxenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0192     hHCAL_D2_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D2_Maxenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0193     hHCAL_D3_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D3_Maxenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0194     hHCAL_D4_Maxenergyvsieta = ibooker.book2D("METTask_HCAL_D4_Maxenergyvsieta", "", 83, -41, 42, 1000, -10, 1990);
0195 
0196     // Integrated over phi
0197     hHCAL_D1_Occvsieta = ibooker.book2D("METTask_HCAL_D1_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0198     hHCAL_D2_Occvsieta = ibooker.book2D("METTask_HCAL_D2_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0199     hHCAL_D3_Occvsieta = ibooker.book2D("METTask_HCAL_D3_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0200     hHCAL_D4_Occvsieta = ibooker.book2D("METTask_HCAL_D4_Occvsieta", "", 83, -41, 42, 73, 0, 73);
0201 
0202     hHCAL_D1_SETvsieta = ibooker.book2D("METTask_HCAL_D1_SETvsieta", "", 83, -41, 42, 1000, 0, 2000);
0203     hHCAL_D2_SETvsieta = ibooker.book2D("METTask_HCAL_D2_SETvsieta", "", 83, -41, 42, 1000, 0, 2000);
0204     hHCAL_D3_SETvsieta = ibooker.book2D("METTask_HCAL_D3_SETvsieta", "", 83, -41, 42, 1000, 0, 2000);
0205     hHCAL_D4_SETvsieta = ibooker.book2D("METTask_HCAL_D4_SETvsieta", "", 83, -41, 42, 1000, 0, 2000);
0206 
0207     hHCAL_D1_METvsieta = ibooker.book2D("METTask_HCAL_D1_METvsieta", "", 83, -41, 42, 1000, 0, 2000);
0208     hHCAL_D2_METvsieta = ibooker.book2D("METTask_HCAL_D2_METvsieta", "", 83, -41, 42, 1000, 0, 2000);
0209     hHCAL_D3_METvsieta = ibooker.book2D("METTask_HCAL_D3_METvsieta", "", 83, -41, 42, 1000, 0, 2000);
0210     hHCAL_D4_METvsieta = ibooker.book2D("METTask_HCAL_D4_METvsieta", "", 83, -41, 42, 1000, 0, 2000);
0211 
0212     hHCAL_D1_METPhivsieta = ibooker.book2D("METTask_HCAL_D1_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0213     hHCAL_D2_METPhivsieta = ibooker.book2D("METTask_HCAL_D2_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0214     hHCAL_D3_METPhivsieta = ibooker.book2D("METTask_HCAL_D3_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0215     hHCAL_D4_METPhivsieta = ibooker.book2D("METTask_HCAL_D4_METPhivsieta", "", 83, -41, 42, 80, -4, 4);
0216 
0217     hHCAL_D1_MExvsieta = ibooker.book2D("METTask_HCAL_D1_MExvsieta", "", 83, -41, 42, 500, -500, 500);
0218     hHCAL_D2_MExvsieta = ibooker.book2D("METTask_HCAL_D2_MExvsieta", "", 83, -41, 42, 500, -500, 500);
0219     hHCAL_D3_MExvsieta = ibooker.book2D("METTask_HCAL_D3_MExvsieta", "", 83, -41, 42, 500, -500, 500);
0220     hHCAL_D4_MExvsieta = ibooker.book2D("METTask_HCAL_D4_MExvsieta", "", 83, -41, 42, 500, -500, 500);
0221 
0222     hHCAL_D1_MEyvsieta = ibooker.book2D("METTask_HCAL_D1_MEyvsieta", "", 83, -41, 42, 500, -500, 500);
0223     hHCAL_D2_MEyvsieta = ibooker.book2D("METTask_HCAL_D2_MEyvsieta", "", 83, -41, 42, 500, -500, 500);
0224     hHCAL_D3_MEyvsieta = ibooker.book2D("METTask_HCAL_D3_MEyvsieta", "", 83, -41, 42, 500, -500, 500);
0225     hHCAL_D4_MEyvsieta = ibooker.book2D("METTask_HCAL_D4_MEyvsieta", "", 83, -41, 42, 500, -500, 500);
0226   }
0227   // Inspect Setup for CaloTower Geometry
0228   //  FillGeometry(iSetup);
0229 }
0230 
0231 void HCALRecHitAnalyzer::FillGeometry(const edm::EventSetup& iSetup) {
0232   // ==========================================================
0233   // Retrieve!
0234   // ==========================================================
0235 
0236   const auto& pG = iSetup.getHandle(caloGeomToken_);
0237 
0238   if (!pG.isValid()) {
0239     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Setup Handle, Aborting Task "
0240                                << "HCALRecHitAnalyzer::FillGeometry!\n";
0241     return;
0242   }
0243 
0244   const CaloGeometry cG = *pG;
0245 
0246   const HcalGeometry* HBgeom = dynamic_cast<const HcalGeometry*>(cG.getSubdetectorGeometry(DetId::Hcal, HcalBarrel));
0247   const HcalGeometry* HEgeom = dynamic_cast<const HcalGeometry*>(cG.getSubdetectorGeometry(DetId::Hcal, HcalEndcap));
0248   const CaloSubdetectorGeometry* HOgeom = cG.getSubdetectorGeometry(DetId::Hcal, HcalOuter);
0249   const CaloSubdetectorGeometry* HFgeom = cG.getSubdetectorGeometry(DetId::Hcal, HcalForward);
0250 
0251   // ==========================================================
0252   // Fill Histograms!
0253   // ==========================================================
0254 
0255   std::vector<DetId>::iterator i;
0256 
0257   int HBmin_ieta = 99, HBmax_ieta = -99;
0258   int HBmin_iphi = 99, HBmax_iphi = -99;
0259 
0260   // Loop Over all Hcal Barrel DetId's
0261   int nHBdetid = 0;
0262   std::vector<DetId> HBids = HBgeom->getValidDetIds(DetId::Hcal, HcalBarrel);
0263 
0264   for (i = HBids.begin(); i != HBids.end(); i++) {
0265     nHBdetid++;
0266 
0267     HcalDetId HcalID(*i);
0268 
0269     int Calo_ieta = 42 + HcalID.ieta();
0270     int Calo_iphi = HcalID.iphi();
0271     double Calo_eta = HBgeom->getPosition(HcalID).eta();
0272     double Calo_phi = HBgeom->getPosition(HcalID).phi();
0273 
0274     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0275       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0276       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0277     }
0278 
0279     if (Calo_ieta > HBmax_ieta)
0280       HBmax_ieta = Calo_ieta;
0281     if (Calo_ieta < HBmin_ieta)
0282       HBmin_ieta = Calo_ieta;
0283     if (Calo_iphi > HBmax_iphi)
0284       HBmax_iphi = Calo_iphi;
0285     if (Calo_iphi > HBmax_iphi)
0286       HBmin_iphi = Calo_iphi;
0287   }
0288 
0289   int HEmin_ieta = 99, HEmax_ieta = -99;
0290   int HEmin_iphi = 99, HEmax_iphi = -99;
0291 
0292   // Loop Over all Hcal Endcap DetId's
0293   int nHEdetid = 0;
0294   std::vector<DetId> HEids = HEgeom->getValidDetIds(DetId::Hcal, HcalEndcap);
0295 
0296   for (i = HEids.begin(); i != HEids.end(); i++) {
0297     nHEdetid++;
0298 
0299     HcalDetId HcalID(*i);
0300 
0301     int Calo_ieta = 42 + HcalID.ieta();
0302     int Calo_iphi = HcalID.iphi();
0303     double Calo_eta = HEgeom->getPosition(HcalID).eta();
0304     double Calo_phi = HEgeom->getPosition(HcalID).phi();
0305 
0306     // HCAL to HE eta, phi map comparison
0307     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0308       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0309       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0310     }
0311 
0312     if (Calo_ieta > HEmax_ieta)
0313       HEmax_ieta = Calo_ieta;
0314     if (Calo_ieta < HEmin_ieta)
0315       HEmin_ieta = Calo_ieta;
0316     if (Calo_iphi > HEmax_iphi)
0317       HEmax_iphi = Calo_iphi;
0318     if (Calo_iphi > HEmax_iphi)
0319       HEmin_iphi = Calo_iphi;
0320   }
0321 
0322   int HFmin_ieta = 99, HFmax_ieta = -99;
0323   int HFmin_iphi = 99, HFmax_iphi = -99;
0324 
0325   // Loop Over all Hcal Forward DetId's
0326   int nHFdetid = 0;
0327   std::vector<DetId> HFids = HFgeom->getValidDetIds(DetId::Hcal, HcalForward);
0328 
0329   for (i = HFids.begin(); i != HFids.end(); i++) {
0330     nHFdetid++;
0331 
0332     auto cell = HFgeom->getGeometry(*i);
0333     HcalDetId HcalID(i->rawId());
0334     //GlobalPoint p = cell->getPosition();
0335 
0336     int Calo_ieta = 42 + HcalID.ieta();
0337     int Calo_iphi = HcalID.iphi();
0338     double Calo_eta = cell->getPosition().eta();
0339     double Calo_phi = cell->getPosition().phi();
0340 
0341     // HCAL to HF eta, phi map comparison
0342     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0343       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0344       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0345     }
0346 
0347     if (Calo_ieta > HFmax_ieta)
0348       HFmax_ieta = Calo_ieta;
0349     if (Calo_ieta < HFmin_ieta)
0350       HFmin_ieta = Calo_ieta;
0351     if (Calo_iphi > HFmax_iphi)
0352       HFmax_iphi = Calo_iphi;
0353     if (Calo_iphi > HFmax_iphi)
0354       HFmin_iphi = Calo_iphi;
0355   }
0356 
0357   int HOmin_ieta = 99, HOmax_ieta = -99;
0358   int HOmin_iphi = 99, HOmax_iphi = -99;
0359 
0360   // Loop Over all Hcal Outer DetId's
0361   int nHOdetid = 0;
0362   std::vector<DetId> HOids = HOgeom->getValidDetIds(DetId::Hcal, HcalOuter);
0363 
0364   for (i = HOids.begin(); i != HOids.end(); i++) {
0365     nHOdetid++;
0366 
0367     auto cell = HOgeom->getGeometry(*i);
0368     HcalDetId HcalID(i->rawId());
0369     //GlobalPoint p = cell->getPosition();
0370 
0371     int Calo_ieta = 42 + HcalID.ieta();
0372     int Calo_iphi = HcalID.iphi();
0373     double Calo_eta = cell->getPosition().eta();
0374     double Calo_phi = cell->getPosition().phi();
0375 
0376     // HCAL to HO eta, phi map comparison
0377     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0378       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0379       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0380     }
0381 
0382     if (Calo_ieta > HOmax_ieta)
0383       HOmax_ieta = Calo_ieta;
0384     if (Calo_ieta < HOmin_ieta)
0385       HOmin_ieta = Calo_ieta;
0386     if (Calo_iphi > HOmax_iphi)
0387       HOmax_iphi = Calo_iphi;
0388     if (Calo_iphi > HOmax_iphi)
0389       HOmin_iphi = Calo_iphi;
0390   }
0391 
0392   // Set the Cell Size for each (ieta, iphi) Bin
0393   double currentLowEdge_eta = 0;  //double currentHighEdge_eta = 0;
0394   for (int ieta = 1; ieta < 42; ieta++) {
0395     int ieta_ = 42 + ieta;
0396     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(ieta_, 3);
0397     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(ieta_, 3);
0398     double deta = 2.0 * (eta - currentLowEdge_eta);
0399     deta = ((float)((int)(1.0E3 * deta + 0.5))) / 1.0E3;
0400     double dphi = 2.0 * phi;
0401     if (ieta == 40 || ieta == 41)
0402       dphi = 20;
0403     if (ieta <= 39 && ieta >= 21)
0404       dphi = 10;
0405     if (ieta <= 20)
0406       dphi = 5;
0407     // BS: This is WRONG...need to correct overlap
0408     if (ieta == 28)
0409       deta = 0.218;
0410     if (ieta == 29)
0411       deta = 0.096;
0412     currentLowEdge_eta += deta;
0413 
0414     // BS: This is WRONG...need to correct overlap
0415     if (ieta == 29)
0416       currentLowEdge_eta = 2.964;
0417 
0418     hHCAL_ieta_detaMap->setBinContent(ieta_, deta);      // positive rings
0419     hHCAL_ieta_dphiMap->setBinContent(ieta_, dphi);      // positive rings
0420     hHCAL_ieta_detaMap->setBinContent(42 - ieta, deta);  // negative rings
0421     hHCAL_ieta_dphiMap->setBinContent(42 - ieta, dphi);  // negative rings
0422 
0423   }  // end loop over ieta
0424 
0425   edm::LogInfo("OutputInfo") << "HB ieta range: " << HBmin_ieta << " " << HBmax_ieta;
0426   edm::LogInfo("OutputInfo") << "HB iphi range: " << HBmin_iphi << " " << HBmax_iphi;
0427   edm::LogInfo("OutputInfo") << "HE ieta range: " << HEmin_ieta << " " << HEmax_ieta;
0428   edm::LogInfo("OutputInfo") << "HE iphi range: " << HEmin_iphi << " " << HEmax_iphi;
0429   edm::LogInfo("OutputInfo") << "HF ieta range: " << HFmin_ieta << " " << HFmax_ieta;
0430   edm::LogInfo("OutputInfo") << "HF iphi range: " << HFmin_iphi << " " << HFmax_iphi;
0431   edm::LogInfo("OutputInfo") << "HO ieta range: " << HOmin_ieta << " " << HOmax_ieta;
0432   edm::LogInfo("OutputInfo") << "HO iphi range: " << HOmin_iphi << " " << HOmax_iphi;
0433 }
0434 
0435 void HCALRecHitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0436   Nevents++;
0437   hHCAL_Nevents->Fill(0);
0438   // ==========================================================
0439   // Retrieve!
0440   // ==========================================================
0441 
0442   const HBHERecHitCollection* HBHERecHits;
0443   const HORecHitCollection* HORecHits;
0444   const HFRecHitCollection* HFRecHits;
0445 
0446   edm::Handle<HBHERecHitCollection> HBHERecHitsHandle;
0447   iEvent.getByToken(hBHERecHitsLabel_, HBHERecHitsHandle);
0448   if (!HBHERecHitsHandle.isValid()) {
0449     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
0450                                << "HCALRecHitAnalyzer::analyze!\n";
0451     return;
0452   } else {
0453     HBHERecHits = HBHERecHitsHandle.product();
0454   }
0455   edm::Handle<HORecHitCollection> HORecHitsHandle;
0456   iEvent.getByToken(hORecHitsLabel_, HORecHitsHandle);
0457   if (!HORecHitsHandle.isValid()) {
0458     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
0459                                << "HCALRecHitAnalyzer::analyze!\n";
0460     return;
0461   } else {
0462     HORecHits = HORecHitsHandle.product();
0463   }
0464   edm::Handle<HFRecHitCollection> HFRecHitsHandle;
0465   iEvent.getByToken(hFRecHitsLabel_, HFRecHitsHandle);
0466   if (!HFRecHitsHandle.isValid()) {
0467     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
0468                                << "HCALRecHitAnalyzer::analyze!\n";
0469     return;
0470   } else {
0471     HFRecHits = HFRecHitsHandle.product();
0472   }
0473 
0474   // ==========================================================
0475   // Fill Histograms!
0476   // ==========================================================
0477 
0478   TLorentzVector vHBHEMET_EtaRing[83][4];
0479   int HBHEActiveRing[83][4];
0480   int HBHENActiveCells[83][4];
0481   double HBHESET_EtaRing[83][4];
0482   double HBHEMinEnergy_EtaRing[83][4];
0483   double HBHEMaxEnergy_EtaRing[83][4];
0484 
0485   for (int i = 0; i < 83; i++) {
0486     for (int j = 0; j < 4; j++) {
0487       HBHEActiveRing[i][j] = 0;
0488       HBHENActiveCells[i][j] = 0;
0489       HBHESET_EtaRing[i][j] = 0;
0490       HBHEMinEnergy_EtaRing[i][j] = 14E3;
0491       HBHEMaxEnergy_EtaRing[i][j] = -999;
0492     }
0493   }
0494 
0495   // Loop over HBHERecHit's
0496   HBHERecHitCollection::const_iterator hbherechit;
0497   int nHBrechit = 0, nHErechit = 0;
0498 
0499   for (hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); hbherechit++) {
0500     HcalDetId det = hbherechit->id();
0501     double Energy = hbherechit->energy();
0502     Int_t depth = det.depth();
0503     Int_t ieta = det.ieta();
0504     Int_t iphi = det.iphi();
0505     int EtaRing = 41 + ieta;  // this counts from 0
0506     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing + 1, iphi);
0507     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing + 1, iphi);
0508     double theta = 2 * TMath::ATan(exp(-1 * eta));
0509     double ET = Energy * TMath::Sin(theta);
0510     HcalSubdetector HcalNum = det.subdet();
0511     TLorentzVector v_;
0512 
0513     if (Energy > 0)  // zero suppress
0514     {
0515       HBHEActiveRing[EtaRing][depth - 1] = 1;
0516       HBHENActiveCells[EtaRing][depth - 1]++;
0517       HBHESET_EtaRing[EtaRing][depth - 1] += ET;
0518       v_.SetPtEtaPhiE(ET, 0, phi, ET);
0519       vHBHEMET_EtaRing[EtaRing][depth - 1] -= v_;
0520 
0521       DEBUG(EtaRing << " " << Energy << " " << ET << " " << theta << " " << eta << " " << phi << " : "
0522                     << vHBHEMET_EtaRing[EtaRing][depth - 1].Phi() << " " << vHBHEMET_EtaRing[EtaRing][depth - 1].Pt());
0523 
0524       switch (depth) {
0525         case 1:
0526           hHCAL_D1_Occ_ieta_iphi->Fill(ieta, iphi);
0527           break;
0528         case 2:
0529           hHCAL_D2_Occ_ieta_iphi->Fill(ieta, iphi);
0530           break;
0531         case 3:
0532           hHCAL_D3_Occ_ieta_iphi->Fill(ieta, iphi);
0533           break;
0534       }  // end switch
0535     }
0536 
0537     if (Energy > HBHEMaxEnergy_EtaRing[EtaRing][depth - 1])
0538       HBHEMaxEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0539     if (Energy < HBHEMinEnergy_EtaRing[EtaRing][depth - 1])
0540       HBHEMinEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0541 
0542     switch (depth) {
0543       case 1:
0544         hHCAL_D1_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0545         hHCAL_D1_energyvsieta->Fill(ieta, Energy);
0546         if (Energy > hHCAL_D1_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0547           hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0548         if (Energy < hHCAL_D1_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0549           hHCAL_D1_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0550 
0551         break;
0552       case 2:
0553         hHCAL_D2_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0554 
0555         hHCAL_D2_energyvsieta->Fill(ieta, Energy);
0556         if (Energy > hHCAL_D2_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0557           hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0558         if (Energy < hHCAL_D2_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0559           hHCAL_D2_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0560         break;
0561       case 3:
0562         hHCAL_D3_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0563 
0564         hHCAL_D3_energyvsieta->Fill(ieta, Energy);
0565         if (Energy > hHCAL_D3_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0566           hHCAL_D3_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0567         if (Energy < hHCAL_D3_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0568           hHCAL_D3_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0569         break;
0570     }  // end switch
0571 
0572     if (HcalNum == HcalBarrel) {
0573       nHBrechit++;
0574     } else {  // HcalEndcap
0575       nHErechit++;
0576     }
0577   }  // end loop over HBHERecHit's
0578 
0579   // Fill eta-ring MET quantities
0580   for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0581     for (int jDepth = 0; jDepth < 3; jDepth++) {
0582       switch (jDepth + 1) {
0583         case 1:
0584           hHCAL_D1_Maxenergyvsieta->Fill(iEtaRing - 41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
0585           hHCAL_D1_Minenergyvsieta->Fill(iEtaRing - 41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
0586           break;
0587         case 2:
0588           hHCAL_D2_Maxenergyvsieta->Fill(iEtaRing - 41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
0589           hHCAL_D2_Minenergyvsieta->Fill(iEtaRing - 41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
0590           break;
0591         case 3:
0592           hHCAL_D3_Maxenergyvsieta->Fill(iEtaRing - 41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
0593           hHCAL_D3_Minenergyvsieta->Fill(iEtaRing - 41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
0594           break;
0595       }
0596 
0597       if (HBHEActiveRing[iEtaRing][jDepth]) {
0598         switch (jDepth + 1) {
0599           case 1:
0600             hHCAL_D1_METPhivsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
0601             hHCAL_D1_MExvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
0602             hHCAL_D1_MEyvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
0603             hHCAL_D1_METvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
0604             hHCAL_D1_SETvsieta->Fill(iEtaRing - 41, HBHESET_EtaRing[iEtaRing][jDepth]);
0605             hHCAL_D1_Occvsieta->Fill(iEtaRing - 41, HBHENActiveCells[iEtaRing][jDepth]);
0606             break;
0607           case 2:
0608             hHCAL_D2_METPhivsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
0609             hHCAL_D2_MExvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
0610             hHCAL_D2_MEyvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
0611             hHCAL_D2_METvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
0612             hHCAL_D2_SETvsieta->Fill(iEtaRing - 41, HBHESET_EtaRing[iEtaRing][jDepth]);
0613             hHCAL_D2_Occvsieta->Fill(iEtaRing - 41, HBHENActiveCells[iEtaRing][jDepth]);
0614             break;
0615           case 3:
0616             hHCAL_D3_METPhivsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
0617             hHCAL_D3_MExvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
0618             hHCAL_D3_MEyvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
0619             hHCAL_D3_METvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
0620             hHCAL_D3_SETvsieta->Fill(iEtaRing - 41, HBHESET_EtaRing[iEtaRing][jDepth]);
0621             hHCAL_D3_Occvsieta->Fill(iEtaRing - 41, HBHENActiveCells[iEtaRing][jDepth]);
0622             break;
0623         }  //switch
0624       }    //activering
0625     }      //depth
0626   }        //etaring
0627 
0628   //-------------------------------------------------HO
0629   TLorentzVector vHOMET_EtaRing[83];
0630   int HOActiveRing[83];
0631   int HONActiveCells[83];
0632   double HOSET_EtaRing[83];
0633   double HOMinEnergy_EtaRing[83];
0634   double HOMaxEnergy_EtaRing[83];
0635 
0636   for (int i = 0; i < 83; i++) {
0637     HOActiveRing[i] = 0;
0638     HONActiveCells[i] = 0;
0639     HOSET_EtaRing[i] = 0;
0640     HOMinEnergy_EtaRing[i] = 14E3;
0641     HOMaxEnergy_EtaRing[i] = -999;
0642   }
0643 
0644   // Loop over HORecHit's
0645   HORecHitCollection::const_iterator horechit;
0646   int nHOrechit = 0;
0647 
0648   for (horechit = HORecHits->begin(); horechit != HORecHits->end(); horechit++) {
0649     nHOrechit++;
0650 
0651     HcalDetId det = horechit->id();
0652     double Energy = horechit->energy();
0653     ///Int_t depth = det.depth(); //always 4
0654     Int_t ieta = det.ieta();
0655     Int_t iphi = det.iphi();
0656     int EtaRing = 41 + ieta;  // this counts from 0
0657     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing + 1, iphi);
0658     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing + 1, iphi);
0659     double theta = 2 * TMath::ATan(exp(-1 * eta));
0660     double ET = Energy * TMath::Sin(theta);
0661     TLorentzVector v_;
0662 
0663     if (Energy > 0)  // zero suppress
0664     {
0665       HOActiveRing[EtaRing] = 1;
0666       HONActiveCells[EtaRing]++;
0667       HOSET_EtaRing[EtaRing] += ET;
0668       v_.SetPtEtaPhiE(ET, 0, phi, ET);
0669       vHOMET_EtaRing[EtaRing] -= v_;
0670 
0671       hHCAL_D4_Occ_ieta_iphi->Fill(ieta, iphi);
0672     }
0673 
0674     if (Energy > HOMaxEnergy_EtaRing[EtaRing])
0675       HOMaxEnergy_EtaRing[EtaRing] = Energy;
0676     if (Energy < HOMinEnergy_EtaRing[EtaRing])
0677       HOMinEnergy_EtaRing[EtaRing] = Energy;
0678 
0679     hHCAL_D4_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0680 
0681     hHCAL_D4_energyvsieta->Fill(ieta, Energy);
0682     if (Energy > hHCAL_D4_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0683       hHCAL_D4_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0684     if (Energy < hHCAL_D4_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0685       hHCAL_D4_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0686 
0687   }  // end loop over HORecHit's
0688 
0689   // Fill eta-ring MET quantities
0690   for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0691     hHCAL_D4_Maxenergyvsieta->Fill(iEtaRing - 41, HOMaxEnergy_EtaRing[iEtaRing]);
0692     hHCAL_D4_Minenergyvsieta->Fill(iEtaRing - 41, HOMinEnergy_EtaRing[iEtaRing]);
0693 
0694     if (HOActiveRing[iEtaRing]) {
0695       hHCAL_D4_METPhivsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Phi());
0696       hHCAL_D4_MExvsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Px());
0697       hHCAL_D4_MEyvsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Py());
0698       hHCAL_D4_METvsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Pt());
0699       hHCAL_D4_SETvsieta->Fill(iEtaRing - 41, HOSET_EtaRing[iEtaRing]);
0700       hHCAL_D4_Occvsieta->Fill(iEtaRing - 41, HONActiveCells[iEtaRing]);
0701     }
0702   }
0703 
0704   //------------------------------------------------HF
0705 
0706   TLorentzVector vHFMET_EtaRing[83][2];
0707   int HFActiveRing[83][2];
0708   int HFNActiveCells[83][2];
0709   double HFSET_EtaRing[83][2];
0710   double HFMinEnergy_EtaRing[83][2];
0711   double HFMaxEnergy_EtaRing[83][2];
0712 
0713   for (int i = 0; i < 83; i++) {
0714     for (int j = 0; j < 2; j++) {
0715       HFActiveRing[i][j] = 0;
0716       HFNActiveCells[i][j] = 0;
0717       HFSET_EtaRing[i][j] = 0;
0718       HFMinEnergy_EtaRing[i][j] = 14E3;
0719       HFMaxEnergy_EtaRing[i][j] = -999;
0720     }
0721   }
0722 
0723   // Loop over HFRecHit's
0724   HFRecHitCollection::const_iterator hfrechit;
0725   int nHFrechit = 0;
0726 
0727   for (hfrechit = HFRecHits->begin(); hfrechit != HFRecHits->end(); hfrechit++) {
0728     nHFrechit++;
0729 
0730     HcalDetId det = hfrechit->id();
0731     double Energy = hfrechit->energy();
0732     Int_t depth = det.depth();
0733     Int_t ieta = det.ieta();
0734     Int_t iphi = det.iphi();
0735     int EtaRing = 41 + ieta;  // this counts from 0
0736     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing + 1, iphi);
0737     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing + 1, iphi);
0738     double theta = 2 * TMath::ATan(exp(-1 * eta));
0739     double ET = Energy * TMath::Sin(theta);
0740 
0741     TLorentzVector v_;
0742     if (Energy > 0)  // zero suppress
0743     {
0744       HFActiveRing[EtaRing][depth - 1] = 1;
0745       HFNActiveCells[EtaRing][depth - 1]++;
0746       HFSET_EtaRing[EtaRing][depth - 1] += ET;
0747       v_.SetPtEtaPhiE(ET, 0, phi, ET);
0748       vHFMET_EtaRing[EtaRing][depth - 1] -= v_;
0749 
0750       switch (depth) {
0751         case 1:
0752           hHCAL_D1_Occ_ieta_iphi->Fill(ieta, iphi);
0753           break;
0754         case 2:
0755           hHCAL_D2_Occ_ieta_iphi->Fill(ieta, iphi);
0756           break;
0757       }
0758     }
0759 
0760     if (Energy > HFMaxEnergy_EtaRing[EtaRing][depth - 1])
0761       HFMaxEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0762     if (Energy < HFMinEnergy_EtaRing[EtaRing][depth - 1])
0763       HFMinEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0764 
0765     switch (depth) {
0766       case 1:
0767         hHCAL_D1_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0768 
0769         hHCAL_D1_energyvsieta->Fill(ieta, Energy);
0770         if (Energy > hHCAL_D1_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0771           hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0772         if (Energy < hHCAL_D1_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0773           hHCAL_D1_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0774         break;
0775       case 2:
0776         hHCAL_D2_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0777 
0778         hHCAL_D2_energyvsieta->Fill(ieta, Energy);
0779         if (Energy > hHCAL_D2_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0780           hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0781         if (Energy < hHCAL_D2_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0782           hHCAL_D2_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0783         break;
0784     }
0785 
0786   }  // end loop over HFRecHit's
0787 
0788   // Fill eta-ring MET quantities
0789   for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0790     for (int jDepth = 0; jDepth < 2; jDepth++) {
0791       switch (jDepth + 1) {
0792         case 1:
0793           hHCAL_D1_Maxenergyvsieta->Fill(iEtaRing - 41, HFMaxEnergy_EtaRing[iEtaRing][jDepth]);
0794           hHCAL_D1_Minenergyvsieta->Fill(iEtaRing - 41, HFMinEnergy_EtaRing[iEtaRing][jDepth]);
0795           break;
0796         case 2:
0797           hHCAL_D2_Maxenergyvsieta->Fill(iEtaRing - 41, HFMaxEnergy_EtaRing[iEtaRing][jDepth]);
0798           hHCAL_D2_Minenergyvsieta->Fill(iEtaRing - 41, HFMinEnergy_EtaRing[iEtaRing][jDepth]);
0799           break;
0800       }
0801 
0802       if (HFActiveRing[iEtaRing][jDepth]) {
0803         switch (jDepth + 1) {
0804           case 1:
0805 
0806             hHCAL_D1_METPhivsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Phi());
0807             hHCAL_D1_MExvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Px());
0808             hHCAL_D1_MEyvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Py());
0809             hHCAL_D1_METvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Pt());
0810             hHCAL_D1_SETvsieta->Fill(iEtaRing - 41, HFSET_EtaRing[iEtaRing][jDepth]);
0811             hHCAL_D1_Occvsieta->Fill(iEtaRing - 41, HFNActiveCells[iEtaRing][jDepth]);
0812             break;
0813 
0814           case 2:
0815 
0816             hHCAL_D2_METPhivsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Phi());
0817             hHCAL_D2_MExvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Px());
0818             hHCAL_D2_MEyvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Py());
0819             hHCAL_D2_METvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Pt());
0820             hHCAL_D2_SETvsieta->Fill(iEtaRing - 41, HFSET_EtaRing[iEtaRing][jDepth]);
0821             hHCAL_D2_Occvsieta->Fill(iEtaRing - 41, HFNActiveCells[iEtaRing][jDepth]);
0822             break;
0823         }
0824       }
0825     }
0826   }
0827 }