Reference Manual

Go to the documentation of this file.
00001 #include "DQMOffline/JetMET/interface/ECALRecHitAnalyzer.h"
00002 #include "DQMServices/Core/interface/DQMStore.h"
00003 #include "Geometry/Records/interface/CaloGeometryRecord.h"
00004 
00005 // author: Bobby Scurlock, University of Florida
00006 // first version 11/20/2006
00007 
00008 #define DEBUG(X) { if (debug_) { std::cout << X << std::endl; } }
00009 
00010 ECALRecHitAnalyzer::ECALRecHitAnalyzer(const edm::ParameterSet& iConfig)
00011 {
00012   
00013   // Retrieve Information from the Configuration File
00014   EBRecHitsLabel_  = iConfig.getParameter<edm::InputTag>("EBRecHitsLabel");
00015   EERecHitsLabel_  = iConfig.getParameter<edm::InputTag>("EERecHitsLabel");
00016   FolderName_      = iConfig.getUntrackedParameter<std::string>("FolderName");
00017   debug_           = iConfig.getParameter<bool>("Debug");
00018 
00019 
00020 }
00021 
00022 void ECALRecHitAnalyzer::endJob() {
00023 
00024 } 
00025 
00026 //void ECALRecHitAnalyzer::beginJob(void){
00027 void ECALRecHitAnalyzer::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup){
00028   CurrentEvent = -1;
00029   // Book the Histograms
00030   // Fill the geometry histograms
00031   BookHistos();
00032   FillGeometry(iSetup);
00033 }
00034 
00035 void ECALRecHitAnalyzer::BookHistos()
00036 {
00037   // get ahold of back-end interface
00038   dbe_ = edm::Service<DQMStore>().operator->();
00039   
00040   if (dbe_) {
00041     
00042     // Book Geometry Histograms
00043     dbe_->setCurrentFolder(FolderName_+"/geometry");
00044     
00045 
00046     // ECAL barrel
00047     hEB_ieta_iphi_etaMap = dbe_->book2D("hEB_ieta_iphi_etaMap","", 171, -85, 86, 360, 1, 361);
00048     hEB_ieta_iphi_phiMap = dbe_->book2D("hEB_ieta_iphi_phiMap","", 171, -85, 86, 360, 1, 361);
00049     hEB_ieta_detaMap = dbe_->book1D("hEB_ieta_detaMap","", 171, -85, 86);
00050     hEB_ieta_dphiMap = dbe_->book1D("hEB_ieta_dphiMap","", 171, -85, 86);
00051     // ECAL +endcap
00052     hEEpZ_ix_iy_irMap = dbe_->book2D("hEEpZ_ix_iy_irMap","", 100,1,101, 100,1,101);
00053     hEEpZ_ix_iy_xMap = dbe_->book2D("hEEpZ_ix_iy_xMap","", 100,1,101, 100,1,101);
00054     hEEpZ_ix_iy_yMap = dbe_->book2D("hEEpZ_ix_iy_yMap","", 100,1,101, 100,1,101);
00055     hEEpZ_ix_iy_zMap = dbe_->book2D("hEEpZ_ix_iy_zMap","", 100,1,101, 100,1,101);
00056     hEEpZ_ix_iy_dxMap = dbe_->book2D("hEEpZ_ix_iy_dxMap","", 100,1,101, 100,1,101);  
00057     hEEpZ_ix_iy_dyMap = dbe_->book2D("hEEpZ_ix_iy_dyMap","", 100,1,101, 100,1,101);
00058     // ECAL -endcap
00059     hEEmZ_ix_iy_irMap = dbe_->book2D("hEEmZ_ix_iy_irMap","", 100,1,101, 100,1,101);
00060     hEEmZ_ix_iy_xMap = dbe_->book2D("hEEmZ_ix_iy_xMap","", 100,1,101, 100,1,101);
00061     hEEmZ_ix_iy_yMap = dbe_->book2D("hEEmZ_ix_iy_yMap","", 100,1,101, 100,1,101);
00062     hEEmZ_ix_iy_zMap = dbe_->book2D("hEEmZ_ix_iy_zMap","", 100,1,101, 100,1,101);
00063     hEEmZ_ix_iy_dxMap = dbe_->book2D("hEEmZ_ix_iy_dxMap","", 100,1,101, 100,1,101);  
00064     hEEmZ_ix_iy_dyMap = dbe_->book2D("hEEmZ_ix_iy_dyMap","", 100,1,101, 100,1,101);
00065 
00066     // Initialize bins for geometry to -999 because z = 0 is a valid entry 
00067     for (int i=1; i<=100; i++)
00068       for (int j=1; j<=100; j++)
00069         {
00070           hEEpZ_ix_iy_irMap->setBinContent(i,j,-999);
00071           hEEpZ_ix_iy_xMap->setBinContent(i,j,-999);
00072           hEEpZ_ix_iy_yMap->setBinContent(i,j,-999);
00073           hEEpZ_ix_iy_zMap->setBinContent(i,j,-999);
00074           hEEpZ_ix_iy_dxMap->setBinContent(i,j,-999);
00075           hEEpZ_ix_iy_dyMap->setBinContent(i,j,-999);
00076 
00077           hEEmZ_ix_iy_irMap->setBinContent(i,j,-999);
00078           hEEmZ_ix_iy_xMap->setBinContent(i,j,-999);
00079           hEEmZ_ix_iy_yMap->setBinContent(i,j,-999);
00080           hEEmZ_ix_iy_zMap->setBinContent(i,j,-999);
00081           hEEmZ_ix_iy_dxMap->setBinContent(i,j,-999);
00082           hEEmZ_ix_iy_dyMap->setBinContent(i,j,-999);
00083         }
00084 
00085     for (int i=1; i<=171; i++)
00086       {
00087         hEB_ieta_detaMap->setBinContent(i,-999);
00088         hEB_ieta_dphiMap->setBinContent(i,-999);
00089         for (int j=1; j<=360; j++)
00090           {
00091             hEB_ieta_iphi_etaMap->setBinContent(i,j,-999);
00092             hEB_ieta_iphi_phiMap->setBinContent(i,j,-999);
00093           }
00094       }
00095 
00096     // Book Data Histograms
00097     dbe_->setCurrentFolder(FolderName_);
00098 
00099     hECAL_Nevents          = dbe_->book1D("hECAL_Nevents","",1,0,1); 
00100   
00101 
00102     // Energy Histograms by logical index
00103     hEEpZ_energy_ix_iy = dbe_->book2D("hEEpZ_energy_ix_iy","", 100,1,101, 100,1,101);
00104     hEEmZ_energy_ix_iy = dbe_->book2D("hEEmZ_energy_ix_iy","", 100,1,101, 100,1,101);
00105     hEB_energy_ieta_iphi = dbe_->book2D("hEB_energy_ieta_iphi","", 171, -85, 86, 360, 1, 361);   
00106 
00107     hEEpZ_Minenergy_ix_iy = dbe_->book2D("hEEpZ_Minenergy_ix_iy","", 100,1,101, 100,1,101);
00108     hEEmZ_Minenergy_ix_iy = dbe_->book2D("hEEmZ_Minenergy_ix_iy","", 100,1,101, 100,1,101);
00109     hEB_Minenergy_ieta_iphi = dbe_->book2D("hEB_Minenergy_ieta_iphi","", 171, -85, 86, 360, 1, 361);   
00110 
00111     hEEpZ_Maxenergy_ix_iy = dbe_->book2D("hEEpZ_Maxenergy_ix_iy","", 100,1,101, 100,1,101);
00112     hEEmZ_Maxenergy_ix_iy = dbe_->book2D("hEEmZ_Maxenergy_ix_iy","", 100,1,101, 100,1,101);
00113     hEB_Maxenergy_ieta_iphi = dbe_->book2D("hEB_Maxenergy_ieta_iphi","", 171, -85, 86, 360, 1, 361);   
00114 
00115     // need to initialize those
00116     for (int i=1; i<=171; i++)
00117       for (int j=1; j<=360; j++)
00118         {
00119           hEB_Maxenergy_ieta_iphi->setBinContent(i,j,-999);
00120           hEB_Minenergy_ieta_iphi->setBinContent(i,j,14000);
00121         }
00122     for (int i=1; i<=100; i++)
00123       for (int j=1; j<=100; j++)
00124         {
00125           hEEpZ_Maxenergy_ix_iy->setBinContent(i,j,-999);
00126           hEEpZ_Minenergy_ix_iy->setBinContent(i,j,14000);
00127           hEEmZ_Maxenergy_ix_iy->setBinContent(i,j,-999);
00128           hEEmZ_Minenergy_ix_iy->setBinContent(i,j,14000);
00129         }
00130   
00131 
00132     // Occupancy Histograms by logical index
00133     hEEpZ_Occ_ix_iy = dbe_->book2D("hEEpZ_Occ_ix_iy","", 100,1,101, 100,1,101);  
00134     hEEmZ_Occ_ix_iy = dbe_->book2D("hEEmZ_Occ_ix_iy","", 100,1,101, 100,1,101);  
00135     hEB_Occ_ieta_iphi = dbe_->book2D("hEB_Occ_ieta_iphi","",171, -85, 86, 360, 1, 361);   
00136 
00137     // Integrated Histograms
00138     if(finebinning_)
00139       {
00140         hEEpZ_energyvsir = dbe_->book2D("hEEpZ_energyvsir","", 100,1,101, 20110,-10,201);
00141         hEEmZ_energyvsir = dbe_->book2D("hEEmZ_energyvsir","", 100,1,101, 20110,-10,201);
00142         hEB_energyvsieta = dbe_->book2D("hEB_energyvsieta","", 171, -85, 86, 20110, -10, 201);   
00143       
00144         hEEpZ_Maxenergyvsir = dbe_->book2D("hEEpZ_Maxenergyvsir","", 100,1,101, 20110,-10,201);
00145         hEEmZ_Maxenergyvsir = dbe_->book2D("hEEmZ_Maxenergyvsir","", 100,1,101, 20110,-10,201);
00146         hEB_Maxenergyvsieta = dbe_->book2D("hEB_Maxenergyvsieta","", 171, -85, 86, 20110, -10, 201);   
00147       
00148         hEEpZ_Minenergyvsir = dbe_->book2D("hEEpZ_Minenergyvsir","", 100,1,101, 20110,-10,201);
00149         hEEmZ_Minenergyvsir = dbe_->book2D("hEEmZ_Minenergyvsir","", 100,1,101, 20110,-10,201);
00150         hEB_Minenergyvsieta = dbe_->book2D("hEB_Minenergyvsieta","", 171, -85, 86, 20110, -10, 201);   
00151       
00152         hEEpZ_SETvsir = dbe_->book2D("hEEpZ_SETvsir","", 50,1,51, 20010,0,201);
00153         hEEmZ_SETvsir = dbe_->book2D("hEEmZ_SETvsir","", 50,1,51, 20010,0,201);
00154         hEB_SETvsieta = dbe_->book2D("hEB_SETvsieta","", 171, -85, 86, 20010, 0, 201);   
00155       
00156         hEEpZ_METvsir = dbe_->book2D("hEEpZ_METvsir","", 50,1,51, 20010,0,201);
00157         hEEmZ_METvsir = dbe_->book2D("hEEmZ_METvsir","", 50,1,51, 20010,0,201);
00158         hEB_METvsieta = dbe_->book2D("hEB_METvsieta","", 171, -85, 86, 20010, 0, 201);   
00159       
00160         hEEpZ_METPhivsir = dbe_->book2D("hEEpZ_METPhivsir","", 50,1,51, 80,-4,4);
00161         hEEmZ_METPhivsir = dbe_->book2D("hEEmZ_METPhivsir","", 50,1,51, 80,-4,4);
00162         hEB_METPhivsieta = dbe_->book2D("hEB_METPhivsieta","", 171, -85, 86, 80,-4,4);   
00163       
00164         hEEpZ_MExvsir = dbe_->book2D("hEEpZ_MExvsir","", 50,1,51, 10010,-50,51);
00165         hEEmZ_MExvsir = dbe_->book2D("hEEmZ_MExvsir","", 50,1,51, 10010,-50,51);
00166         hEB_MExvsieta = dbe_->book2D("hEB_MExvsieta","", 171, -85, 86, 10010,-50,51);   
00167       
00168         hEEpZ_MEyvsir = dbe_->book2D("hEEpZ_MEyvsir","", 50,1,51, 10010,-50,51);
00169         hEEmZ_MEyvsir = dbe_->book2D("hEEmZ_MEyvsir","", 50,1,51, 10010,-50,51);
00170         hEB_MEyvsieta = dbe_->book2D("hEB_MEyvsieta","", 171, -85, 86, 10010,-50,51);   
00171       
00172         hEEpZ_Occvsir = dbe_->book2D("hEEpZ_Occvsir","", 50,1,51, 1000,0,1000);
00173         hEEmZ_Occvsir = dbe_->book2D("hEEmZ_Occvsir","", 50,1,51, 1000,0,1000);
00174         hEB_Occvsieta = dbe_->book2D("hEB_Occvsieta","", 171, -85, 86, 400,0,400);   
00175       }
00176     else 
00177       {
00178         hEEpZ_energyvsir = dbe_->book2D("hEEpZ_energyvsir","", 100,1,101, 510,-10,100);
00179         hEEmZ_energyvsir = dbe_->book2D("hEEmZ_energyvsir","", 100,1,101, 510,-10,100);
00180         hEB_energyvsieta = dbe_->book2D("hEB_energyvsieta","", 171, -85, 86, 510, -10, 100);
00181       
00182         hEEpZ_Maxenergyvsir = dbe_->book2D("hEEpZ_Maxenergyvsir","", 100,1,101, 510,-10,100);
00183         hEEmZ_Maxenergyvsir = dbe_->book2D("hEEmZ_Maxenergyvsir","", 100,1,101, 510,-10,100);
00184         hEB_Maxenergyvsieta = dbe_->book2D("hEB_Maxenergyvsieta","", 171, -85, 86, 510, -10, 100);
00185 
00186         hEEpZ_Minenergyvsir = dbe_->book2D("hEEpZ_Minenergyvsir","", 100,1,101, 510,-10,100);
00187         hEEmZ_Minenergyvsir = dbe_->book2D("hEEmZ_Minenergyvsir","", 100,1,101, 510,-10,100);
00188         hEB_Minenergyvsieta = dbe_->book2D("hEB_Minenergyvsieta","", 171, -85, 86, 510, -10, 100);
00189 
00190         hEEpZ_SETvsir = dbe_->book2D("hEEpZ_SETvsir","", 50,1,51, 510,0,100);
00191         hEEmZ_SETvsir = dbe_->book2D("hEEmZ_SETvsir","", 50,1,51, 510,0,100);
00192         hEB_SETvsieta = dbe_->book2D("hEB_SETvsieta","", 171, -85, 86, 510, 0, 100);
00193 
00194         hEEpZ_METvsir = dbe_->book2D("hEEpZ_METvsir","", 50,1,51, 510,0,100);
00195         hEEmZ_METvsir = dbe_->book2D("hEEmZ_METvsir","", 50,1,51, 510,0,100);
00196         hEB_METvsieta = dbe_->book2D("hEB_METvsieta","", 171, -85, 86, 510, 0, 100);
00197 
00198         hEEpZ_METPhivsir = dbe_->book2D("hEEpZ_METPhivsir","", 50,1,51, 80,-4,4);
00199         hEEmZ_METPhivsir = dbe_->book2D("hEEmZ_METPhivsir","", 50,1,51, 80,-4,4);
00200         hEB_METPhivsieta = dbe_->book2D("hEB_METPhivsieta","", 171, -85, 86, 80,-4,4);
00201 
00202         hEEpZ_MExvsir = dbe_->book2D("hEEpZ_MExvsir","", 50,1,51, 510,-50,51);
00203         hEEmZ_MExvsir = dbe_->book2D("hEEmZ_MExvsir","", 50,1,51, 510,-50,51);
00204         hEB_MExvsieta = dbe_->book2D("hEB_MExvsieta","", 171, -85, 86, 510,-50,51);
00205 
00206         hEEpZ_MEyvsir = dbe_->book2D("hEEpZ_MEyvsir","", 50,1,51, 510,-50,51);
00207         hEEmZ_MEyvsir = dbe_->book2D("hEEmZ_MEyvsir","", 50,1,51, 510,-50,51);
00208         hEB_MEyvsieta = dbe_->book2D("hEB_MEyvsieta","", 171, -85, 86, 510,-50,51);
00209 
00210         hEEpZ_Occvsir = dbe_->book2D("hEEpZ_Occvsir","", 50,1,51, 1000,0,1000);
00211         hEEmZ_Occvsir = dbe_->book2D("hEEmZ_Occvsir","", 50,1,51, 1000,0,1000);
00212         hEB_Occvsieta = dbe_->book2D("hEB_Occvsieta","", 171, -85, 86, 400,0,400);
00213 
00214       }
00215 
00216 
00217 
00218   }
00219 }
00220 
00221 void ECALRecHitAnalyzer::FillGeometry(const edm::EventSetup& iSetup)
00222 {
00223   // Fill geometry histograms
00224   using namespace edm;
00225   //int b=0;
00226   edm::ESHandle<CaloGeometry> pG;
00227   iSetup.get<CaloGeometryRecord>().get(pG);
00228   const CaloGeometry cG = *pG;
00229   
00230   //----Fill Ecal Barrel----//
00231   const CaloSubdetectorGeometry* EBgeom=cG.getSubdetectorGeometry(DetId::Ecal,1);
00232   int n=0;
00233   std::vector<DetId> EBids=EBgeom->getValidDetIds(DetId::Ecal, 1);
00234   for (std::vector<DetId>::iterator i=EBids.begin(); i!=EBids.end(); i++) {
00235     n++;
00236     const CaloCellGeometry* cell=EBgeom->getGeometry(*i);
00237     //GlobalPoint p = cell->getPosition();
00238     
00239     EBDetId EcalID(i->rawId());
00240     
00241     int Crystal_ieta = EcalID.ieta();
00242     int Crystal_iphi = EcalID.iphi();
00243     double Crystal_eta = cell->getPosition().eta();
00244     double Crystal_phi = cell->getPosition().phi();
00245     hEB_ieta_iphi_etaMap->setBinContent(Crystal_ieta+86, Crystal_iphi, Crystal_eta);
00246     hEB_ieta_iphi_phiMap->setBinContent(Crystal_ieta+86, Crystal_iphi, (Crystal_phi*180/M_PI) );
00247     
00248     DEBUG( " Crystal " << n );
00249     DEBUG( "  ieta, iphi = " << Crystal_ieta << ", " << Crystal_iphi);
00250     DEBUG( "   eta,  phi = " << cell->getPosition().eta() << ", " << cell->getPosition().phi());
00251     DEBUG( " " );
00252     
00253   }
00254   //----Fill Ecal Endcap----------//
00255   const CaloSubdetectorGeometry* EEgeom=cG.getSubdetectorGeometry(DetId::Ecal,2);
00256   n=0;
00257   std::vector<DetId> EEids=EEgeom->getValidDetIds(DetId::Ecal, 2);
00258   for (std::vector<DetId>::iterator i=EEids.begin(); i!=EEids.end(); i++) {
00259     n++;
00260     const CaloCellGeometry* cell=EEgeom->getGeometry(*i);
00261     //GlobalPoint p = cell->getPosition();
00262     EEDetId EcalID(i->rawId());
00263     int Crystal_zside = EcalID.zside();
00264     int Crystal_ix = EcalID.ix();
00265     int Crystal_iy = EcalID.iy();
00266     Float_t ix_ = Crystal_ix-50.5;
00267     Float_t iy_ = Crystal_iy-50.5;
00268     Int_t ir = (Int_t)sqrt(ix_*ix_ + iy_*iy_);
00269 
00270     //double Crystal_eta = cell->getPosition().eta();
00271     //double Crystal_phi = cell->getPosition().phi();
00272     double Crystal_x = cell->getPosition().x();
00273     double Crystal_y = cell->getPosition().y();
00274     double Crystal_z = cell->getPosition().z();
00275     // ECAL -endcap
00276     if (Crystal_zside == -1)
00277       {
00278         hEEmZ_ix_iy_irMap->setBinContent(Crystal_ix, Crystal_iy, ir);
00279         hEEmZ_ix_iy_xMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_x);
00280         hEEmZ_ix_iy_yMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_y);
00281         hEEmZ_ix_iy_zMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_z);
00282       }
00283     // ECAL +endcap
00284     if (Crystal_zside == 1)
00285       {
00286         hEEpZ_ix_iy_irMap->setBinContent(Crystal_ix, Crystal_iy, ir);
00287         hEEpZ_ix_iy_xMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_x);
00288         hEEpZ_ix_iy_yMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_y);
00289         hEEpZ_ix_iy_zMap->setBinContent(Crystal_ix, Crystal_iy, Crystal_z);
00290       }
00291 
00292     DEBUG( " Crystal " << n );
00293     DEBUG( "  side = " << Crystal_zside );
00294     DEBUG("   ix, iy = " << Crystal_ix << ", " << Crystal_iy);
00295     DEBUG("    x,  y = " << Crystal_x << ", " << Crystal_y);;
00296     DEBUG( " " );
00297 
00298   }
00299  
00300   //-------Set the cell size for each (ieta, iphi) bin-------//
00301   double currentLowEdge_eta = 0;
00302   //double currentHighEdge_eta = 0;
00303   for (int ieta=1; ieta<=85 ; ieta++)
00304     {
00305       int ieta_ = 86 + ieta;
00306       
00307       double eta = hEB_ieta_iphi_etaMap->getBinContent(ieta_, 1);
00308       double etam1 = -999;
00309       
00310       if (ieta==1) 
00311         etam1 = hEB_ieta_iphi_etaMap->getBinContent(85, 1);
00312       else 
00313         etam1 = hEB_ieta_iphi_etaMap->getBinContent(ieta_ - 1, 1);
00314 
00315       //double phi = hEB_ieta_iphi_phiMap->getBinContent(ieta_, 1);
00316       double deta = fabs( eta - etam1 );
00317       double dphi = fabs( hEB_ieta_iphi_phiMap->getBinContent(ieta_, 1) - hEB_ieta_iphi_phiMap->getBinContent(ieta_, 2) );
00318           
00319       currentLowEdge_eta += deta;
00320       hEB_ieta_detaMap->setBinContent(ieta_, deta); // positive rings
00321       hEB_ieta_dphiMap->setBinContent(ieta_, dphi); // positive rings
00322       hEB_ieta_detaMap->setBinContent(86-ieta, deta); // negative rings
00323       hEB_ieta_dphiMap->setBinContent(86-ieta, dphi); // negative rings
00324     }
00325 }
00326 
00327 
00328 
00329 void ECALRecHitAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00330 {
00331   CurrentEvent++;
00332   DEBUG( "Event: " << CurrentEvent);
00333   WriteECALRecHits( iEvent, iSetup );
00334   hECAL_Nevents->Fill(0.5);
00335 }
00336 
00337 void ECALRecHitAnalyzer::WriteECALRecHits(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00338 {
00339   edm::Handle<EBRecHitCollection> EBRecHits;
00340   edm::Handle<EERecHitCollection> EERecHits;
00341   iEvent.getByLabel( EBRecHitsLabel_, EBRecHits );
00342   iEvent.getByLabel( EERecHitsLabel_, EERecHits );
00343   DEBUG( "Got ECALRecHits");
00344 
00345   /*
00346     edm::Handle<reco::CandidateCollection> to;
00347     iEvent.getByLabel( "caloTowers", to );
00348     const CandidateCollection *towers = (CandidateCollection *)to.product();
00349     reco::CandidateCollection::const_iterator tower = towers->begin();
00350     edm::Ref<CaloTowerCollection> towerRef = tower->get<CaloTowerRef>();
00351     const CaloTowerCollection *towerCollection = towerRef.product();
00352     CaloTowerCollection::const_iterator calotower = towerCollection->begin();
00353     DEBUG( "Got Towers");    
00354     DEBUG( "tower size = " << towerCollection->size());
00355   */
00356 
00357   edm::ESHandle<CaloGeometry> pG;
00358   iSetup.get<CaloGeometryRecord>().get(pG);
00359   const CaloGeometry cG = *pG;
00360   //const CaloSubdetectorGeometry* EBgeom=cG.getSubdetectorGeometry(DetId::Ecal,1);
00361   //const CaloSubdetectorGeometry* EEgeom=cG.getSubdetectorGeometry(DetId::Ecal,2);
00362   DEBUG( "Got Geometry");
00363 
00364   TLorentzVector vEBMET_EtaRing[171];
00365   int EBActiveRing[171];
00366   int EBNActiveCells[171];
00367   double EBSET_EtaRing[171];
00368   double EBMaxEnergy_EtaRing[171];
00369   double EBMinEnergy_EtaRing[171];
00370   double EBenergy_EtaRing[171];
00371 
00372   for (int i=0; i<171; i++)
00373     {
00374       EBActiveRing[i] = 0;
00375       EBNActiveCells[i] = 0;
00376       EBSET_EtaRing[i] = 0.0;
00377       EBMaxEnergy_EtaRing[i] = -999; 
00378       EBMinEnergy_EtaRing[i] = 14E3; 
00379       EBenergy_EtaRing[i] = 0.0;
00380     }
00381 
00382   edm::LogInfo("OutputInfo") << "Looping over EB" << std::endl;
00383 
00384   EBRecHitCollection::const_iterator ebrechit;
00385   //int nEBrechit = 0;
00386 
00387   for (ebrechit = EBRecHits->begin(); ebrechit != EBRecHits->end(); ebrechit++) {
00388     
00389     EBDetId det = ebrechit->id();
00390     double Energy = ebrechit->energy();
00391     Int_t ieta = det.ieta();
00392     Int_t iphi = det.iphi();
00393     int EtaRing = 85 + ieta; // this counts from 0
00394     double eta = hEB_ieta_iphi_etaMap->getBinContent(EtaRing+1,iphi);
00395     double phi = hEB_ieta_iphi_phiMap->getBinContent(EtaRing+1,iphi);
00396     double theta = 2*TMath::ATan(exp(-1*eta));
00397     double ET = Energy*TMath::Sin(theta);
00398     TLorentzVector v_;
00399 
00400     if (Energy>EBMaxEnergy_EtaRing[EtaRing]) EBMaxEnergy_EtaRing[EtaRing] = Energy;
00401     if (Energy<EBMinEnergy_EtaRing[EtaRing]) EBMinEnergy_EtaRing[EtaRing] = Energy;
00402     
00403     if (Energy>0)
00404       {
00405         EBActiveRing[EtaRing] = 1;
00406         EBNActiveCells[EtaRing]++;
00407         EBSET_EtaRing[EtaRing]+=ET;
00408         v_.SetPtEtaPhiE(ET, 0, phi, ET);
00409         vEBMET_EtaRing[EtaRing]-=v_;
00410         EBenergy_EtaRing[EtaRing]+=Energy;
00411         hEB_Occ_ieta_iphi->Fill(ieta, iphi);
00412 
00413       }
00414     
00415     hEB_energy_ieta_iphi->Fill(ieta, iphi, Energy);
00416     if (Energy>hEB_Maxenergy_ieta_iphi->getBinContent(EtaRing+1, iphi))
00417       hEB_Maxenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00418     if (Energy<hEB_Minenergy_ieta_iphi->getBinContent(EtaRing+1, iphi))
00419       hEB_Minenergy_ieta_iphi->setBinContent(EtaRing+1, iphi, Energy);
00420 
00421 
00422   } // loop over EB
00423 
00424   for (int iEtaRing = 0; iEtaRing < 171; iEtaRing++)
00425     {
00426       hEB_Minenergyvsieta->Fill(iEtaRing-85, EBMinEnergy_EtaRing[iEtaRing]);
00427       hEB_Maxenergyvsieta->Fill(iEtaRing-85, EBMaxEnergy_EtaRing[iEtaRing]);
00428 
00429       if (EBActiveRing[iEtaRing])
00430         {
00431           hEB_METvsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Pt());
00432           hEB_METPhivsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Phi());
00433           hEB_MExvsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Px());
00434           hEB_MEyvsieta->Fill(iEtaRing-85, vEBMET_EtaRing[iEtaRing].Py());
00435           hEB_SETvsieta->Fill(iEtaRing-85, EBSET_EtaRing[iEtaRing]);
00436           hEB_Occvsieta->Fill(iEtaRing-85, EBNActiveCells[iEtaRing]);
00437           hEB_energyvsieta->Fill(iEtaRing-85, EBenergy_EtaRing[iEtaRing]);
00438         }
00439     }
00440 
00441 
00442   TLorentzVector vEEpZMET_EtaRing[101];
00443   int EEpZActiveRing[101];
00444   int EEpZNActiveCells[101];
00445   double EEpZSET_EtaRing[101];
00446   double EEpZMaxEnergy_EtaRing[101];
00447   double EEpZMinEnergy_EtaRing[101];
00448 
00449   TLorentzVector vEEmZMET_EtaRing[101];
00450   int EEmZActiveRing[101];
00451   int EEmZNActiveCells[101];
00452   double EEmZSET_EtaRing[101];
00453   double EEmZMaxEnergy_EtaRing[101];
00454   double EEmZMinEnergy_EtaRing[101];
00455 
00456   for (int i=0;i<101; i++)
00457     {
00458       EEpZActiveRing[i] = 0;
00459       EEpZNActiveCells[i] = 0;
00460       EEpZSET_EtaRing[i] = 0.0;
00461       EEpZMaxEnergy_EtaRing[i] = -999;
00462       EEpZMinEnergy_EtaRing[i] = 14E3;
00463 
00464       EEmZActiveRing[i] = 0;
00465       EEmZNActiveCells[i] = 0;
00466       EEmZSET_EtaRing[i] = 0.0;
00467       EEmZMaxEnergy_EtaRing[i] = -999;
00468       EEmZMinEnergy_EtaRing[i] = 14E3;
00469     }
00470 
00471   edm::LogInfo("OutputInfo") << "Looping over EE" << std::endl;
00472   EERecHitCollection::const_iterator eerechit;
00473   //int nEErechit = 0;
00474   for (eerechit = EERecHits->begin(); eerechit != EERecHits->end(); eerechit++) {
00475     
00476     EEDetId det = eerechit->id();
00477     double Energy = eerechit->energy();
00478     Int_t ix = det.ix();
00479     Int_t iy = det.iy();
00480     //Float_t ix_ = (Float_t)-999;
00481     //Float_t iy_ = (Float_t)-999;
00482     Int_t ir = -999;
00483     //    edm::LogInfo("OutputInfo") << ix << " " << iy << " " << ix_ << " " << iy_ << " " << ir << std::endl;
00484 
00485     double x = -999;
00486     double y = -999;
00487     double z = -999;
00488     double eta = -999;
00489     double theta = -999;
00490     double phi = -999;
00491 
00492     int Crystal_zside = det.zside();
00493 
00494     if (Crystal_zside == -1)
00495       {
00496         ir = (Int_t)hEEmZ_ix_iy_irMap->getBinContent(ix,iy);
00497         x = hEEmZ_ix_iy_xMap->getBinContent(ix,iy);
00498         y = hEEmZ_ix_iy_yMap->getBinContent(ix,iy);
00499         z = hEEmZ_ix_iy_zMap->getBinContent(ix,iy);
00500       }
00501     if (Crystal_zside == 1)
00502       {
00503         ir = (Int_t)hEEpZ_ix_iy_irMap->getBinContent(ix,iy);
00504         x = hEEpZ_ix_iy_xMap->getBinContent(ix,iy);
00505         y = hEEpZ_ix_iy_yMap->getBinContent(ix,iy);
00506         z = hEEpZ_ix_iy_zMap->getBinContent(ix,iy);
00507       }
00508     TVector3 pos_vector(x,y,z);
00509     phi = pos_vector.Phi();
00510     theta = pos_vector.Theta();
00511     eta = pos_vector.Eta();
00512     double ET = Energy*TMath::Sin(theta);
00513     TLorentzVector v_;
00514 
00515 
00516     if (Crystal_zside == -1)
00517       {
00518         if (Energy>0)
00519           {
00520             EEmZActiveRing[ir] = 1;
00521             EEmZNActiveCells[ir]++;
00522             EEmZSET_EtaRing[ir]+=ET;
00523             v_.SetPtEtaPhiE(ET,0,phi,ET);
00524             vEEmZMET_EtaRing[ir]-=v_;
00525             hEEmZ_Occ_ix_iy->Fill(ix, iy);
00526           }
00527         hEEmZ_energyvsir->Fill(ir, Energy);
00528         hEEmZ_energy_ix_iy->Fill(ix, iy, Energy);
00529 
00530         if (Energy>EEmZMaxEnergy_EtaRing[ir]) EEmZMaxEnergy_EtaRing[ir] = Energy;
00531         if (Energy<EEmZMinEnergy_EtaRing[ir]) EEmZMinEnergy_EtaRing[ir] = Energy;
00532 
00533         if (Energy>hEEmZ_Maxenergy_ix_iy->getBinContent(ix,iy))
00534           hEEmZ_Maxenergy_ix_iy->setBinContent(ix,iy, Energy);
00535         if (Energy<hEEmZ_Minenergy_ix_iy->getBinContent(ix,iy))
00536           hEEmZ_Minenergy_ix_iy->setBinContent(ix,iy, Energy);
00537       }
00538     if (Crystal_zside == 1)
00539       {
00540         if (Energy>0)
00541           {
00542             EEpZActiveRing[ir] = 1;
00543             EEpZNActiveCells[ir]++;
00544             EEpZSET_EtaRing[ir]+=ET;
00545             v_.SetPtEtaPhiE(ET,0,phi,ET);
00546             vEEpZMET_EtaRing[ir]-=v_;
00547             hEEpZ_Occ_ix_iy->Fill(ix, iy);
00548           }
00549         hEEpZ_energyvsir->Fill(ir, Energy);
00550         hEEpZ_energy_ix_iy->Fill(ix, iy, Energy);
00551 
00552         if (Energy>EEpZMaxEnergy_EtaRing[ir]) EEpZMaxEnergy_EtaRing[ir] = Energy;
00553         if (Energy<EEpZMinEnergy_EtaRing[ir]) EEpZMinEnergy_EtaRing[ir] = Energy;
00554         if (Energy>hEEpZ_Maxenergy_ix_iy->getBinContent(ix,iy))
00555           hEEpZ_Maxenergy_ix_iy->setBinContent(ix,iy, Energy);
00556         if (Energy<hEEpZ_Minenergy_ix_iy->getBinContent(ix,iy))
00557           hEEpZ_Minenergy_ix_iy->setBinContent(ix,iy, Energy);
00558       }
00559   } // loop over EE
00560   edm::LogInfo("OutputInfo") << "Done Looping over EE" << std::endl;
00561   for (int iEtaRing = 0; iEtaRing<101; iEtaRing++)
00562     {
00563       hEEpZ_Maxenergyvsir->Fill(iEtaRing, EEpZMaxEnergy_EtaRing[iEtaRing]);
00564       hEEpZ_Minenergyvsir->Fill(iEtaRing, EEpZMinEnergy_EtaRing[iEtaRing]);
00565       hEEmZ_Maxenergyvsir->Fill(iEtaRing, EEmZMaxEnergy_EtaRing[iEtaRing]);
00566       hEEmZ_Minenergyvsir->Fill(iEtaRing, EEmZMinEnergy_EtaRing[iEtaRing]);
00567 
00568       if (EEpZActiveRing[iEtaRing])
00569         {
00570           hEEpZ_METvsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Pt());
00571           hEEpZ_METPhivsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Phi());
00572           hEEpZ_MExvsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Px());
00573           hEEpZ_MEyvsir->Fill(iEtaRing, vEEpZMET_EtaRing[iEtaRing].Py());
00574           hEEpZ_SETvsir->Fill(iEtaRing, EEpZSET_EtaRing[iEtaRing]);
00575           hEEpZ_Occvsir->Fill(iEtaRing, EEpZNActiveCells[iEtaRing]);
00576         }
00577 
00578       if (EEmZActiveRing[iEtaRing])
00579         {
00580           hEEmZ_METvsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Pt());
00581           hEEmZ_METPhivsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Phi());
00582           hEEmZ_MExvsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Px());
00583           hEEmZ_MEyvsir->Fill(iEtaRing, vEEmZMET_EtaRing[iEtaRing].Py());
00584           hEEmZ_SETvsir->Fill(iEtaRing, EEmZSET_EtaRing[iEtaRing]);
00585           hEEmZ_Occvsir->Fill(iEtaRing, EEmZNActiveCells[iEtaRing]);
00586         }
00587     }
00588   edm::LogInfo("OutputInfo") << "Done ..." << std::endl;
00589 } // loop over RecHits
00590 
00591 
00592
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/DQMOffline/JetMET/src/ECALRecHitAnalyzer.cc
