Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-11 04:32:46 

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   std::vector<DetId> HBids = HBgeom->getValidDetIds(DetId::Hcal, HcalBarrel);
0262 
0263   for (i = HBids.begin(); i != HBids.end(); i++) {
0264     HcalDetId HcalID(*i);
0265 
0266     int Calo_ieta = 42 + HcalID.ieta();
0267     int Calo_iphi = HcalID.iphi();
0268     double Calo_eta = HBgeom->getPosition(HcalID).eta();
0269     double Calo_phi = HBgeom->getPosition(HcalID).phi();
0270 
0271     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0272       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0273       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0274     }
0275 
0276     if (Calo_ieta > HBmax_ieta)
0277       HBmax_ieta = Calo_ieta;
0278     if (Calo_ieta < HBmin_ieta)
0279       HBmin_ieta = Calo_ieta;
0280     if (Calo_iphi > HBmax_iphi)
0281       HBmax_iphi = Calo_iphi;
0282     if (Calo_iphi > HBmax_iphi)
0283       HBmin_iphi = Calo_iphi;
0284   }
0285 
0286   int HEmin_ieta = 99, HEmax_ieta = -99;
0287   int HEmin_iphi = 99, HEmax_iphi = -99;
0288 
0289   // Loop Over all Hcal Endcap DetId's
0290   std::vector<DetId> HEids = HEgeom->getValidDetIds(DetId::Hcal, HcalEndcap);
0291 
0292   for (i = HEids.begin(); i != HEids.end(); i++) {
0293     HcalDetId HcalID(*i);
0294 
0295     int Calo_ieta = 42 + HcalID.ieta();
0296     int Calo_iphi = HcalID.iphi();
0297     double Calo_eta = HEgeom->getPosition(HcalID).eta();
0298     double Calo_phi = HEgeom->getPosition(HcalID).phi();
0299 
0300     // HCAL to HE eta, phi map comparison
0301     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0302       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0303       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0304     }
0305 
0306     if (Calo_ieta > HEmax_ieta)
0307       HEmax_ieta = Calo_ieta;
0308     if (Calo_ieta < HEmin_ieta)
0309       HEmin_ieta = Calo_ieta;
0310     if (Calo_iphi > HEmax_iphi)
0311       HEmax_iphi = Calo_iphi;
0312     if (Calo_iphi > HEmax_iphi)
0313       HEmin_iphi = Calo_iphi;
0314   }
0315 
0316   int HFmin_ieta = 99, HFmax_ieta = -99;
0317   int HFmin_iphi = 99, HFmax_iphi = -99;
0318 
0319   // Loop Over all Hcal Forward DetId's
0320   std::vector<DetId> HFids = HFgeom->getValidDetIds(DetId::Hcal, HcalForward);
0321 
0322   for (i = HFids.begin(); i != HFids.end(); i++) {
0323     auto cell = HFgeom->getGeometry(*i);
0324     HcalDetId HcalID(i->rawId());
0325     //GlobalPoint p = cell->getPosition();
0326 
0327     int Calo_ieta = 42 + HcalID.ieta();
0328     int Calo_iphi = HcalID.iphi();
0329     double Calo_eta = cell->getPosition().eta();
0330     double Calo_phi = cell->getPosition().phi();
0331 
0332     // HCAL to HF eta, phi map comparison
0333     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0334       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0335       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0336     }
0337 
0338     if (Calo_ieta > HFmax_ieta)
0339       HFmax_ieta = Calo_ieta;
0340     if (Calo_ieta < HFmin_ieta)
0341       HFmin_ieta = Calo_ieta;
0342     if (Calo_iphi > HFmax_iphi)
0343       HFmax_iphi = Calo_iphi;
0344     if (Calo_iphi > HFmax_iphi)
0345       HFmin_iphi = Calo_iphi;
0346   }
0347 
0348   int HOmin_ieta = 99, HOmax_ieta = -99;
0349   int HOmin_iphi = 99, HOmax_iphi = -99;
0350 
0351   // Loop Over all Hcal Outer DetId's
0352   std::vector<DetId> HOids = HOgeom->getValidDetIds(DetId::Hcal, HcalOuter);
0353 
0354   for (i = HOids.begin(); i != HOids.end(); i++) {
0355     auto cell = HOgeom->getGeometry(*i);
0356     HcalDetId HcalID(i->rawId());
0357     //GlobalPoint p = cell->getPosition();
0358 
0359     int Calo_ieta = 42 + HcalID.ieta();
0360     int Calo_iphi = HcalID.iphi();
0361     double Calo_eta = cell->getPosition().eta();
0362     double Calo_phi = cell->getPosition().phi();
0363 
0364     // HCAL to HO eta, phi map comparison
0365     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta, Calo_iphi) == -999) {
0366       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta, Calo_iphi, Calo_eta);
0367       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta, Calo_iphi, Calo_phi * 180.0 / M_PI);
0368     }
0369 
0370     if (Calo_ieta > HOmax_ieta)
0371       HOmax_ieta = Calo_ieta;
0372     if (Calo_ieta < HOmin_ieta)
0373       HOmin_ieta = Calo_ieta;
0374     if (Calo_iphi > HOmax_iphi)
0375       HOmax_iphi = Calo_iphi;
0376     if (Calo_iphi > HOmax_iphi)
0377       HOmin_iphi = Calo_iphi;
0378   }
0379 
0380   // Set the Cell Size for each (ieta, iphi) Bin
0381   double currentLowEdge_eta = 0;  //double currentHighEdge_eta = 0;
0382   for (int ieta = 1; ieta < 42; ieta++) {
0383     int ieta_ = 42 + ieta;
0384     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(ieta_, 3);
0385     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(ieta_, 3);
0386     double deta = 2.0 * (eta - currentLowEdge_eta);
0387     deta = ((float)((int)(1.0E3 * deta + 0.5))) / 1.0E3;
0388     double dphi = 2.0 * phi;
0389     if (ieta == 40 || ieta == 41)
0390       dphi = 20;
0391     if (ieta <= 39 && ieta >= 21)
0392       dphi = 10;
0393     if (ieta <= 20)
0394       dphi = 5;
0395     // BS: This is WRONG...need to correct overlap
0396     if (ieta == 28)
0397       deta = 0.218;
0398     if (ieta == 29)
0399       deta = 0.096;
0400     currentLowEdge_eta += deta;
0401 
0402     // BS: This is WRONG...need to correct overlap
0403     if (ieta == 29)
0404       currentLowEdge_eta = 2.964;
0405 
0406     hHCAL_ieta_detaMap->setBinContent(ieta_, deta);      // positive rings
0407     hHCAL_ieta_dphiMap->setBinContent(ieta_, dphi);      // positive rings
0408     hHCAL_ieta_detaMap->setBinContent(42 - ieta, deta);  // negative rings
0409     hHCAL_ieta_dphiMap->setBinContent(42 - ieta, dphi);  // negative rings
0410 
0411   }  // end loop over ieta
0412 
0413   edm::LogInfo("OutputInfo") << "HB ieta range: " << HBmin_ieta << " " << HBmax_ieta;
0414   edm::LogInfo("OutputInfo") << "HB iphi range: " << HBmin_iphi << " " << HBmax_iphi;
0415   edm::LogInfo("OutputInfo") << "HE ieta range: " << HEmin_ieta << " " << HEmax_ieta;
0416   edm::LogInfo("OutputInfo") << "HE iphi range: " << HEmin_iphi << " " << HEmax_iphi;
0417   edm::LogInfo("OutputInfo") << "HF ieta range: " << HFmin_ieta << " " << HFmax_ieta;
0418   edm::LogInfo("OutputInfo") << "HF iphi range: " << HFmin_iphi << " " << HFmax_iphi;
0419   edm::LogInfo("OutputInfo") << "HO ieta range: " << HOmin_ieta << " " << HOmax_ieta;
0420   edm::LogInfo("OutputInfo") << "HO iphi range: " << HOmin_iphi << " " << HOmax_iphi;
0421 }
0422 
0423 void HCALRecHitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0424   Nevents++;
0425   hHCAL_Nevents->Fill(0);
0426   // ==========================================================
0427   // Retrieve!
0428   // ==========================================================
0429 
0430   const HBHERecHitCollection* HBHERecHits;
0431   const HORecHitCollection* HORecHits;
0432   const HFRecHitCollection* HFRecHits;
0433 
0434   edm::Handle<HBHERecHitCollection> HBHERecHitsHandle;
0435   iEvent.getByToken(hBHERecHitsLabel_, HBHERecHitsHandle);
0436   if (!HBHERecHitsHandle.isValid()) {
0437     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
0438                                << "HCALRecHitAnalyzer::analyze!\n";
0439     return;
0440   } else {
0441     HBHERecHits = HBHERecHitsHandle.product();
0442   }
0443   edm::Handle<HORecHitCollection> HORecHitsHandle;
0444   iEvent.getByToken(hORecHitsLabel_, HORecHitsHandle);
0445   if (!HORecHitsHandle.isValid()) {
0446     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
0447                                << "HCALRecHitAnalyzer::analyze!\n";
0448     return;
0449   } else {
0450     HORecHits = HORecHitsHandle.product();
0451   }
0452   edm::Handle<HFRecHitCollection> HFRecHitsHandle;
0453   iEvent.getByToken(hFRecHitsLabel_, HFRecHitsHandle);
0454   if (!HFRecHitsHandle.isValid()) {
0455     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
0456                                << "HCALRecHitAnalyzer::analyze!\n";
0457     return;
0458   } else {
0459     HFRecHits = HFRecHitsHandle.product();
0460   }
0461 
0462   // ==========================================================
0463   // Fill Histograms!
0464   // ==========================================================
0465 
0466   TLorentzVector vHBHEMET_EtaRing[83][4];
0467   int HBHEActiveRing[83][4];
0468   int HBHENActiveCells[83][4];
0469   double HBHESET_EtaRing[83][4];
0470   double HBHEMinEnergy_EtaRing[83][4];
0471   double HBHEMaxEnergy_EtaRing[83][4];
0472 
0473   for (int i = 0; i < 83; i++) {
0474     for (int j = 0; j < 4; j++) {
0475       HBHEActiveRing[i][j] = 0;
0476       HBHENActiveCells[i][j] = 0;
0477       HBHESET_EtaRing[i][j] = 0;
0478       HBHEMinEnergy_EtaRing[i][j] = 14E3;
0479       HBHEMaxEnergy_EtaRing[i][j] = -999;
0480     }
0481   }
0482 
0483   // Loop over HBHERecHit's
0484   HBHERecHitCollection::const_iterator hbherechit;
0485 
0486   for (hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); hbherechit++) {
0487     HcalDetId det = hbherechit->id();
0488     double Energy = hbherechit->energy();
0489     Int_t depth = det.depth();
0490     Int_t ieta = det.ieta();
0491     Int_t iphi = det.iphi();
0492     int EtaRing = 41 + ieta;  // this counts from 0
0493     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing + 1, iphi);
0494     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing + 1, iphi);
0495     double theta = 2 * TMath::ATan(exp(-1 * eta));
0496     double ET = Energy * TMath::Sin(theta);
0497     TLorentzVector v_;
0498 
0499     if (Energy > 0)  // zero suppress
0500     {
0501       HBHEActiveRing[EtaRing][depth - 1] = 1;
0502       HBHENActiveCells[EtaRing][depth - 1]++;
0503       HBHESET_EtaRing[EtaRing][depth - 1] += ET;
0504       v_.SetPtEtaPhiE(ET, 0, phi, ET);
0505       vHBHEMET_EtaRing[EtaRing][depth - 1] -= v_;
0506 
0507       DEBUG(EtaRing << " " << Energy << " " << ET << " " << theta << " " << eta << " " << phi << " : "
0508                     << vHBHEMET_EtaRing[EtaRing][depth - 1].Phi() << " " << vHBHEMET_EtaRing[EtaRing][depth - 1].Pt());
0509 
0510       switch (depth) {
0511         case 1:
0512           hHCAL_D1_Occ_ieta_iphi->Fill(ieta, iphi);
0513           break;
0514         case 2:
0515           hHCAL_D2_Occ_ieta_iphi->Fill(ieta, iphi);
0516           break;
0517         case 3:
0518           hHCAL_D3_Occ_ieta_iphi->Fill(ieta, iphi);
0519           break;
0520       }  // end switch
0521     }
0522 
0523     if (Energy > HBHEMaxEnergy_EtaRing[EtaRing][depth - 1])
0524       HBHEMaxEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0525     if (Energy < HBHEMinEnergy_EtaRing[EtaRing][depth - 1])
0526       HBHEMinEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0527 
0528     switch (depth) {
0529       case 1:
0530         hHCAL_D1_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0531         hHCAL_D1_energyvsieta->Fill(ieta, Energy);
0532         if (Energy > hHCAL_D1_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0533           hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0534         if (Energy < hHCAL_D1_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0535           hHCAL_D1_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0536 
0537         break;
0538       case 2:
0539         hHCAL_D2_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0540 
0541         hHCAL_D2_energyvsieta->Fill(ieta, Energy);
0542         if (Energy > hHCAL_D2_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0543           hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0544         if (Energy < hHCAL_D2_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0545           hHCAL_D2_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0546         break;
0547       case 3:
0548         hHCAL_D3_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0549 
0550         hHCAL_D3_energyvsieta->Fill(ieta, Energy);
0551         if (Energy > hHCAL_D3_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0552           hHCAL_D3_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0553         if (Energy < hHCAL_D3_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0554           hHCAL_D3_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0555         break;
0556     }  // end switch
0557   }  // end loop over HBHERecHit's
0558 
0559   // Fill eta-ring MET quantities
0560   for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0561     for (int jDepth = 0; jDepth < 3; jDepth++) {
0562       switch (jDepth + 1) {
0563         case 1:
0564           hHCAL_D1_Maxenergyvsieta->Fill(iEtaRing - 41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
0565           hHCAL_D1_Minenergyvsieta->Fill(iEtaRing - 41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
0566           break;
0567         case 2:
0568           hHCAL_D2_Maxenergyvsieta->Fill(iEtaRing - 41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
0569           hHCAL_D2_Minenergyvsieta->Fill(iEtaRing - 41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
0570           break;
0571         case 3:
0572           hHCAL_D3_Maxenergyvsieta->Fill(iEtaRing - 41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
0573           hHCAL_D3_Minenergyvsieta->Fill(iEtaRing - 41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
0574           break;
0575       }
0576 
0577       if (HBHEActiveRing[iEtaRing][jDepth]) {
0578         switch (jDepth + 1) {
0579           case 1:
0580             hHCAL_D1_METPhivsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
0581             hHCAL_D1_MExvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
0582             hHCAL_D1_MEyvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
0583             hHCAL_D1_METvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
0584             hHCAL_D1_SETvsieta->Fill(iEtaRing - 41, HBHESET_EtaRing[iEtaRing][jDepth]);
0585             hHCAL_D1_Occvsieta->Fill(iEtaRing - 41, HBHENActiveCells[iEtaRing][jDepth]);
0586             break;
0587           case 2:
0588             hHCAL_D2_METPhivsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
0589             hHCAL_D2_MExvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
0590             hHCAL_D2_MEyvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
0591             hHCAL_D2_METvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
0592             hHCAL_D2_SETvsieta->Fill(iEtaRing - 41, HBHESET_EtaRing[iEtaRing][jDepth]);
0593             hHCAL_D2_Occvsieta->Fill(iEtaRing - 41, HBHENActiveCells[iEtaRing][jDepth]);
0594             break;
0595           case 3:
0596             hHCAL_D3_METPhivsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
0597             hHCAL_D3_MExvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
0598             hHCAL_D3_MEyvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
0599             hHCAL_D3_METvsieta->Fill(iEtaRing - 41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
0600             hHCAL_D3_SETvsieta->Fill(iEtaRing - 41, HBHESET_EtaRing[iEtaRing][jDepth]);
0601             hHCAL_D3_Occvsieta->Fill(iEtaRing - 41, HBHENActiveCells[iEtaRing][jDepth]);
0602             break;
0603         }  //switch
0604       }  //activering
0605     }  //depth
0606   }  //etaring
0607 
0608   //-------------------------------------------------HO
0609   TLorentzVector vHOMET_EtaRing[83];
0610   int HOActiveRing[83];
0611   int HONActiveCells[83];
0612   double HOSET_EtaRing[83];
0613   double HOMinEnergy_EtaRing[83];
0614   double HOMaxEnergy_EtaRing[83];
0615 
0616   for (int i = 0; i < 83; i++) {
0617     HOActiveRing[i] = 0;
0618     HONActiveCells[i] = 0;
0619     HOSET_EtaRing[i] = 0;
0620     HOMinEnergy_EtaRing[i] = 14E3;
0621     HOMaxEnergy_EtaRing[i] = -999;
0622   }
0623 
0624   // Loop over HORecHit's
0625   HORecHitCollection::const_iterator horechit;
0626 
0627   for (horechit = HORecHits->begin(); horechit != HORecHits->end(); horechit++) {
0628     HcalDetId det = horechit->id();
0629     double Energy = horechit->energy();
0630     ///Int_t depth = det.depth(); //always 4
0631     Int_t ieta = det.ieta();
0632     Int_t iphi = det.iphi();
0633     int EtaRing = 41 + ieta;  // this counts from 0
0634     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing + 1, iphi);
0635     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing + 1, iphi);
0636     double theta = 2 * TMath::ATan(exp(-1 * eta));
0637     double ET = Energy * TMath::Sin(theta);
0638     TLorentzVector v_;
0639 
0640     if (Energy > 0)  // zero suppress
0641     {
0642       HOActiveRing[EtaRing] = 1;
0643       HONActiveCells[EtaRing]++;
0644       HOSET_EtaRing[EtaRing] += ET;
0645       v_.SetPtEtaPhiE(ET, 0, phi, ET);
0646       vHOMET_EtaRing[EtaRing] -= v_;
0647 
0648       hHCAL_D4_Occ_ieta_iphi->Fill(ieta, iphi);
0649     }
0650 
0651     if (Energy > HOMaxEnergy_EtaRing[EtaRing])
0652       HOMaxEnergy_EtaRing[EtaRing] = Energy;
0653     if (Energy < HOMinEnergy_EtaRing[EtaRing])
0654       HOMinEnergy_EtaRing[EtaRing] = Energy;
0655 
0656     hHCAL_D4_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0657 
0658     hHCAL_D4_energyvsieta->Fill(ieta, Energy);
0659     if (Energy > hHCAL_D4_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0660       hHCAL_D4_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0661     if (Energy < hHCAL_D4_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0662       hHCAL_D4_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0663 
0664   }  // end loop over HORecHit's
0665 
0666   // Fill eta-ring MET quantities
0667   for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0668     hHCAL_D4_Maxenergyvsieta->Fill(iEtaRing - 41, HOMaxEnergy_EtaRing[iEtaRing]);
0669     hHCAL_D4_Minenergyvsieta->Fill(iEtaRing - 41, HOMinEnergy_EtaRing[iEtaRing]);
0670 
0671     if (HOActiveRing[iEtaRing]) {
0672       hHCAL_D4_METPhivsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Phi());
0673       hHCAL_D4_MExvsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Px());
0674       hHCAL_D4_MEyvsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Py());
0675       hHCAL_D4_METvsieta->Fill(iEtaRing - 41, vHOMET_EtaRing[iEtaRing].Pt());
0676       hHCAL_D4_SETvsieta->Fill(iEtaRing - 41, HOSET_EtaRing[iEtaRing]);
0677       hHCAL_D4_Occvsieta->Fill(iEtaRing - 41, HONActiveCells[iEtaRing]);
0678     }
0679   }
0680 
0681   //------------------------------------------------HF
0682 
0683   TLorentzVector vHFMET_EtaRing[83][2];
0684   int HFActiveRing[83][2];
0685   int HFNActiveCells[83][2];
0686   double HFSET_EtaRing[83][2];
0687   double HFMinEnergy_EtaRing[83][2];
0688   double HFMaxEnergy_EtaRing[83][2];
0689 
0690   for (int i = 0; i < 83; i++) {
0691     for (int j = 0; j < 2; j++) {
0692       HFActiveRing[i][j] = 0;
0693       HFNActiveCells[i][j] = 0;
0694       HFSET_EtaRing[i][j] = 0;
0695       HFMinEnergy_EtaRing[i][j] = 14E3;
0696       HFMaxEnergy_EtaRing[i][j] = -999;
0697     }
0698   }
0699 
0700   // Loop over HFRecHit's
0701   HFRecHitCollection::const_iterator hfrechit;
0702 
0703   for (hfrechit = HFRecHits->begin(); hfrechit != HFRecHits->end(); hfrechit++) {
0704     HcalDetId det = hfrechit->id();
0705     double Energy = hfrechit->energy();
0706     Int_t depth = det.depth();
0707     Int_t ieta = det.ieta();
0708     Int_t iphi = det.iphi();
0709     int EtaRing = 41 + ieta;  // this counts from 0
0710     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing + 1, iphi);
0711     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing + 1, iphi);
0712     double theta = 2 * TMath::ATan(exp(-1 * eta));
0713     double ET = Energy * TMath::Sin(theta);
0714 
0715     TLorentzVector v_;
0716     if (Energy > 0)  // zero suppress
0717     {
0718       HFActiveRing[EtaRing][depth - 1] = 1;
0719       HFNActiveCells[EtaRing][depth - 1]++;
0720       HFSET_EtaRing[EtaRing][depth - 1] += ET;
0721       v_.SetPtEtaPhiE(ET, 0, phi, ET);
0722       vHFMET_EtaRing[EtaRing][depth - 1] -= v_;
0723 
0724       switch (depth) {
0725         case 1:
0726           hHCAL_D1_Occ_ieta_iphi->Fill(ieta, iphi);
0727           break;
0728         case 2:
0729           hHCAL_D2_Occ_ieta_iphi->Fill(ieta, iphi);
0730           break;
0731       }
0732     }
0733 
0734     if (Energy > HFMaxEnergy_EtaRing[EtaRing][depth - 1])
0735       HFMaxEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0736     if (Energy < HFMinEnergy_EtaRing[EtaRing][depth - 1])
0737       HFMinEnergy_EtaRing[EtaRing][depth - 1] = Energy;
0738 
0739     switch (depth) {
0740       case 1:
0741         hHCAL_D1_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0742 
0743         hHCAL_D1_energyvsieta->Fill(ieta, Energy);
0744         if (Energy > hHCAL_D1_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0745           hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0746         if (Energy < hHCAL_D1_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0747           hHCAL_D1_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0748         break;
0749       case 2:
0750         hHCAL_D2_energy_ieta_iphi->Fill(ieta, iphi, Energy);
0751 
0752         hHCAL_D2_energyvsieta->Fill(ieta, Energy);
0753         if (Energy > hHCAL_D2_Maxenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0754           hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0755         if (Energy < hHCAL_D2_Minenergy_ieta_iphi->getBinContent(EtaRing + 1, iphi))
0756           hHCAL_D2_Minenergy_ieta_iphi->setBinContent(EtaRing + 1, iphi, Energy);
0757         break;
0758     }
0759 
0760   }  // end loop over HFRecHit's
0761 
0762   // Fill eta-ring MET quantities
0763   for (int iEtaRing = 0; iEtaRing < 83; iEtaRing++) {
0764     for (int jDepth = 0; jDepth < 2; jDepth++) {
0765       switch (jDepth + 1) {
0766         case 1:
0767           hHCAL_D1_Maxenergyvsieta->Fill(iEtaRing - 41, HFMaxEnergy_EtaRing[iEtaRing][jDepth]);
0768           hHCAL_D1_Minenergyvsieta->Fill(iEtaRing - 41, HFMinEnergy_EtaRing[iEtaRing][jDepth]);
0769           break;
0770         case 2:
0771           hHCAL_D2_Maxenergyvsieta->Fill(iEtaRing - 41, HFMaxEnergy_EtaRing[iEtaRing][jDepth]);
0772           hHCAL_D2_Minenergyvsieta->Fill(iEtaRing - 41, HFMinEnergy_EtaRing[iEtaRing][jDepth]);
0773           break;
0774       }
0775 
0776       if (HFActiveRing[iEtaRing][jDepth]) {
0777         switch (jDepth + 1) {
0778           case 1:
0779 
0780             hHCAL_D1_METPhivsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Phi());
0781             hHCAL_D1_MExvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Px());
0782             hHCAL_D1_MEyvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Py());
0783             hHCAL_D1_METvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Pt());
0784             hHCAL_D1_SETvsieta->Fill(iEtaRing - 41, HFSET_EtaRing[iEtaRing][jDepth]);
0785             hHCAL_D1_Occvsieta->Fill(iEtaRing - 41, HFNActiveCells[iEtaRing][jDepth]);
0786             break;
0787 
0788           case 2:
0789 
0790             hHCAL_D2_METPhivsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Phi());
0791             hHCAL_D2_MExvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Px());
0792             hHCAL_D2_MEyvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Py());
0793             hHCAL_D2_METvsieta->Fill(iEtaRing - 41, vHFMET_EtaRing[iEtaRing][jDepth].Pt());
0794             hHCAL_D2_SETvsieta->Fill(iEtaRing - 41, HFSET_EtaRing[iEtaRing][jDepth]);
0795             hHCAL_D2_Occvsieta->Fill(iEtaRing - 41, HFNActiveCells[iEtaRing][jDepth]);
0796             break;
0797         }
0798       }
0799     }
0800   }
0801 }
[ 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!