Reference Manual

Go to the documentation of this file.
00001 #include "DQMOffline/JetMET/interface/HCALRecHitAnalyzer.h"
00002 // author: Bobby Scurlock, University of Florida
00003 // first version 12/7/2006
00004 // modified: Mike Schmitt
00005 // date: 03.05.2007
00006 // note: 1) code rewrite. 2.) changed to loop over all hcal detids;
00007 //       not only those within calotowers.
00008 
00009 #include "FWCore/Framework/interface/MakerMacros.h"
00010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00011 #include "FWCore/ServiceRegistry/interface/Service.h"
00012                                                                            
00013 //#include "PluginManager/ModuleDef.h"
00014 #include "FWCore/PluginManager/interface/ModuleDef.h"
00015 #include "FWCore/Framework/interface/Event.h"
00016 #include "FWCore/Framework/interface/ESHandle.h"
00017 #include "DataFormats/Common/interface/Handle.h"
00018 //#include "FWCore/Framework/interface/Handle.h"
00019 #include "FWCore/Framework/interface/EventSetup.h"
00020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00021 
00022 //#include "Geometry/Vector/interface/GlobalPoint.h"
00023 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00024 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
00025 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
00026 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
00027 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00028 
00029 #include "DataFormats/DetId/interface/DetId.h"
00030 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00031 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00032 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
00033 #include "DataFormats/HcalRecHit/interface/HFRecHit.h"
00034 #include "DataFormats/HcalRecHit/interface/HORecHit.h"
00035 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00036 
00037 #include <memory>
00038 #include <vector>
00039 #include <utility>
00040 #include <ostream>
00041 #include <fstream>
00042 #include <string>
00043 #include <algorithm>
00044 #include <cmath>
00045 #include <TLorentzVector.h>
00046 #include "DQMServices/Core/interface/DQMStore.h"
00047 
00048 #define DEBUG(X) { if (debug_) { std::cout << X << std::endl; } }
00049 
00050 HCALRecHitAnalyzer::HCALRecHitAnalyzer(const edm::ParameterSet& iConfig)
00051 {
00052 
00053   // Retrieve Information from the Configuration File
00054   hBHERecHitsLabel_  = iConfig.getParameter<edm::InputTag>("HBHERecHitsLabel");
00055   hORecHitsLabel_    = iConfig.getParameter<edm::InputTag>("HORecHitsLabel");
00056   hFRecHitsLabel_    = iConfig.getParameter<edm::InputTag>("HFRecHitsLabel");
00057 
00058   debug_             = iConfig.getParameter<bool>("Debug");
00059   finebinning_ = iConfig.getUntrackedParameter<bool>("FineBinning");
00060   FolderName_ = iConfig.getUntrackedParameter<std::string>("FolderName");
00061 
00062 }
00063 
00064 void HCALRecHitAnalyzer::beginRun(const edm::Run& iRun,  const edm::EventSetup& iSetup)
00065 {
00066   Nevents = 0;
00067   BookHistos();
00068   FillGeometry(iSetup);
00069 
00070 }
00071 
00072 void HCALRecHitAnalyzer::BookHistos(){
00073   // get ahold of back-end interface
00074   dbe_ = edm::Service<DQMStore>().operator->();
00075 
00076 
00077   if (dbe_) {
00078 
00079     dbe_->setCurrentFolder(FolderName_+"/geometry");
00080     hHCAL_ieta_iphi_HBMap = dbe_->book2D("METTask_HCAL_ieta_iphi_HBMap","",83,-41,42,72,1,73); 
00081     hHCAL_ieta_iphi_HEMap = dbe_->book2D("METTask_HCAL_ieta_iphi_HEMap","",83,-41,42,72,1,73); 
00082     hHCAL_ieta_iphi_HFMap = dbe_->book2D("METTask_HCAL_ieta_iphi_HFMap","",83,-41,42,72,1,73); 
00083     hHCAL_ieta_iphi_HOMap = dbe_->book2D("METTask_HCAL_ieta_iphi_HOMap","",83,-41,42,72,1,73); 
00084     hHCAL_ieta_iphi_etaMap = dbe_->book2D("METTask_HCAL_ieta_iphi_etaMap","",83,-41,42,72,1,73);
00085     hHCAL_ieta_iphi_phiMap = dbe_->book2D("METTask_HCAL_ieta_iphi_phiMap","",83,-41,42,72,1,73);
00086     hHCAL_ieta_detaMap = dbe_->book1D("METTask_HCAL_ieta_detaMap","",83,-41,42);
00087     hHCAL_ieta_dphiMap = dbe_->book1D("METTask_HCAL_ieta_dphiMap","",83,-41,42);  
00088 
00089     // Initialize bins for geometry to -999 because z = 0 is a valid entry 
00090     for (int i = 1; i <= 83; i++) {
00091 
00092       hHCAL_ieta_detaMap->setBinContent(i,-999);
00093       hHCAL_ieta_dphiMap->setBinContent(i,-999);
00094 
00095       for (int j = 1; j <= 72; j++) {
00096 
00097         hHCAL_ieta_iphi_HBMap->setBinContent(i,j,0);
00098         hHCAL_ieta_iphi_HEMap->setBinContent(i,j,0);
00099         hHCAL_ieta_iphi_HFMap->setBinContent(i,j,0);
00100         hHCAL_ieta_iphi_HOMap->setBinContent(i,j,0);
00101         hHCAL_ieta_iphi_etaMap->setBinContent(i,j,-999);
00102         hHCAL_ieta_iphi_phiMap->setBinContent(i,j,-999);
00103 
00104       }
00105 
00106     }
00107   
00108     //    dbe_->setCurrentFolder("RecoMETV/MET_HCAL/data");
00109     dbe_->setCurrentFolder(FolderName_);
00110     //--Store number of events used
00111     hHCAL_Nevents          = dbe_->book1D("METTask_HCAL_Nevents","",1,0,1);  
00112     //--Data integrated over all events and stored by HCAL(ieta,iphi) 
00113 
00114     hHCAL_D1_energy_ieta_iphi = dbe_->book2D("METTask_HCAL_D1_energy_ieta_iphi","",83,-41,42,72,1,73);  
00115     hHCAL_D2_energy_ieta_iphi = dbe_->book2D("METTask_HCAL_D2_energy_ieta_iphi","",83,-41,42,72,1,73);  
00116     hHCAL_D3_energy_ieta_iphi = dbe_->book2D("METTask_HCAL_D3_energy_ieta_iphi","",83,-41,42,72,1,73);  
00117     hHCAL_D4_energy_ieta_iphi = dbe_->book2D("METTask_HCAL_D4_energy_ieta_iphi","",83,-41,42,72,1,73);  
00118     
00119     hHCAL_D1_Minenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D1_Minenergy_ieta_iphi","",83,-41,42,72,1,73);  
00120     hHCAL_D2_Minenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D2_Minenergy_ieta_iphi","",83,-41,42,72,1,73);  
00121     hHCAL_D3_Minenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D3_Minenergy_ieta_iphi","",83,-41,42,72,1,73);  
00122     hHCAL_D4_Minenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D4_Minenergy_ieta_iphi","",83,-41,42,72,1,73); 
00123     
00124     hHCAL_D1_Maxenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D1_Maxenergy_ieta_iphi","",83,-41,42,72,1,73);  
00125     hHCAL_D2_Maxenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D2_Maxenergy_ieta_iphi","",83,-41,42,72,1,73);  
00126     hHCAL_D3_Maxenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D3_Maxenergy_ieta_iphi","",83,-41,42,72,1,73);  
00127     hHCAL_D4_Maxenergy_ieta_iphi = dbe_->book2D("METTask_HCAL_D4_Maxenergy_ieta_iphi","",83,-41,42,72,1,73); 
00128     
00129     // need to initialize those
00130     for (int i=1; i<=83; i++)
00131       for (int j=1; j<=73; j++)
00132         {
00133           hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(i,j,-999);
00134           hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(i,j,-999);
00135           hHCAL_D3_Maxenergy_ieta_iphi->setBinContent(i,j,-999);
00136           hHCAL_D4_Maxenergy_ieta_iphi->setBinContent(i,j,-999);
00137 
00138           hHCAL_D1_Minenergy_ieta_iphi->setBinContent(i,j, 14000);
00139           hHCAL_D2_Minenergy_ieta_iphi->setBinContent(i,j, 14000);
00140           hHCAL_D3_Minenergy_ieta_iphi->setBinContent(i,j, 14000);
00141           hHCAL_D4_Minenergy_ieta_iphi->setBinContent(i,j, 14000);
00142         }
00143 
00144     hHCAL_D1_Occ_ieta_iphi = dbe_->book2D("METTask_HCAL_D1_Occ_ieta_iphi","",83,-41,42,72,1,73);  
00145     hHCAL_D2_Occ_ieta_iphi = dbe_->book2D("METTask_HCAL_D2_Occ_ieta_iphi","",83,-41,42,72,1,73);  
00146     hHCAL_D3_Occ_ieta_iphi = dbe_->book2D("METTask_HCAL_D3_Occ_ieta_iphi","",83,-41,42,72,1,73);  
00147     hHCAL_D4_Occ_ieta_iphi = dbe_->book2D("METTask_HCAL_D4_Occ_ieta_iphi","",83,-41,42,72,1,73);  
00148     //--Data over eta-rings
00149     
00150     
00151     // CaloTower values
00152 
00153     if(finebinning_)
00154       {
00155         hHCAL_D1_energyvsieta = dbe_->book2D("METTask_HCAL_D1_energyvsieta","",83,-41,42,20101,-100,2001);  
00156         hHCAL_D2_energyvsieta = dbe_->book2D("METTask_HCAL_D2_energyvsieta","",83,-41,42,20101,-100,2001);
00157         hHCAL_D3_energyvsieta = dbe_->book2D("METTask_HCAL_D3_energyvsieta","",83,-41,42,20101,-100,2001);
00158         hHCAL_D4_energyvsieta = dbe_->book2D("METTask_HCAL_D4_energyvsieta","",83,-41,42,20101,-100,2001);
00159         
00160         hHCAL_D1_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D1_Minenergyvsieta","",83,-41,42,20101,-100,2001);  
00161         hHCAL_D2_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D2_Minenergyvsieta","",83,-41,42,20101,-100,2001);
00162         hHCAL_D3_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D3_Minenergyvsieta","",83,-41,42,20101,-100,2001);
00163         hHCAL_D4_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D4_Minenergyvsieta","",83,-41,42,20101,-100,2001);
00164         
00165         hHCAL_D1_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D1_Maxenergyvsieta","",83,-41,42,20101,-100,2001);  
00166         hHCAL_D2_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D2_Maxenergyvsieta","",83,-41,42,20101,-100,2001);
00167         hHCAL_D3_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D3_Maxenergyvsieta","",83,-41,42,20101,-100,2001);
00168         hHCAL_D4_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D4_Maxenergyvsieta","",83,-41,42,20101,-100,2001);
00169         
00170         // Integrated over phi
00171         hHCAL_D1_Occvsieta = dbe_->book2D("METTask_HCAL_D1_Occvsieta","",83,-41,42,73,0,73);
00172         hHCAL_D2_Occvsieta = dbe_->book2D("METTask_HCAL_D2_Occvsieta","",83,-41,42,73,0,73);
00173         hHCAL_D3_Occvsieta = dbe_->book2D("METTask_HCAL_D3_Occvsieta","",83,-41,42,73,0,73);
00174         hHCAL_D4_Occvsieta = dbe_->book2D("METTask_HCAL_D4_Occvsieta","",83,-41,42,73,0,73);
00175         
00176         hHCAL_D1_SETvsieta = dbe_->book2D("METTask_HCAL_D1_SETvsieta","",83,-41,42,20001,0,2001); 
00177         hHCAL_D2_SETvsieta = dbe_->book2D("METTask_HCAL_D2_SETvsieta","",83,-41,42,20001,0,2001); 
00178         hHCAL_D3_SETvsieta = dbe_->book2D("METTask_HCAL_D3_SETvsieta","",83,-41,42,20001,0,2001); 
00179         hHCAL_D4_SETvsieta = dbe_->book2D("METTask_HCAL_D4_SETvsieta","",83,-41,42,20001,0,2001); 
00180         
00181         hHCAL_D1_METvsieta = dbe_->book2D("METTask_HCAL_D1_METvsieta","",83,-41,42,20001,0,2001); 
00182         hHCAL_D2_METvsieta = dbe_->book2D("METTask_HCAL_D2_METvsieta","",83,-41,42,20001,0,2001); 
00183         hHCAL_D3_METvsieta = dbe_->book2D("METTask_HCAL_D3_METvsieta","",83,-41,42,20001,0,2001); 
00184         hHCAL_D4_METvsieta = dbe_->book2D("METTask_HCAL_D4_METvsieta","",83,-41,42,20001,0,2001); 
00185         
00186         hHCAL_D1_METPhivsieta      = dbe_->book2D("METTask_HCAL_D1_METPhivsieta","",83,-41,42, 80, -4,4);  
00187         hHCAL_D2_METPhivsieta      = dbe_->book2D("METTask_HCAL_D2_METPhivsieta","",83,-41,42, 80, -4,4);  
00188         hHCAL_D3_METPhivsieta      = dbe_->book2D("METTask_HCAL_D3_METPhivsieta","",83,-41,42, 80, -4,4);  
00189         hHCAL_D4_METPhivsieta      = dbe_->book2D("METTask_HCAL_D4_METPhivsieta","",83,-41,42, 80, -4,4);  
00190         
00191         hHCAL_D1_MExvsieta         = dbe_->book2D("METTask_HCAL_D1_MExvsieta","",83,-41,42, 10001,-500,501);  
00192         hHCAL_D2_MExvsieta         = dbe_->book2D("METTask_HCAL_D2_MExvsieta","",83,-41,42, 10001,-500,501);  
00193         hHCAL_D3_MExvsieta         = dbe_->book2D("METTask_HCAL_D3_MExvsieta","",83,-41,42, 10001,-500,501);  
00194         hHCAL_D4_MExvsieta         = dbe_->book2D("METTask_HCAL_D4_MExvsieta","",83,-41,42, 10001,-500,501);  
00195         
00196         hHCAL_D1_MEyvsieta         = dbe_->book2D("METTask_HCAL_D1_MEyvsieta","",83,-41,42, 10001,-500,501);  
00197         hHCAL_D2_MEyvsieta         = dbe_->book2D("METTask_HCAL_D2_MEyvsieta","",83,-41,42, 10001,-500,501);  
00198         hHCAL_D3_MEyvsieta         = dbe_->book2D("METTask_HCAL_D3_MEyvsieta","",83,-41,42, 10001,-500,501);  
00199         hHCAL_D4_MEyvsieta         = dbe_->book2D("METTask_HCAL_D4_MEyvsieta","",83,-41,42, 10001,-500,501);  
00200       }
00201     else
00202       {
00203         hHCAL_D1_energyvsieta = dbe_->book2D("METTask_HCAL_D1_energyvsieta","",83,-41,42,1000,-10,1990);
00204         hHCAL_D2_energyvsieta = dbe_->book2D("METTask_HCAL_D2_energyvsieta","",83,-41,42,1000,-10,1990);
00205         hHCAL_D3_energyvsieta = dbe_->book2D("METTask_HCAL_D3_energyvsieta","",83,-41,42,1000,-10,1990);
00206         hHCAL_D4_energyvsieta = dbe_->book2D("METTask_HCAL_D4_energyvsieta","",83,-41,42,1000,-10,1990);
00207 
00208         hHCAL_D1_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D1_Minenergyvsieta","",83,-41,42,1000,-10,1990);
00209         hHCAL_D2_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D2_Minenergyvsieta","",83,-41,42,1000,-10,1990);
00210         hHCAL_D3_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D3_Minenergyvsieta","",83,-41,42,1000,-10,1990);
00211         hHCAL_D4_Minenergyvsieta = dbe_->book2D("METTask_HCAL_D4_Minenergyvsieta","",83,-41,42,1000,-10,1990);
00212 
00213         hHCAL_D1_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D1_Maxenergyvsieta","",83,-41,42,1000,-10,1990);
00214         hHCAL_D2_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D2_Maxenergyvsieta","",83,-41,42,1000,-10,1990);
00215         hHCAL_D3_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D3_Maxenergyvsieta","",83,-41,42,1000,-10,1990);
00216         hHCAL_D4_Maxenergyvsieta = dbe_->book2D("METTask_HCAL_D4_Maxenergyvsieta","",83,-41,42,1000,-10,1990);
00217 
00218         // Integrated over phi                                                                                                                                                                                 
00219         hHCAL_D1_Occvsieta = dbe_->book2D("METTask_HCAL_D1_Occvsieta","",83,-41,42,73,0,73);
00220         hHCAL_D2_Occvsieta = dbe_->book2D("METTask_HCAL_D2_Occvsieta","",83,-41,42,73,0,73);
00221         hHCAL_D3_Occvsieta = dbe_->book2D("METTask_HCAL_D3_Occvsieta","",83,-41,42,73,0,73);
00222         hHCAL_D4_Occvsieta = dbe_->book2D("METTask_HCAL_D4_Occvsieta","",83,-41,42,73,0,73);
00223 
00224         hHCAL_D1_SETvsieta = dbe_->book2D("METTask_HCAL_D1_SETvsieta","",83,-41,42,1000,0,2000);
00225         hHCAL_D2_SETvsieta = dbe_->book2D("METTask_HCAL_D2_SETvsieta","",83,-41,42,1000,0,2000);
00226         hHCAL_D3_SETvsieta = dbe_->book2D("METTask_HCAL_D3_SETvsieta","",83,-41,42,1000,0,2000);
00227         hHCAL_D4_SETvsieta = dbe_->book2D("METTask_HCAL_D4_SETvsieta","",83,-41,42,1000,0,2000);
00228 
00229         hHCAL_D1_METvsieta = dbe_->book2D("METTask_HCAL_D1_METvsieta","",83,-41,42,1000,0,2000);
00230         hHCAL_D2_METvsieta = dbe_->book2D("METTask_HCAL_D2_METvsieta","",83,-41,42,1000,0,2000);
00231         hHCAL_D3_METvsieta = dbe_->book2D("METTask_HCAL_D3_METvsieta","",83,-41,42,1000,0,2000);
00232         hHCAL_D4_METvsieta = dbe_->book2D("METTask_HCAL_D4_METvsieta","",83,-41,42,1000,0,2000);
00233 
00234         hHCAL_D1_METPhivsieta      = dbe_->book2D("METTask_HCAL_D1_METPhivsieta","",83,-41,42, 80, -4,4);
00235         hHCAL_D2_METPhivsieta      = dbe_->book2D("METTask_HCAL_D2_METPhivsieta","",83,-41,42, 80, -4,4);
00236         hHCAL_D3_METPhivsieta      = dbe_->book2D("METTask_HCAL_D3_METPhivsieta","",83,-41,42, 80, -4,4);
00237         hHCAL_D4_METPhivsieta      = dbe_->book2D("METTask_HCAL_D4_METPhivsieta","",83,-41,42, 80, -4,4);
00238 
00239         hHCAL_D1_MExvsieta         = dbe_->book2D("METTask_HCAL_D1_MExvsieta","",83,-41,42, 500,-500,500);
00240         hHCAL_D2_MExvsieta         = dbe_->book2D("METTask_HCAL_D2_MExvsieta","",83,-41,42, 500,-500,500);
00241         hHCAL_D3_MExvsieta         = dbe_->book2D("METTask_HCAL_D3_MExvsieta","",83,-41,42, 500,-500,500);
00242         hHCAL_D4_MExvsieta         = dbe_->book2D("METTask_HCAL_D4_MExvsieta","",83,-41,42, 500,-500,500);
00243 
00244         hHCAL_D1_MEyvsieta         = dbe_->book2D("METTask_HCAL_D1_MEyvsieta","",83,-41,42, 500,-500,500);
00245         hHCAL_D2_MEyvsieta         = dbe_->book2D("METTask_HCAL_D2_MEyvsieta","",83,-41,42, 500,-500,500);
00246         hHCAL_D3_MEyvsieta         = dbe_->book2D("METTask_HCAL_D3_MEyvsieta","",83,-41,42, 500,-500,500);
00247         hHCAL_D4_MEyvsieta         = dbe_->book2D("METTask_HCAL_D4_MEyvsieta","",83,-41,42, 500,-500,500);
00248         
00249       }
00250     
00251   }
00252   
00253   // Inspect Setup for CaloTower Geometry
00254   //  FillGeometry(iSetup);
00255   
00256 }
00257 
00258 void HCALRecHitAnalyzer::FillGeometry(const edm::EventSetup& iSetup)
00259 {
00260 
00261   // ==========================================================
00262   // Retrieve!
00263   // ==========================================================
00264 
00265   const CaloSubdetectorGeometry* HBgeom;
00266   const CaloSubdetectorGeometry* HEgeom;
00267   const CaloSubdetectorGeometry* HOgeom;
00268   const CaloSubdetectorGeometry* HFgeom;
00269 
00270   edm::ESHandle<CaloGeometry> pG;
00271   iSetup.get<CaloGeometryRecord>().get(pG);
00272   
00273   if (!pG.isValid() )
00274     {
00275       edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Setup Handle, Aborting Task "
00276                                  << "HCALRecHitAnalyzer::FillGeometry!\n"; 
00277       return;
00278     }
00279   
00280   const CaloGeometry cG = *pG;
00281   
00282   HBgeom = cG.getSubdetectorGeometry(DetId::Hcal,HcalBarrel);
00283   HEgeom = cG.getSubdetectorGeometry(DetId::Hcal,HcalEndcap);
00284   HOgeom = cG.getSubdetectorGeometry(DetId::Hcal,HcalOuter);
00285   HFgeom = cG.getSubdetectorGeometry(DetId::Hcal,HcalForward);
00286     
00287   
00288   // ==========================================================
00289   // Fill Histograms!
00290   // ==========================================================
00291 
00292   std::vector<DetId>::iterator i;
00293 
00294   int HBmin_ieta = 99, HBmax_ieta = -99;
00295   int HBmin_iphi = 99, HBmax_iphi = -99;
00296 
00297   // Loop Over all Hcal Barrel DetId's
00298   int nHBdetid = 0;
00299   std::vector<DetId> HBids = HBgeom->getValidDetIds(DetId::Hcal,HcalBarrel);
00300 
00301   for (i = HBids.begin(); i != HBids.end(); i++) {
00302 
00303     nHBdetid++;
00304 
00305     const CaloCellGeometry* cell = HBgeom->getGeometry(*i);
00306     HcalDetId HcalID(i->rawId());
00307     //GlobalPoint p = cell->getPosition();
00308 
00309     int Calo_ieta = 42 + HcalID.ieta();
00310     int Calo_iphi = HcalID.iphi();
00311     double Calo_eta = cell->getPosition().eta();
00312     double Calo_phi = cell->getPosition().phi();
00313 
00314     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta,Calo_iphi) == -999) {
00315 
00316       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta,Calo_iphi,Calo_eta);
00317       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta,Calo_iphi,Calo_phi*180.0/M_PI);
00318 
00319     }
00320 
00321     if (Calo_ieta > HBmax_ieta) HBmax_ieta = Calo_ieta;
00322     if (Calo_ieta < HBmin_ieta) HBmin_ieta = Calo_ieta;
00323     if (Calo_iphi > HBmax_iphi) HBmax_iphi = Calo_iphi;
00324     if (Calo_iphi > HBmax_iphi) HBmin_iphi = Calo_iphi;
00325 
00326   }
00327 
00328   int HEmin_ieta = 99, HEmax_ieta = -99;
00329   int HEmin_iphi = 99, HEmax_iphi = -99;
00330 
00331   // Loop Over all Hcal Endcap DetId's
00332   int nHEdetid = 0;
00333   std::vector<DetId> HEids = HEgeom->getValidDetIds(DetId::Hcal,HcalEndcap);
00334 
00335   for (i = HEids.begin(); i != HEids.end(); i++) {
00336 
00337     nHEdetid++;
00338 
00339     const CaloCellGeometry* cell = HEgeom->getGeometry(*i);
00340     HcalDetId HcalID(i->rawId());
00341     //GlobalPoint p = cell->getPosition();
00342 
00343     int Calo_ieta = 42 + HcalID.ieta();
00344     int Calo_iphi = HcalID.iphi();
00345     double Calo_eta = cell->getPosition().eta();
00346     double Calo_phi = cell->getPosition().phi();
00347 
00348     // HCAL to HE eta, phi map comparison
00349     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta,Calo_iphi) == -999) {
00350       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta,Calo_iphi,Calo_eta);
00351       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta,Calo_iphi,Calo_phi*180.0/M_PI);
00352     } 
00353 
00354     if (Calo_ieta > HEmax_ieta) HEmax_ieta = Calo_ieta;
00355     if (Calo_ieta < HEmin_ieta) HEmin_ieta = Calo_ieta;
00356     if (Calo_iphi > HEmax_iphi) HEmax_iphi = Calo_iphi;
00357     if (Calo_iphi > HEmax_iphi) HEmin_iphi = Calo_iphi;
00358 
00359   }
00360 
00361   int HFmin_ieta = 99, HFmax_ieta = -99;
00362   int HFmin_iphi = 99, HFmax_iphi = -99;
00363 
00364   // Loop Over all Hcal Forward DetId's
00365   int nHFdetid = 0;
00366   std::vector<DetId> HFids = HFgeom->getValidDetIds(DetId::Hcal,HcalForward);
00367 
00368   for (i = HFids.begin(); i != HFids.end(); i++) {
00369 
00370     nHFdetid++;
00371 
00372     const CaloCellGeometry* cell = HFgeom->getGeometry(*i);
00373     HcalDetId HcalID(i->rawId());
00374     //GlobalPoint p = cell->getPosition();
00375 
00376     int Calo_ieta = 42 + HcalID.ieta();
00377     int Calo_iphi = HcalID.iphi();
00378     double Calo_eta = cell->getPosition().eta();
00379     double Calo_phi = cell->getPosition().phi();
00380 
00381     // HCAL to HF eta, phi map comparison
00382     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta,Calo_iphi) == -999) {
00383       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta,Calo_iphi,Calo_eta);
00384       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta,Calo_iphi,Calo_phi*180.0/M_PI);
00385     } 
00386 
00387     if (Calo_ieta > HFmax_ieta) HFmax_ieta = Calo_ieta;
00388     if (Calo_ieta < HFmin_ieta) HFmin_ieta = Calo_ieta;
00389     if (Calo_iphi > HFmax_iphi) HFmax_iphi = Calo_iphi;
00390     if (Calo_iphi > HFmax_iphi) HFmin_iphi = Calo_iphi;
00391 
00392   }
00393 
00394   int HOmin_ieta = 99, HOmax_ieta = -99;
00395   int HOmin_iphi = 99, HOmax_iphi = -99;
00396 
00397   // Loop Over all Hcal Outer DetId's
00398   int nHOdetid = 0;
00399   std::vector<DetId> HOids = HOgeom->getValidDetIds(DetId::Hcal,HcalOuter);
00400 
00401   for (i = HOids.begin(); i != HOids.end(); i++) {
00402 
00403     nHOdetid++;
00404 
00405     const CaloCellGeometry* cell = HOgeom->getGeometry(*i);
00406     HcalDetId HcalID(i->rawId());
00407     //GlobalPoint p = cell->getPosition();
00408 
00409     int Calo_ieta = 42 + HcalID.ieta();
00410     int Calo_iphi = HcalID.iphi();
00411     double Calo_eta = cell->getPosition().eta();
00412     double Calo_phi = cell->getPosition().phi();
00413 
00414     // HCAL to HO eta, phi map comparison
00415     if (hHCAL_ieta_iphi_etaMap->getBinContent(Calo_ieta,Calo_iphi) == -999) {
00416       hHCAL_ieta_iphi_etaMap->setBinContent(Calo_ieta,Calo_iphi,Calo_eta);
00417       hHCAL_ieta_iphi_phiMap->setBinContent(Calo_ieta,Calo_iphi,Calo_phi*180.0/M_PI);
00418     } 
00419 
00420     if (Calo_ieta > HOmax_ieta) HOmax_ieta = Calo_ieta;
00421     if (Calo_ieta < HOmin_ieta) HOmin_ieta = Calo_ieta;
00422     if (Calo_iphi > HOmax_iphi) HOmax_iphi = Calo_iphi;
00423     if (Calo_iphi > HOmax_iphi) HOmin_iphi = Calo_iphi;
00424 
00425   }
00426 
00427   // Set the Cell Size for each (ieta, iphi) Bin
00428   double currentLowEdge_eta = 0; //double currentHighEdge_eta = 0;
00429   for (int ieta = 1; ieta < 42 ; ieta++) {
00430     
00431     int ieta_ = 42 + ieta;
00432     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(ieta_,3);
00433     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(ieta_,3);
00434     double deta = 2.0*(eta-currentLowEdge_eta);
00435     deta = ((float)((int)(1.0E3*deta + 0.5)))/1.0E3;
00436     double dphi = 2.0*phi;
00437     if (ieta==40 || ieta==41) dphi = 20;
00438     if (ieta<=39 && ieta>=21) dphi = 10;
00439     if (ieta<=20) dphi = 5;
00440     // BS: This is WRONG...need to correct overlap 
00441     if (ieta == 28) deta = 0.218;
00442     if (ieta == 29) deta= 0.096;      
00443     currentLowEdge_eta += deta;
00444 
00445     // BS: This is WRONG...need to correct overlap 
00446     if (ieta == 29) currentLowEdge_eta = 2.964;
00447 
00448     hHCAL_ieta_detaMap->setBinContent(ieta_,deta); // positive rings
00449     hHCAL_ieta_dphiMap->setBinContent(ieta_,dphi); // positive rings
00450     hHCAL_ieta_detaMap->setBinContent(42-ieta,deta); // negative rings
00451     hHCAL_ieta_dphiMap->setBinContent(42-ieta,dphi); // negative rings
00452 
00453   } // end loop over ieta
00454 
00455   edm::LogInfo("OutputInfo") << "HB ieta range: " << HBmin_ieta << " " << HBmax_ieta;
00456   edm::LogInfo("OutputInfo") << "HB iphi range: " << HBmin_iphi << " " << HBmax_iphi;
00457   edm::LogInfo("OutputInfo") << "HE ieta range: " << HEmin_ieta << " " << HEmax_ieta;
00458   edm::LogInfo("OutputInfo") << "HE iphi range: " << HEmin_iphi << " " << HEmax_iphi;
00459   edm::LogInfo("OutputInfo") << "HF ieta range: " << HFmin_ieta << " " << HFmax_ieta;
00460   edm::LogInfo("OutputInfo") << "HF iphi range: " << HFmin_iphi << " " << HFmax_iphi;
00461   edm::LogInfo("OutputInfo") << "HO ieta range: " << HOmin_ieta << " " << HOmax_ieta;
00462   edm::LogInfo("OutputInfo") << "HO iphi range: " << HOmin_iphi << " " << HOmax_iphi;
00463 
00464 }
00465 
00466 void HCALRecHitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00467 {
00468     Nevents++;
00469     hHCAL_Nevents->Fill(0);
00470   // ==========================================================
00471   // Retrieve!
00472   // ==========================================================
00473                                                                                                                                                              
00474   const HBHERecHitCollection *HBHERecHits;
00475   const HORecHitCollection *HORecHits;
00476   const HFRecHitCollection *HFRecHits;
00477 
00478   edm::Handle<HBHERecHitCollection> HBHERecHitsHandle;
00479   iEvent.getByLabel(hBHERecHitsLabel_,HBHERecHitsHandle);
00480   if (!HBHERecHitsHandle.isValid()) {
00481     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
00482                                << "HCALRecHitAnalyzer::analyze!\n"; return;
00483   } else {
00484     HBHERecHits = HBHERecHitsHandle.product();
00485   }
00486   edm::Handle<HORecHitCollection> HORecHitsHandle;
00487   iEvent.getByLabel(hORecHitsLabel_,HORecHitsHandle);
00488   if (!HORecHitsHandle.isValid()) {
00489     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
00490                                << "HCALRecHitAnalyzer::analyze!\n"; return;
00491   } else {
00492     HORecHits = HORecHitsHandle.product();
00493   }
00494   edm::Handle<HFRecHitCollection> HFRecHitsHandle;
00495   iEvent.getByLabel(hFRecHitsLabel_,HFRecHitsHandle);
00496   if (!HFRecHitsHandle.isValid()) {
00497     edm::LogInfo("OutputInfo") << "Failed to retrieve an Event Handle, Aborting Task "
00498                                << "HCALRecHitAnalyzer::analyze!\n"; return;
00499   } else {
00500     HFRecHits = HFRecHitsHandle.product();
00501   }
00502 
00503   // ==========================================================
00504   // Fill Histograms!
00505   // ==========================================================
00506                                                    
00507   TLorentzVector vHBHEMET_EtaRing[83][4];
00508   int HBHEActiveRing[83][4];
00509   int HBHENActiveCells[83][4];
00510   double HBHESET_EtaRing[83][4];
00511   double HBHEMinEnergy_EtaRing[83][4];
00512   double HBHEMaxEnergy_EtaRing[83][4];
00513 
00514   for (int i = 0; i < 83; i++) {
00515     for (int j = 0; j < 4; j++) {
00516 
00517       HBHEActiveRing[i][j] = 0;
00518       HBHENActiveCells[i][j] = 0;
00519       HBHESET_EtaRing[i][j] = 0; 
00520       HBHEMinEnergy_EtaRing[i][j] = 14E3; 
00521       HBHEMaxEnergy_EtaRing[i][j] = -999; 
00522     }
00523   }
00524   
00525   // Loop over HBHERecHit's
00526   HBHERecHitCollection::const_iterator hbherechit;
00527   int nHBrechit = 0, nHErechit = 0;
00528   
00529   for (hbherechit = HBHERecHits->begin(); hbherechit != HBHERecHits->end(); hbherechit++) {
00530     
00531     HcalDetId det = hbherechit->id();
00532     double Energy = hbherechit->energy();
00533     Int_t depth = det.depth();
00534     Int_t ieta = det.ieta();
00535     Int_t iphi = det.iphi();
00536     int EtaRing = 41 + ieta; // this counts from 0    
00537     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing+1,iphi);
00538     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing+1,iphi);
00539     double theta = 2*TMath::ATan(exp(-1*eta));
00540     double ET = Energy*TMath::Sin(theta);
00541     HcalSubdetector HcalNum = det.subdet();
00542     TLorentzVector v_;
00543     
00544     
00545     
00546     if (Energy>0) // zero suppress
00547       {
00548         HBHEActiveRing[EtaRing][depth-1] = 1;
00549         HBHENActiveCells[EtaRing][depth-1]++;
00550         HBHESET_EtaRing[EtaRing][depth-1]+=ET;
00551         v_.SetPtEtaPhiE(ET, 0, phi, ET);
00552         vHBHEMET_EtaRing[EtaRing][depth-1]-=v_;
00553         
00554         
00555         DEBUG( EtaRing << " " << Energy << " " << ET << " " << theta << " " << eta << " " << phi << " : " << vHBHEMET_EtaRing[EtaRing][depth-1].Phi() << " " << vHBHEMET_EtaRing[EtaRing][depth-1].Pt() );  
00556         
00557         switch (depth) {
00558         case 1:
00559           hHCAL_D1_Occ_ieta_iphi->Fill(ieta,iphi);
00560           break;
00561         case 2:
00562           hHCAL_D2_Occ_ieta_iphi->Fill(ieta,iphi);
00563           break;
00564         case 3:
00565           hHCAL_D3_Occ_ieta_iphi->Fill(ieta,iphi);
00566           break;
00567         } // end switch
00568       }
00569     
00570     if (Energy > HBHEMaxEnergy_EtaRing[EtaRing][depth-1])
00571       HBHEMaxEnergy_EtaRing[EtaRing][depth-1] = Energy;
00572     if (Energy < HBHEMinEnergy_EtaRing[EtaRing][depth-1])
00573       HBHEMinEnergy_EtaRing[EtaRing][depth-1] = Energy;
00574     
00575     switch (depth) {
00576     case 1:
00577       hHCAL_D1_energy_ieta_iphi->Fill(ieta,iphi,Energy);           
00578       hHCAL_D1_energyvsieta->Fill(ieta,Energy);    
00579       if (Energy>hHCAL_D1_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00580         hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00581       if (Energy<hHCAL_D1_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00582         hHCAL_D1_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00583       
00584       break;
00585     case 2:
00586       hHCAL_D2_energy_ieta_iphi->Fill(ieta,iphi,Energy);
00587       
00588       hHCAL_D2_energyvsieta->Fill(ieta,Energy);
00589       if (Energy>hHCAL_D2_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00590         hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00591       if (Energy<hHCAL_D2_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00592         hHCAL_D2_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00593       break;
00594     case 3:
00595       hHCAL_D3_energy_ieta_iphi->Fill(ieta,iphi,Energy);
00596       
00597       hHCAL_D3_energyvsieta->Fill(ieta,Energy);
00598       if (Energy>hHCAL_D3_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00599         hHCAL_D3_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00600       if (Energy<hHCAL_D3_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00601         hHCAL_D3_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00602       break;
00603     } // end switch
00604     
00605     
00606     if (HcalNum == HcalBarrel) {
00607       nHBrechit++;
00608     } else { // HcalEndcap
00609       nHErechit++;
00610     }
00611   } // end loop over HBHERecHit's
00612   
00613   // Fill eta-ring MET quantities
00614   for (int iEtaRing=0; iEtaRing < 83; iEtaRing++) {
00615     for (int jDepth=0; jDepth < 3; jDepth++) {
00616       
00617       switch (jDepth+1) {
00618       case 1:
00619         hHCAL_D1_Maxenergyvsieta->Fill(iEtaRing-41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
00620         hHCAL_D1_Minenergyvsieta->Fill(iEtaRing-41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
00621         break;
00622       case 2:
00623         hHCAL_D2_Maxenergyvsieta->Fill(iEtaRing-41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
00624         hHCAL_D2_Minenergyvsieta->Fill(iEtaRing-41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
00625         break;
00626       case 3:
00627         hHCAL_D3_Maxenergyvsieta->Fill(iEtaRing-41, HBHEMaxEnergy_EtaRing[iEtaRing][jDepth]);
00628         hHCAL_D3_Minenergyvsieta->Fill(iEtaRing-41, HBHEMinEnergy_EtaRing[iEtaRing][jDepth]);
00629         break;
00630       }
00631       
00632       if (HBHEActiveRing[iEtaRing][jDepth]) {
00633         
00634         switch (jDepth+1) {
00635         case 1:
00636           hHCAL_D1_METPhivsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
00637           hHCAL_D1_MExvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
00638           hHCAL_D1_MEyvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
00639           hHCAL_D1_METvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
00640           hHCAL_D1_SETvsieta->Fill(iEtaRing-41, HBHESET_EtaRing[iEtaRing][jDepth]);
00641           hHCAL_D1_Occvsieta->Fill(iEtaRing-41, HBHENActiveCells[iEtaRing][jDepth]);
00642           break;
00643         case 2:
00644           hHCAL_D2_METPhivsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
00645           hHCAL_D2_MExvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
00646           hHCAL_D2_MEyvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
00647           hHCAL_D2_METvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
00648           hHCAL_D2_SETvsieta->Fill(iEtaRing-41, HBHESET_EtaRing[iEtaRing][jDepth]);
00649           hHCAL_D2_Occvsieta->Fill(iEtaRing-41, HBHENActiveCells[iEtaRing][jDepth]);
00650           break;
00651         case 3:
00652           hHCAL_D3_METPhivsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Phi());
00653           hHCAL_D3_MExvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Px());
00654           hHCAL_D3_MEyvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Py());
00655           hHCAL_D3_METvsieta->Fill(iEtaRing-41, vHBHEMET_EtaRing[iEtaRing][jDepth].Pt());
00656           hHCAL_D3_SETvsieta->Fill(iEtaRing-41, HBHESET_EtaRing[iEtaRing][jDepth]);
00657           hHCAL_D3_Occvsieta->Fill(iEtaRing-41, HBHENActiveCells[iEtaRing][jDepth]);
00658           break;
00659         } //switch
00660       } //activering
00661     }//depth
00662   } //etaring    
00663 
00664   
00665   //-------------------------------------------------HO
00666   TLorentzVector vHOMET_EtaRing[83];
00667   int HOActiveRing[83];
00668   int HONActiveCells[83];
00669   double HOSET_EtaRing[83];
00670   double HOMinEnergy_EtaRing[83];
00671   double HOMaxEnergy_EtaRing[83];
00672   
00673   for (int i=0;i<83; i++) {
00674 
00675     HOActiveRing[i] = 0;
00676     HONActiveCells[i] = 0;
00677     HOSET_EtaRing[i] = 0; 
00678     HOMinEnergy_EtaRing[i] = 14E3;
00679     HOMaxEnergy_EtaRing[i] = -999;
00680   }
00681 
00682   // Loop over HORecHit's
00683   HORecHitCollection::const_iterator horechit;
00684   int nHOrechit = 0;
00685 
00686   for (horechit = HORecHits->begin(); horechit != HORecHits->end(); horechit++) {
00687                                                                                                                                                              
00688     nHOrechit++;
00689                                                                                                                                                              
00690     HcalDetId det = horechit->id();
00691     double Energy = horechit->energy();
00693     Int_t ieta = det.ieta();
00694     Int_t iphi = det.iphi();
00695     int EtaRing = 41+ieta; // this counts from 0    
00696     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing+1,iphi);
00697     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing+1,iphi);
00698     double theta = 2*TMath::ATan(exp(-1*eta));
00699     double ET = Energy*TMath::Sin(theta);
00700     TLorentzVector v_; 
00701 
00702 
00703     if (Energy>0) // zero suppress
00704       {
00705         HOActiveRing[EtaRing] = 1;
00706         HONActiveCells[EtaRing]++;
00707         HOSET_EtaRing[EtaRing]+=ET;
00708         v_.SetPtEtaPhiE(ET, 0, phi, ET);
00709         vHOMET_EtaRing[EtaRing]-=v_;
00710       
00711         hHCAL_D4_Occ_ieta_iphi->Fill(ieta,iphi);
00712       }
00713 
00714     if (Energy > HOMaxEnergy_EtaRing[EtaRing])
00715       HOMaxEnergy_EtaRing[EtaRing] = Energy;
00716     if (Energy < HOMinEnergy_EtaRing[EtaRing])
00717       HOMinEnergy_EtaRing[EtaRing] = Energy;
00718     
00719     hHCAL_D4_energy_ieta_iphi->Fill(ieta,iphi,Energy);
00720     
00721     hHCAL_D4_energyvsieta->Fill(ieta,Energy);     
00722     if (Energy>hHCAL_D4_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00723         hHCAL_D4_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00724       if (Energy<hHCAL_D4_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00725         hHCAL_D4_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00726                                                                                                                                                    
00727   } // end loop over HORecHit's
00728 
00729   // Fill eta-ring MET quantities
00730   for (int iEtaRing=0; iEtaRing<83; iEtaRing++) {
00731     hHCAL_D4_Maxenergyvsieta->Fill(iEtaRing-41, HOMaxEnergy_EtaRing[iEtaRing]);
00732     hHCAL_D4_Minenergyvsieta->Fill(iEtaRing-41, HOMinEnergy_EtaRing[iEtaRing]);
00733 
00734     if (HOActiveRing[iEtaRing]) {
00735 
00736       hHCAL_D4_METPhivsieta->Fill(iEtaRing-41, vHOMET_EtaRing[iEtaRing].Phi());
00737       hHCAL_D4_MExvsieta->Fill(iEtaRing-41, vHOMET_EtaRing[iEtaRing].Px());
00738       hHCAL_D4_MEyvsieta->Fill(iEtaRing-41, vHOMET_EtaRing[iEtaRing].Py());
00739       hHCAL_D4_METvsieta->Fill(iEtaRing-41, vHOMET_EtaRing[iEtaRing].Pt());
00740       hHCAL_D4_SETvsieta->Fill(iEtaRing-41, HOSET_EtaRing[iEtaRing]);
00741       hHCAL_D4_Occvsieta->Fill(iEtaRing-41, HONActiveCells[iEtaRing]);
00742 
00743     }
00744 
00745   }
00746 
00747   //------------------------------------------------HF
00748 
00749   TLorentzVector vHFMET_EtaRing[83][2];
00750   int HFActiveRing[83][2];
00751   int HFNActiveCells[83][2];
00752   double HFSET_EtaRing[83][2];
00753   double HFMinEnergy_EtaRing[83][2];
00754   double HFMaxEnergy_EtaRing[83][2];
00755 
00756   for (int i=0;i<83; i++) {
00757     for (int j=0;j<2; j++) {
00758 
00759       HFActiveRing[i][j] = 0;
00760       HFNActiveCells[i][j] = 0;
00761       HFSET_EtaRing[i][j] = 0; 
00762       HFMinEnergy_EtaRing[i][j] = 14E3; 
00763       HFMaxEnergy_EtaRing[i][j] = -999; 
00764       
00765     }
00766   }
00767 
00768   // Loop over HFRecHit's
00769   HFRecHitCollection::const_iterator hfrechit;
00770   int nHFrechit = 0;
00771 
00772   for (hfrechit = HFRecHits->begin(); hfrechit != HFRecHits->end(); hfrechit++) {
00773                                                                                                                                                              
00774     nHFrechit++;
00775                                                                                                                                                              
00776     HcalDetId det = hfrechit->id();
00777     double Energy = hfrechit->energy();
00778     Int_t depth = det.depth();
00779     Int_t ieta = det.ieta();
00780     Int_t iphi = det.iphi();
00781     int EtaRing = 41+ieta; // this counts from 0    
00782     double eta = hHCAL_ieta_iphi_etaMap->getBinContent(EtaRing+1,iphi);
00783     double phi = hHCAL_ieta_iphi_phiMap->getBinContent(EtaRing+1,iphi);
00784     double theta = 2*TMath::ATan(exp(-1*eta));
00785     double ET = Energy*TMath::Sin(theta);
00786 
00787     TLorentzVector v_;
00788     if (Energy>0) // zero suppress
00789       {
00790         HFActiveRing[EtaRing][depth-1] = 1;
00791         HFNActiveCells[EtaRing][depth-1]++;
00792         HFSET_EtaRing[EtaRing][depth-1]+=ET;
00793         v_.SetPtEtaPhiE(ET, 0, phi, ET);
00794         vHFMET_EtaRing[EtaRing][depth-1]-=v_;                                                                         
00795         
00796         switch (depth) {
00797         case 1:
00798           hHCAL_D1_Occ_ieta_iphi->Fill(ieta,iphi);
00799           break;
00800         case 2:
00801           hHCAL_D2_Occ_ieta_iphi->Fill(ieta,iphi);
00802           break;
00803           
00804         }
00805       }
00806     
00807     if (Energy > HFMaxEnergy_EtaRing[EtaRing][depth-1])
00808       HFMaxEnergy_EtaRing[EtaRing][depth-1] = Energy;
00809     if (Energy < HFMinEnergy_EtaRing[EtaRing][depth-1])
00810       HFMinEnergy_EtaRing[EtaRing][depth-1] = Energy;
00811     
00812 
00813     switch (depth) {
00814     case 1:
00815       hHCAL_D1_energy_ieta_iphi->Fill(ieta,iphi,Energy);
00816       
00817       hHCAL_D1_energyvsieta->Fill(ieta,Energy);  
00818       if (Energy>hHCAL_D1_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00819         hHCAL_D1_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00820       if (Energy<hHCAL_D1_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00821         hHCAL_D1_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy); 
00822       break;
00823     case 2:
00824       hHCAL_D2_energy_ieta_iphi->Fill(ieta,iphi,Energy);
00825       
00826       hHCAL_D2_energyvsieta->Fill(ieta,Energy);  
00827       if (Energy>hHCAL_D2_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00828         hHCAL_D2_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00829       if (Energy<hHCAL_D2_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi)) 
00830         hHCAL_D2_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00831       break;
00832     }
00833 
00834   } // end loop over HFRecHit's
00835 
00836   // Fill eta-ring MET quantities
00837   for (int iEtaRing=0; iEtaRing<83; iEtaRing++) {
00838     for (int jDepth=0; jDepth<2; jDepth++) {
00839       switch (jDepth+1) {
00840       case 1:
00841         hHCAL_D1_Maxenergyvsieta->Fill(iEtaRing-41, HFMaxEnergy_EtaRing[iEtaRing][jDepth]);
00842         hHCAL_D1_Minenergyvsieta->Fill(iEtaRing-41, HFMinEnergy_EtaRing[iEtaRing][jDepth]);
00843         break;
00844       case 2:
00845         hHCAL_D2_Maxenergyvsieta->Fill(iEtaRing-41, HFMaxEnergy_EtaRing[iEtaRing][jDepth]);
00846         hHCAL_D2_Minenergyvsieta->Fill(iEtaRing-41, HFMinEnergy_EtaRing[iEtaRing][jDepth]);
00847         break;
00848       }
00849 
00850 
00851       if (HFActiveRing[iEtaRing][jDepth]) {
00852 
00853         switch (jDepth+1) {
00854         case 1:
00855 
00856           hHCAL_D1_METPhivsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Phi());
00857           hHCAL_D1_MExvsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Px());
00858           hHCAL_D1_MEyvsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Py());
00859           hHCAL_D1_METvsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Pt());
00860           hHCAL_D1_SETvsieta->Fill(iEtaRing-41, HFSET_EtaRing[iEtaRing][jDepth]);
00861           hHCAL_D1_Occvsieta->Fill(iEtaRing-41, HFNActiveCells[iEtaRing][jDepth]);
00862           break;
00863 
00864         case 2:
00865 
00866           hHCAL_D2_METPhivsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Phi());
00867           hHCAL_D2_MExvsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Px());
00868           hHCAL_D2_MEyvsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Py());
00869           hHCAL_D2_METvsieta->Fill(iEtaRing-41, vHFMET_EtaRing[iEtaRing][jDepth].Pt());
00870           hHCAL_D2_SETvsieta->Fill(iEtaRing-41, HFSET_EtaRing[iEtaRing][jDepth]);
00871           hHCAL_D2_Occvsieta->Fill(iEtaRing-41, HFNActiveCells[iEtaRing][jDepth]);
00872           break;
00873 
00874         }
00875       }
00876     }
00877   }
00878 
00879 
00880 
00881 }
00882 
00883 
00884 void HCALRecHitAnalyzer::endJob()
00885 {
00886 
00887 
00888 }
/data/refman/pasoursint/CMSSW_5_3_10_patch2/src/DQMOffline/JetMET/src/HCALRecHitAnalyzer.cc
