/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/DQM/HcalMonitorTasks/src/HcalBeamMonitor.cc

Go to the documentation of this file.

#include "DQM/HcalMonitorTasks/interface/HcalBeamMonitor.h"
#include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
#include "CondFormats/DataRecord/interface/HcalChannelQualityRcd.h"
#include "Geometry/HcalTowerAlgo/src/HcalHardcodeGeometryData.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/LuminosityBlock.h"
#include "FWCore/Framework/interface/Run.h"

#include <iomanip>
#include <cmath>

// define sizes of ieta arrays for each subdetector

#define PI        3.1415926535897932
#define HBETASIZE 34  // one more bin than needed, I think
#define HEETASIZE 60  // ""
#define HOETASIZE 32  // ""
#define HFETASIZE 84  // ""

/*  Task calculates various moments of Hcal recHits 

    v1.0
    16 August 2008
    by Jeff Temple

*/

// constructor
HcalBeamMonitor::HcalBeamMonitor(const edm::ParameterSet& ps):
  ETA_OFFSET_HB(16),
  ETA_OFFSET_HE(29),
  ETA_BOUND_HE(17),
  ETA_OFFSET_HO(15),
  ETA_OFFSET_HF(41),
  ETA_BOUND_HF(29)
{
  Online_                = ps.getUntrackedParameter<bool>("online",false);
  mergeRuns_             = ps.getUntrackedParameter<bool>("mergeRuns",false);
  enableCleanup_         = ps.getUntrackedParameter<bool>("enableCleanup",false);
  debug_                 = ps.getUntrackedParameter<int>("debug",0);
  prefixME_              = ps.getUntrackedParameter<std::string>("subSystemFolder","Hcal/");
  if (prefixME_.substr(prefixME_.size()-1,prefixME_.size())!="/")
    prefixME_.append("/");
  subdir_                = ps.getUntrackedParameter<std::string>("TaskFolder","BeamMonitor_Hcal");
  if (subdir_.size()>0 && subdir_.substr(subdir_.size()-1,subdir_.size())!="/")
    subdir_.append("/");
  subdir_=prefixME_+subdir_;
  AllowedCalibTypes_     = ps.getUntrackedParameter<std::vector<int> > ("AllowedCalibTypes");
  skipOutOfOrderLS_      = ps.getUntrackedParameter<bool>("skipOutOfOrderLS",true);
  NLumiBlocks_           = ps.getUntrackedParameter<int>("NLumiBlocks",4000);
  makeDiagnostics_       = ps.getUntrackedParameter<bool>("makeDiagnostics",false);

  // Beam Monitor-specific stuff

  // Collection type info
  digiLabel_             =ps.getUntrackedParameter<edm::InputTag>("digiLabel");
  hbheRechitLabel_       = ps.getUntrackedParameter<edm::InputTag>("hbheRechitLabel");
  hoRechitLabel_         = ps.getUntrackedParameter<edm::InputTag>("hoRechitLabel");
  hfRechitLabel_         = ps.getUntrackedParameter<edm::InputTag>("hfRechitLabel");

  // minimum events required in lumi block for tests to be processed
  minEvents_       = ps.getUntrackedParameter<int>("minEvents",500); 
  lumiqualitydir_ = ps.getUntrackedParameter<std::string>("lumiqualitydir","");
  if (lumiqualitydir_.size()>0 && lumiqualitydir_.substr(lumiqualitydir_.size()-1,lumiqualitydir_.size())!="/")
    lumiqualitydir_.append("/");
  occThresh_ = ps.getUntrackedParameter<double>("occupancyThresh",0.0625);  // energy required to be counted by dead/hot checks
  hotrate_        = ps.getUntrackedParameter<double>("hotrate",0.25);
  minBadCells_    = ps.getUntrackedParameter<int>("minBadCells",10);
  Overwrite_      = ps.getUntrackedParameter<bool>("Overwrite",false);
  setupDone_      = false;
}



HcalBeamMonitor::~HcalBeamMonitor() {}

void HcalBeamMonitor::reset() 
{
  CenterOfEnergyRadius->Reset();
  CenterOfEnergy->Reset();
  COEradiusVSeta->Reset();

  HBCenterOfEnergyRadius->Reset();
  HBCenterOfEnergy->Reset();
  HECenterOfEnergyRadius->Reset();
  HECenterOfEnergy->Reset();
  HOCenterOfEnergyRadius->Reset();
  HOCenterOfEnergy->Reset();
  HFCenterOfEnergyRadius->Reset();
  HFCenterOfEnergy->Reset();

  Etsum_eta_L->Reset();
  Etsum_eta_S->Reset();
  Etsum_phi_L->Reset();
  Etsum_phi_S->Reset();
  Etsum_ratio_p->Reset();
  Etsum_ratio_m->Reset();
  Etsum_map_L->Reset();
  Etsum_map_S->Reset();
  Etsum_ratio_map->Reset();
  Etsum_rphi_L->Reset();
  Etsum_rphi_S->Reset();
  Energy_Occ->Reset();

  Occ_rphi_L->Reset();
  Occ_rphi_S->Reset();
  Occ_eta_L->Reset();
  Occ_eta_S->Reset();
  Occ_phi_L->Reset();
  Occ_phi_S->Reset();
  Occ_map_L->Reset();
  Occ_map_S->Reset();
  
  HFlumi_ETsum_perwedge->Reset();
  HFlumi_Occupancy_above_thr_r1->Reset();
  HFlumi_Occupancy_between_thrs_r1->Reset();
  HFlumi_Occupancy_below_thr_r1->Reset();
  HFlumi_Occupancy_above_thr_r2->Reset();
  HFlumi_Occupancy_between_thrs_r2->Reset();
  HFlumi_Occupancy_below_thr_r2->Reset();

  HFlumi_Occupancy_per_channel_vs_lumiblock_RING1->Reset();
  HFlumi_Occupancy_per_channel_vs_lumiblock_RING2->Reset();
  HFlumi_Occupancy_per_channel_vs_BX_RING1->Reset();
  HFlumi_Occupancy_per_channel_vs_BX_RING2->Reset();
  HFlumi_ETsum_vs_BX->Reset();
  HFlumi_Et_per_channel_vs_lumiblock->Reset();

  HFlumi_occ_LS->Reset();
  HFlumi_total_hotcells->Reset();
  HFlumi_total_deadcells->Reset();
  HFlumi_diag_hotcells->Reset();
  HFlumi_diag_deadcells->Reset();


  HFlumi_Ring1Status_vs_LS->Reset();
  HFlumi_Ring2Status_vs_LS->Reset();
}

void HcalBeamMonitor::cleanup()
{
  if (dbe_) 
    {
      dbe_->setCurrentFolder(subdir_);
      dbe_->removeContents();
      dbe_->setCurrentFolder(subdir_+"LSvalues");
      dbe_->removeContents();
    } // if (dbe_)
} // void HcalBeamMonitor::cleanup()


void HcalBeamMonitor::setup()
{
  if (setupDone_)
    return;
  setupDone_ = true;
   if (debug_>0) std::cout <<"<HcalBeamMonitor::setup> Setup in progress..."<<std::endl;
  HcalBaseDQMonitor::setup();
  if (!dbe_) return;

  //jason's
  dbe_->setCurrentFolder(subdir_);
  CenterOfEnergyRadius = dbe_->book1D("CenterOfEnergyRadius",
                                       "Center Of Energy radius",
                                       200,0,1);
      
  CenterOfEnergyRadius->setAxisTitle("(normalized) radius",1);
      
  CenterOfEnergy = dbe_->book2D("CenterOfEnergy",
                                 "Center of Energy;normalized x coordinate;normalize y coordinate",
                                 40,-1,1,
                                 40,-1,1);

  COEradiusVSeta = dbe_->bookProfile("COEradiusVSeta",
                                      "Center of Energy radius vs i#eta",
                                      172,-43,43,
                                      20,0,1);
  COEradiusVSeta->setAxisTitle("i#eta",1);
  COEradiusVSeta->setAxisTitle("(normalized) radius",2);
      
  std::stringstream histname;
  std::stringstream histtitle;
  dbe_->setCurrentFolder(subdir_+"HB");
  HBCenterOfEnergyRadius = dbe_->book1D("HBCenterOfEnergyRadius",
                                         "HB Center Of Energy radius",
                                         200,0,1);
  HBCenterOfEnergy = dbe_->book2D("HBCenterOfEnergy",
                                   "HB Center of Energy",
                                   40,-1,1,
                                   40,-1,1);
  if (makeDiagnostics_)
    {
      for (int i=-16;i<=16;++i)
        {
          if (i==0) continue;
          histname.str("");
          histtitle.str("");
          histname<<"HB_CenterOfEnergyRadius_ieta"<<i;
          histtitle<<"HB Center Of Energy ieta = "<<i;
          HB_CenterOfEnergyRadius[i+ETA_OFFSET_HB]=dbe_->book1D(histname.str().c_str(),
                                                                 histtitle.str().c_str(),
                                                                 200,0,1);
        } // end of HB loop
    }
  dbe_->setCurrentFolder(subdir_+"HE");
  HECenterOfEnergyRadius = dbe_->book1D("HECenterOfEnergyRadius",
                                         "HE Center Of Energy radius",
                                         200,0,1);
  HECenterOfEnergy = dbe_->book2D("HECenterOfEnergy",
                                   "HE Center of Energy",
                                   40,-1,1,
                                   40,-1,1);

  if (makeDiagnostics_)
    {
      for (int i=-29;i<=29;++i)
        {
          if (abs(i)<ETA_BOUND_HE) continue;
          histname.str("");
          histtitle.str("");
          histname<<"HE_CenterOfEnergyRadius_ieta"<<i;
          histtitle<<"HE Center Of Energy ieta = "<<i;
          HE_CenterOfEnergyRadius[i+ETA_OFFSET_HE]=dbe_->book1D(histname.str().c_str(),
                                                                 histtitle.str().c_str(),
                                                                 200,0,1);
        } // end of HE loop
    }
  dbe_->setCurrentFolder(subdir_+"HO");
  HOCenterOfEnergyRadius = dbe_->book1D("HOCenterOfEnergyRadius",
                                         "HO Center Of Energy radius",
                                         200,0,1);
  HOCenterOfEnergy = dbe_->book2D("HOCenterOfEnergy",
                                   "HO Center of Energy",
                                   40,-1,1,
                                   40,-1,1);
  if (makeDiagnostics_)
    {
      for (int i=-15;i<=15;++i)
        {
          if (i==0) continue;
          histname.str("");
          histtitle.str("");
          histname<<"HO_CenterOfEnergyRadius_ieta"<<i;
          histtitle<<"HO Center Of Energy radius ieta = "<<i;
          HO_CenterOfEnergyRadius[i+ETA_OFFSET_HO]=dbe_->book1D(histname.str().c_str(),
                                                                 histtitle.str().c_str(),
                                                                 200,0,1);
        } // end of HO loop
    }
  dbe_->setCurrentFolder(subdir_+"HF");
  HFCenterOfEnergyRadius = dbe_->book1D("HFCenterOfEnergyRadius",
                                         "HF Center Of Energy radius",
                                         200,0,1);
  HFCenterOfEnergy = dbe_->book2D("HFCenterOfEnergy",
                                   "HF Center of Energy",
                                   40,-1,1,
                                   40,-1,1);
  if (makeDiagnostics_)
    {
      for (int i=-41;i<=41;++i)
        {
          if (abs(i)<ETA_BOUND_HF) continue;
          histname.str("");
          histtitle.str("");
          histname<<"HF_CenterOfEnergyRadius_ieta"<<i;
          histtitle<<"HF Center Of Energy radius ieta = "<<i;
          HF_CenterOfEnergyRadius[i+ETA_OFFSET_HF]=dbe_->book1D(histname.str().c_str(),
                                                                 histtitle.str().c_str(),
                                                                 200,0,1);
        } // end of HF loop
    }
      
  dbe_->setCurrentFolder(subdir_+"Lumi");
  // Wenhan's 
  // reducing bins from ",200,0,2000" to ",40,0,800"
      
  float radiusbins[13]={169,201,240,286,340,406,483,576,686,818,975,1162,1300};
  float phibins[71]={-3.5,-3.4,-3.3,-3.2,-3.1,
                     -3.0,-2.9,-2.8,-2.7,-2.6,-2.5,-2.4,-2.3,-2.2,-2.1,
                     -2.0,-1.9,-1.8,-1.7,-1.6,-1.5,-1.4,-1.3,-1.2,-1.1,
                     -1.0,-0.9,-0.8,-0.7,-0.6,-0.5,-0.4,-0.3,-0.2,-0.1,
                     0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
                     1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9,
                     2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9,
                     3.0, 3.1, 3.2, 3.3, 3.4, 3.5};
  Etsum_eta_L=dbe_->bookProfile("Et Sum vs Eta Long Fiber","Et Sum per Area vs Eta Long Fiber",27,0,27,100,0,100);
  Etsum_eta_S=dbe_->bookProfile("Et Sum vs Eta Short Fiber","Et Sum per Area vs Eta Short Fiber",27,0,27,100,0,100);
  Etsum_phi_L=dbe_->bookProfile("Et Sum vs Phi Long Fiber","Et Sum per Area vs Phi Long Fiber",36,0.5,72.5,100,0,100);
  Etsum_phi_S=dbe_->bookProfile("Et Sum vs Phi Short Fiber","Et Sum per Area crossing vs Phi Short Fiber",36,0.5,72.5,100,0,100);

  Etsum_ratio_p=dbe_->book1D("Occ vs PMT events HF+","Energy difference of Long and Short Fiber HF+ in PMT events",105,0.,1.05);
  Energy_Occ=dbe_->book1D("Occ vs Energy","Occupancy vs Energy",200,0,2000);
  Etsum_ratio_m=dbe_->book1D("Occ vs PMT events HF-","Energy difference of Long and Short Fiber HF- in PMT events",105,0.,1.05);
  Etsum_map_L=dbe_->book2D("EtSum 2D phi and eta Long Fiber","Et Sum 2D phi and eta Long Fiber",27,0,27,36,0.5,72.5);
  Etsum_map_S=dbe_->book2D("EtSum 2D phi and eta Short Fiber","Et Sum 2D phi and eta Short Fiber",27,0,27,36,0.5,72.5);

  Etsum_rphi_S=dbe_->book2D("EtSum 2D phi and radius Short Fiber","Et Sum 2D phi and radius Short Fiber",12, radiusbins, 70, phibins);
  Etsum_rphi_L=dbe_->book2D("EtSum 2D phi and radius Long Fiber","Et Sum 2D phi and radius Long Fiber",12, radiusbins, 70, phibins);

  Etsum_ratio_map=dbe_->book2D("Abnormal PMT events","Abnormal PMT events",
                                8,0,8,36, 0.5,72.5);
  SetEtaLabels(Etsum_ratio_map);

  HFlumi_occ_LS = dbe_->book2D("HFlumi_occ_LS","HFlumi occupancy for current LS",
                                8,0,8,36, 0.5,72.5);
  SetEtaLabels(HFlumi_occ_LS);
      
  HFlumi_total_deadcells = dbe_->book2D("HFlumi_total_deadcells","Number of dead lumi channels for LS with at least 10 bad channels",
                                         8,0,8,36,0.5,72.5);
  SetEtaLabels(HFlumi_total_deadcells);
  HFlumi_total_hotcells = dbe_->book2D("HFlumi_total_hotcells","Number of hot lumi channels for LS with at least 10 bad channels",
                                        8,0,8,36,0.5,72.5);
  SetEtaLabels(HFlumi_total_hotcells);

  HFlumi_diag_deadcells = dbe_->book2D("HFlumi_diag_deadcells","Channels that had no hit for at least one LS",
                                       8,0,8,36,0.5,72.5);
  SetEtaLabels(HFlumi_diag_deadcells);
  HFlumi_diag_hotcells = dbe_->book2D("HFlumi_diag_hotcells","Channels that appeared hot for at least one LS",
                                      8,0,8,36,0.5,72.5);
  SetEtaLabels(HFlumi_diag_hotcells);



  Occ_rphi_S=dbe_->book2D("Occ 2D phi and radius Short Fiber","Occupancy 2D phi and radius Short Fiber",12, radiusbins, 70, phibins);
  Occ_rphi_L=dbe_->book2D("Occ 2D phi and radius Long Fiber","Occupancy 2D phi and radius Long Fiber",12, radiusbins, 70, phibins);
  Occ_eta_S=dbe_->bookProfile("Occ vs iEta Short Fiber","Occ per Bunch crossing vs iEta Short Fiber",27,0,27,40,0,800);
  Occ_eta_L=dbe_->bookProfile("Occ vs iEta Long Fiber","Occ per Bunch crossing vs iEta Long Fiber",27,0,27,40,0,800);
      
  Occ_phi_L=dbe_->bookProfile("Occ vs iPhi Long Fiber","Occ per Bunch crossing vs iPhi Long Fiber",36,0.5,72.5,40,0,800);
      
  Occ_phi_S=dbe_->bookProfile("Occ vs iPhi Short Fiber","Occ per Bunch crossing vs iPhi Short Fiber",36,0.5,72.5,40,0,800);
      
  Occ_map_L=dbe_->book2D("Occ_map Long Fiber","Occ Map long Fiber (above threshold)",27,0,27,36,0.5,72.5);
  Occ_map_S=dbe_->book2D("Occ_map Short Fiber","Occ Map Short Fiber (above threshold)",27,0,27,36,0.5,72.5);

  std::stringstream binlabel;
  for (int zz=0;zz<27;++zz)
    {
      if (zz<13)
        binlabel<<zz-41;
      else if (zz==13)
        binlabel<<"NULL";
      else
        binlabel<<zz+15;
      Occ_eta_S->setBinLabel(zz+1,binlabel.str().c_str());
      Occ_eta_L->setBinLabel(zz+1,binlabel.str().c_str());
      Occ_map_S->setBinLabel(zz+1,binlabel.str().c_str());
      Occ_map_L->setBinLabel(zz+1,binlabel.str().c_str());
      Etsum_eta_S->setBinLabel(zz+1,binlabel.str().c_str());
      Etsum_eta_L->setBinLabel(zz+1,binlabel.str().c_str());
      Etsum_map_S->setBinLabel(zz+1,binlabel.str().c_str());
      Etsum_map_L->setBinLabel(zz+1,binlabel.str().c_str());
      binlabel.str("");
    }

  //HFlumi plots
  HFlumi_ETsum_perwedge =  dbe_->book1D("HF lumi ET-sum per wedge","HF lumi ET-sum per wedge;wedge",36,1,37);
  HFlumi_ETsum_perwedge->getTH1F()->SetMinimum(0); 
      
  HFlumi_Occupancy_above_thr_r1 =  dbe_->book1D("HF lumi Occupancy above threshold ring1","HF lumi Occupancy above threshold ring1;wedge",36,1,37);
  HFlumi_Occupancy_between_thrs_r1 = dbe_->book1D("HF lumi Occupancy between thresholds ring1","HF lumi Occupancy between thresholds ring1;wedge",36,1,37);
  HFlumi_Occupancy_below_thr_r1 = dbe_->book1D("HF lumi Occupancy below threshold ring1","HF lumi Occupancy below threshold ring1;wedge",36,1,37);
  HFlumi_Occupancy_above_thr_r2 = dbe_->book1D("HF lumi Occupancy above threshold ring2","HF lumi Occupancy above threshold ring2;wedge",36,1,37);
  HFlumi_Occupancy_between_thrs_r2 = dbe_->book1D("HF lumi Occupancy between thresholds ring2","HF lumi Occupancy between thresholds ring2;wedge",36,1,37);
  HFlumi_Occupancy_below_thr_r2 = dbe_->book1D("HF lumi Occupancy below threshold ring2","HF lumi Occupancy below threshold ring2;wedge",36,1,37);

  HFlumi_Occupancy_above_thr_r1->getTH1F()->SetMinimum(0);
  HFlumi_Occupancy_between_thrs_r1->getTH1F()->SetMinimum(0);
  HFlumi_Occupancy_below_thr_r1->getTH1F()->SetMinimum(0);
  HFlumi_Occupancy_above_thr_r2->getTH1F()->SetMinimum(0);
  HFlumi_Occupancy_between_thrs_r2->getTH1F()->SetMinimum(0);
  HFlumi_Occupancy_below_thr_r2->getTH1F()->SetMinimum(0);
 
  HFlumi_Occupancy_per_channel_vs_lumiblock_RING1 = dbe_->bookProfile("HFlumiRing1OccupancyPerChannelVsLB",
                                                                      "HFlumi Occupancy per channel vs lumi-block (RING 1);LS; -ln(empty fraction)",
                                                                      NLumiBlocks_,0.5,NLumiBlocks_+0.5,100,0,10000);
  HFlumi_Occupancy_per_channel_vs_lumiblock_RING2 = dbe_->bookProfile("HFlumiRing2OccupancyPerChannelVsLB","HFlumi Occupancy per channel vs lumi-block (RING 2);LS; -ln(empty fraction)",NLumiBlocks_,0.5,NLumiBlocks_+0.5,100,0,10000);

  HFlumi_Occupancy_per_channel_vs_BX_RING1 = dbe_->bookProfile("HFlumi Occupancy per channel vs BX (RING 1)","HFlumi Occupancy per channel vs BX (RING 1);BX; -ln(empty fraction)",4000,0,4000,100,0,10000);
  HFlumi_Occupancy_per_channel_vs_BX_RING2 = dbe_->bookProfile("HFlumi Occupancy per channel vs BX (RING 2)","HFlumi Occupancy per channel vs BX (RING 2);BX; -ln(empty fraction)",4000,0,4000,100,0,10000);
  HFlumi_ETsum_vs_BX = dbe_->bookProfile("HFlumi_ETsum_vs_BX","HFlumi ETsum vs BX; BX; ETsum",4000,0,4000,100,0,10000);

  HFlumi_Et_per_channel_vs_lumiblock = dbe_->bookProfile("HFlumi Et per channel vs lumi-block","HFlumi Et per channel vs lumi-block;LS;ET",NLumiBlocks_,0.5,NLumiBlocks_+0.5,100,0,10000);

  HFlumi_Occupancy_per_channel_vs_lumiblock_RING1->getTProfile()->SetMarkerStyle(20);
  HFlumi_Occupancy_per_channel_vs_lumiblock_RING2->getTProfile()->SetMarkerStyle(20);
  HFlumi_Et_per_channel_vs_lumiblock->getTProfile()->SetMarkerStyle(20);

  HFlumi_Ring1Status_vs_LS = dbe_->bookProfile("HFlumi_Ring1Status_vs_LS","Fraction of good Ring 1 channels vs LS;LS; Fraction of Good Channels",NLumiBlocks_,0.5,NLumiBlocks_+0.5,100,0,10000);
  HFlumi_Ring2Status_vs_LS = dbe_->bookProfile("HFlumi_Ring2Status_vs_LS","Fraction of good Ring 2 channels vs LS;LS; Fraction of Good Channels",NLumiBlocks_,0.5,NLumiBlocks_+0.5,100,0,10000);
  HFlumi_Ring1Status_vs_LS->getTProfile()->SetMarkerStyle(20);
  HFlumi_Ring2Status_vs_LS->getTProfile()->SetMarkerStyle(20);

  return;
}

void HcalBeamMonitor::beginRun(const edm::Run& run, const edm::EventSetup& c)
{
  HcalBaseDQMonitor::beginRun(run,c);  

  if (debug_>1) std::cout <<"HcalBeamMonitor::beginRun"<<std::endl;
  HcalBaseDQMonitor::beginRun(run,c);

  lastProcessedLS_=0;
  runNumber_=run.id().run();
  if (lumiqualitydir_.size()>0 && Online_==true)
    {
      if (Overwrite_==false)
        outfile_ <<lumiqualitydir_<<"HcalHFLumistatus_"<<runNumber_<<".txt";
      else
        outfile_ <<lumiqualitydir_<<"HcalHFLumistatus.txt";
      std::ofstream outStream(outfile_.str().c_str()); // recreate the file, rather than appending to it
      outStream<<"## Run "<<runNumber_<<std::endl;
      outStream<<"## LumiBlock\tRing1Status\t\tRing2Status\t\tGlobalStatus\tNentries"<<std::endl;
      outStream.close();
    }

  // Get expected good channels in run according to channel quality database
  // Get channel quality status info for each run

  // Default number of expected good channels in the run
  ring1totalchannels_=144;
  ring2totalchannels_=144;
  BadCells_.clear(); // remove any old maps
  // Get Channel quality info for the run
  // Exclude bad channels from overall calculation
  edm::ESHandle<HcalChannelQuality> p;
  c.get<HcalChannelQualityRcd>().get(p);
  HcalChannelQuality* chanquality = new HcalChannelQuality(*p.product());
  std::vector<DetId> mydetids = chanquality->getAllChannels();
  
  for (unsigned int i=0;i<mydetids.size();++i)
    {
      if (mydetids[i].det()!=DetId::Hcal) continue;
      HcalDetId id=mydetids[i];
      
      if (id.subdet()!=HcalForward) continue;
      if ((id.depth()==1 && (abs(id.ieta())==33 || abs(id.ieta())==34)) ||
          (id.depth()==2 && (abs(id.ieta())==35 || abs(id.ieta())==36)))
        {
          const HcalChannelStatus* origstatus=chanquality->getValues(id);
          HcalChannelStatus* mystatus=new HcalChannelStatus(origstatus->rawId(),origstatus->getValue());
          if (mystatus->isBitSet(HcalChannelStatus::HcalCellHot)) 
            BadCells_[id]=HcalChannelStatus::HcalCellHot;
          
          else if (mystatus->isBitSet(HcalChannelStatus::HcalCellDead))
            BadCells_[id]=HcalChannelStatus::HcalCellDead;
          
          if (mystatus->isBitSet(HcalChannelStatus::HcalCellHot) || 
              mystatus->isBitSet(HcalChannelStatus::HcalCellDead))
            {
              if (id.depth()==1) --ring1totalchannels_;
              else if (id.depth()==2) --ring2totalchannels_;
            }
          delete mystatus;
        } // if ((id.depth()==1) ...
    } // for (unsigned int i=0;...)
    
  if (tevt_==0) this->setup(); // create all histograms; not necessary if merging runs together
  if (mergeRuns_==false) this->reset(); // call reset at start of all runs

  return;

} // void HcalBeamMonitor::beginRun(const edm::Run& run, const edm::EventSetup& c)


void HcalBeamMonitor::analyze(const edm::Event& e, const edm::EventSetup& c)
{
  if (!IsAllowedCalibType()) return;
  if (LumiInOrder(e.luminosityBlock())==false) return;

  // try to get rechits and digis
  edm::Handle<HFDigiCollection> hf_digi;

  edm::Handle<HBHERecHitCollection> hbhe_rechit;
  edm::Handle<HORecHitCollection> ho_rechit;
  edm::Handle<HFRecHitCollection> hf_rechit;
  
  if (!(e.getByLabel(digiLabel_,hf_digi)))
    {
      edm::LogWarning("HcalBeamMonitor")<< digiLabel_<<" hf_digi not available";
      return;
    }

  if (!(e.getByLabel(hbheRechitLabel_,hbhe_rechit)))
    {
      edm::LogWarning("HcalBeamMonitor")<< hbheRechitLabel_<<" hbhe_rechit not available";
      return;
    }

  if (!(e.getByLabel(hfRechitLabel_,hf_rechit)))
    {
      edm::LogWarning("HcalBeamMonitor")<< hfRechitLabel_<<" hf_rechit not available";
      return;
    }
  if (!(e.getByLabel(hoRechitLabel_,ho_rechit)))
    {
      edm::LogWarning("HcalBeamMonitor")<< hoRechitLabel_<<" ho_rechit not available";
      return;
    }

  //good event; increment counters and process
  HcalBaseDQMonitor::analyze(e,c);
  processEvent(*hbhe_rechit, *ho_rechit, *hf_rechit, *hf_digi, e.bunchCrossing());

} //void HcalBeamMonitor::analyze(const edm::Event& e, const edm::EventSetup& c)


void HcalBeamMonitor::processEvent(const HBHERecHitCollection& hbheHits,
                                   const HORecHitCollection& hoHits,
                                   const HFRecHitCollection& hfHits,
                                   const HFDigiCollection& hf,
                                   int   bunchCrossing
                                   )
  
{ 
  //processEvent loop
  if (!dbe_)
    {
      if (debug_>0) std::cout <<"HcalBeamMonitor::processEvent   DQMStore not instantiated!!!"<<std::endl;
      return;
    }

  HBHERecHitCollection::const_iterator HBHEiter;
  HORecHitCollection::const_iterator HOiter;
  HFRecHitCollection::const_iterator HFiter;

  double totalX=0;
  double totalY=0;
  double totalE=0;

  double HBtotalX=0;
  double HBtotalY=0;
  double HBtotalE=0;
  double HEtotalX=0;
  double HEtotalY=0;
  double HEtotalE=0;
  double HOtotalX=0;
  double HOtotalY=0;
  double HOtotalE=0;
  double HFtotalX=0;
  double HFtotalY=0;
  double HFtotalE=0;
     
  float hitsp[13][36][2];
  float hitsm[13][36][2];
  float hitsp_Et[13][36][2];
  float hitsm_Et[13][36][2];
  
  for(int m=0;m<13;m++){
    for(int n=0;n<36;n++){
      hitsp[m][n][0]=0;
      hitsp[m][n][1]=0; 
      hitsm[m][n][0]=0;
      hitsm[m][n][1]=0;

      hitsp_Et[m][n][0]=0;
      hitsp_Et[m][n][1]=0; 
      hitsm_Et[m][n][0]=0;
      hitsm_Et[m][n][1]=0;
    }
  }

  if(hbheHits.size()>0)
    {
      double HB_weightedX[HBETASIZE]={0.};
      double HB_weightedY[HBETASIZE]={0.};
      double HB_energy[HBETASIZE]={0.};
      
      double HE_weightedX[HEETASIZE]={0.};
      double HE_weightedY[HEETASIZE]={0.};
      double HE_energy[HEETASIZE]={0.};
      
      int ieta, iphi;
      
      for (HBHEiter=hbheHits.begin(); 
           HBHEiter!=hbheHits.end(); 
           ++HBHEiter) 
        { 
          
          // loop over all hits
          if (HBHEiter->energy()<0) continue; // don't consider negative-energy cells
          HcalDetId id(HBHEiter->detid().rawId());
          ieta=id.ieta();
          iphi=id.iphi();
          
          int index=-1;
          if ((HcalSubdetector)(id.subdet())==HcalBarrel)
            {
              HBtotalX+=HBHEiter->energy()*cos(PI*iphi/36.);
              HBtotalY+=HBHEiter->energy()*sin(PI*iphi/36.);
              HBtotalE+=HBHEiter->energy();
              
              index=ieta+ETA_OFFSET_HB;
              if (index<0 || index>= HBETASIZE) continue;
              HB_weightedX[index]+=HBHEiter->energy()*cos(PI*iphi/36.);
              HB_weightedY[index]+=HBHEiter->energy()*sin(PI*iphi/36.);
              HB_energy[index]+=HBHEiter->energy();
            } // if id.subdet()==HcalBarrel
          
          else
            {
              HEtotalX+=HBHEiter->energy()*cos(PI*iphi/36.);
              HEtotalY+=HBHEiter->energy()*sin(PI*iphi/36.);
              HEtotalE+=HBHEiter->energy();
              
              index=ieta+ETA_OFFSET_HE;
              if (index<0 || index>= HEETASIZE) continue;
              HE_weightedX[index]+=HBHEiter->energy()*cos(PI*iphi/36.);
              HE_weightedY[index]+=HBHEiter->energy()*sin(PI*iphi/36.);
              HE_energy[index]+=HBHEiter->energy();
            }
        } // for (HBHEiter=hbheHits.begin()...
          // Fill each histogram
      
      int hbeta=ETA_OFFSET_HB;
      for (int i=-1*hbeta;i<=hbeta;++i)
        {
          if (i==0) continue;
          int index = i+ETA_OFFSET_HB;
          if (index<0 || index>= HBETASIZE) continue;
          if (HB_energy[index]==0) continue;
          double moment=pow(HB_weightedX[index],2)+pow(HB_weightedY[index],2);
          moment=pow(moment,0.5);
          moment/=HB_energy[index];
          if (moment!=0)
            {
              if (makeDiagnostics_) HB_CenterOfEnergyRadius[index]->Fill(moment);
              COEradiusVSeta->Fill(i,moment);
            }
        } // for (int i=-1*hbeta;i<=hbeta;++i)

      int heeta=ETA_OFFSET_HE;
      for (int i=-1*heeta;i<=heeta;++i)
        {
          if (i==0) continue;
          if (i>-1*ETA_BOUND_HE && i <ETA_BOUND_HE) continue;
          int index = i + ETA_OFFSET_HE;
          if (index<0 || index>= HEETASIZE) continue;
          if (HE_energy[index]==0) continue;
          double moment=pow(HE_weightedX[index],2)+pow(HE_weightedY[index],2);
          moment=pow(moment,0.5);
          moment/=HE_energy[index];
          if (moment!=0)
            {
              if (makeDiagnostics_) HE_CenterOfEnergyRadius[index]->Fill(moment);
              COEradiusVSeta->Fill(i,moment);
            }
        } // for (int i=-1*heeta;i<=heeta;++i)

    } // if (hbheHits.size()>0)

  
  // HO loop
  if(hoHits.size()>0)
    {
      double HO_weightedX[HOETASIZE]={0.};
      double HO_weightedY[HOETASIZE]={0.};
      double HO_energy[HOETASIZE]={0.};
      double offset;

      int ieta, iphi;
      for (HOiter=hoHits.begin(); 
           HOiter!=hoHits.end(); 
           ++HOiter) 
        { 
          // loop over all cells
          if (HOiter->energy()<0) continue;  // don't include negative-energy cells?
          HcalDetId id(HOiter->detid().rawId());
          ieta=id.ieta();
          iphi=id.iphi();

          HOtotalX+=HOiter->energy()*cos(PI*iphi/36.);
          HOtotalY+=HOiter->energy()*sin(PI*iphi/36.);
          HOtotalE+=HOiter->energy();

          int index=ieta+ETA_OFFSET_HO;
          if (index<0 || index>= HOETASIZE) continue;
          HO_weightedX[index]+=HOiter->energy()*cos(PI*iphi/36.);
          HO_weightedY[index]+=HOiter->energy()*sin(PI*iphi/36.);
          HO_energy[index]+=HOiter->energy();
        } // for (HOiter=hoHits.begin();...)
          
      for (int i=-1*ETA_OFFSET_HO;i<=ETA_OFFSET_HO;++i)
        {
          if (i==0) continue;
          int index = i + ETA_OFFSET_HO;
          if (index < 0 || index>= HOETASIZE) continue;
          if (HO_energy[index]==0) continue;
          double moment=pow(HO_weightedX[index],2)+pow(HO_weightedY[index],2);
          moment=pow(moment,0.5);
          moment/=HO_energy[index];
          // Shift HO values by 0.5 units in eta relative to HB
          offset = (i>0 ? 0.5: -0.5);
          if (moment!=0)
            {
              if (makeDiagnostics_) HO_CenterOfEnergyRadius[index]->Fill(moment);
              COEradiusVSeta->Fill(i+offset,moment);
            }
        } // for (int i=-1*hoeta;i<=hoeta;++i)
    } // if (hoHits.size()>0)
    
  // HF loop

  Etsum_ratio_map->Fill(-1,-1,1); // fill underflow bin with number of events
  {
    if(hfHits.size()>0)
      {
        double HF_weightedX[HFETASIZE]={0.};
        double HF_weightedY[HFETASIZE]={0.};
        double HF_energy[HFETASIZE]={0.};
        double offset;
        
        // Assume ZS until shown otherwise
        double emptytowersRing1 = ring1totalchannels_;
        double emptytowersRing2 = ring2totalchannels_;
        double ZStowersRing1 = ring1totalchannels_;
        double ZStowersRing2 = ring2totalchannels_;
        
        int ieta, iphi;
        float et,eta,phi,r;

        HFlumi_occ_LS->Fill(-1,-1,1 ); // event counter in occupancy histogram underflow bin
        // set maximum to HFlumi_occ_LS->getBinContent(0,0)?  
        // that won't work -- offline will add multiple histograms, and maximum will get screwed up?
        // No, we can add it here, but we also need a call to setMaximum in the client as well.
        HFlumi_occ_LS->getTH2F()->SetMaximum(HFlumi_occ_LS->getBinContent(0,0));

        double etx=0, ety=0;

        for (HFiter=hfHits.begin(); 
             HFiter!=hfHits.end(); 
             ++HFiter) 
          {  // loop on hfHits
            // If hit present, don't count it as ZS any more
            ieta = HFiter->id().ieta();
            iphi = HFiter->id().iphi();

            int binieta=ieta;
            if (ieta<0) binieta+=41;
            else if (ieta>0) binieta-=15;

            // Count that hit was found in one of the rings used for luminosity calculation.
            // If so, decrease the number of empty channels per ring by 1
            if (abs(ieta)>=33 && abs(ieta)<=36) // luminosity ring check
              {
                // don't subtract away cells that have already been removed as bad
                if (BadCells_.find(HFiter->id())==BadCells_.end()) // bad cell not found
                  {
                    if ((abs(ieta)<35) && HFiter->id().depth()==1) --ZStowersRing1;
                    else if ((abs(ieta)>34) && HFiter->id().depth()==2) -- ZStowersRing2;
                  }
              }

            if (HFiter->energy()<0) continue;  // don't include negative-energy cells?

            eta=theHFEtaBounds[abs(ieta)-29];
            et=HFiter->energy()/cosh(eta)/area[abs(ieta)-29];
            if (abs(ieta)>=33 && abs(ieta)<=36) // Luminosity ring check
              {
                // don't count cells that are below threshold, or that have been marked bad in Chan Stat DB
                if (et>=occThresh_ && BadCells_.find(HFiter->id())==BadCells_.end() ) // minimum ET threshold
                  {
                    if ((abs(ieta)<35) && HFiter->id().depth()==1) --emptytowersRing1;
                    else if ((abs(ieta)>34) && HFiter->id().depth()==2) -- emptytowersRing2;
                  }
              }
            r=radius[abs(ieta)-29];
            if(HFiter->id().iphi()<37)
              phi=HFiter->id().iphi()*0.087266;
            else phi=(HFiter->id().iphi()-72)*0.087266;
           
            
            if (HFiter->id().depth()==1)
              {
                Etsum_eta_L->Fill(binieta,et);
                Etsum_phi_L->Fill(iphi,et);
                Etsum_map_L->Fill(binieta,iphi,et);
                Etsum_rphi_L->Fill(r,phi,et);
              
                if(ieta>0) {
                  hitsp[ieta-29][(HFiter->id().iphi()-1)/2][0]=HFiter->energy();
                  hitsp_Et[ieta-29][(HFiter->id().iphi()-1)/2][0]=et;
                }
                else if(ieta<0) {
                  hitsm[-ieta-29][(HFiter->id().iphi()-1)/2][0]=HFiter->energy(); 
                  hitsm_Et[-ieta-29][(HFiter->id().iphi()-1)/2][0]=et; 
                }
              } // if (HFiter->id().depth()==1)
         
            //Fill 3 histos for Short Fibers :
            if (HFiter->id().depth()==2)
              {
                Etsum_eta_S->Fill(binieta,et);
                Etsum_phi_S->Fill(iphi,et);
                Etsum_rphi_S->Fill(r,phi,et); 
                Etsum_map_S->Fill(binieta,iphi,et);
                if(ieta>0)  
                  {
                    hitsp[ieta-29][(HFiter->id().iphi()-1)/2][1]=HFiter->energy();
                    hitsp_Et[ieta-29][(HFiter->id().iphi()-1)/2][1]=et;
                  }
                else if(ieta<0)  { 
                  hitsm[-ieta-29][(HFiter->id().iphi()-1)/2][1]=HFiter->energy();
                  hitsm_Et[-ieta-29][(HFiter->id().iphi()-1)/2][1]=et;
                }
          
              } // depth()==2
            Energy_Occ->Fill(HFiter->energy()); 
            
            //HF: no non-threshold occupancy map is filled?

            if ((abs(ieta) == 33 || abs(ieta) == 34) && HFiter->id().depth() == 1)
              { 
                etx+=et*cos(PI*iphi/36.);
                ety+=et*sin(PI*iphi/36.);

                HFlumi_Et_per_channel_vs_lumiblock->Fill(currentLS,et);
                if (et>occThresh_)
                  {
                    int etabin=0;
                    if (ieta<0)
                      etabin=36+ieta; // bins 0-3 correspond to ieta = -36, -35, -34, -33
                    else
                      etabin=ieta-29; // bins 4-7 correspond to ieta = 33, 34, 35, 36
                    HFlumi_occ_LS->Fill(etabin,HFiter->id().iphi());
                  }
              }

            else if ((abs(ieta) == 35 || abs(ieta) == 36) && HFiter->id().depth() == 2)
              { 
                etx+=et*cos(PI*iphi/36.);
                ety+=et*sin(PI*iphi/36.);

                HFlumi_Et_per_channel_vs_lumiblock->Fill(currentLS,et);
                if (et>occThresh_)
                  {
                    int etabin=0;
                    if (ieta<0)
                      etabin=36+ieta; // bins 0-3 correspond to ieta = -36, -35, -34, -33
                    else
                      etabin=ieta-29; // bins 4-7 correspond to ieta = 33, 34, 35, 36
                    HFlumi_occ_LS->Fill(etabin,HFiter->id().iphi());
                  }
              }

            // Fill occupancy plots.
            
            int value=0;
            if(et>occThresh_) value=1;

            if (HFiter->id().depth()==1)
              {
                Occ_eta_L->Fill(binieta,value);
                Occ_phi_L->Fill(iphi,value);
                Occ_map_L->Fill(binieta,iphi,value);
                Occ_rphi_L->Fill(r,phi,value);
              }
              
            else if (HFiter->id().depth()==2)
              {
                Occ_eta_S->Fill(binieta,value);
                Occ_phi_S->Fill(iphi,value);
                Occ_map_S->Fill(binieta,iphi,value);
                Occ_rphi_S->Fill(r,phi,value);
              }  
            HcalDetId id(HFiter->detid().rawId());

            HFtotalX+=HFiter->energy()*cos(PI*iphi/36.);
            HFtotalY+=HFiter->energy()*sin(PI*iphi/36.);
            HFtotalE+=HFiter->energy();

            int index=ieta+ETA_OFFSET_HF;
            if (index<0 || index>= HFETASIZE) continue;
            HF_weightedX[index]+=HFiter->energy()*cos(PI*iphi/36.);
            HF_weightedY[index]+=HFiter->energy()*sin(PI*iphi/36.);
            HF_energy[index]+=HFiter->energy();
            
          } // for (HFiter=hfHits.begin();...)
        
        // looped on all HF hits; calculate empty fraction
        //  empty towers  = # of cells with ET < 0.0625 GeV, or cells missing because of ZS
        //  Calculated as :  144 - (# of cells with ET >= 0.0625 GeV)
        //  At some point, allow for calculations when channels are masked (and less than 144 channels expected)

        // Check Ring 1
        double logvalue=-1;
        if (ring1totalchannels_>0)
          {
            if (emptytowersRing1>0)
              logvalue=-1.*log(emptytowersRing1/ring1totalchannels_);
            HFlumi_Occupancy_per_channel_vs_lumiblock_RING1->Fill(currentLS,logvalue);
            HFlumi_Occupancy_per_channel_vs_BX_RING1->Fill(bunchCrossing,logvalue);
          }
        // Check Ring 2
        logvalue=-1;
        if (ring2totalchannels_>0)
          {
            if (emptytowersRing2>0)
              logvalue=-1.*log(emptytowersRing2/ring2totalchannels_);
            HFlumi_Occupancy_per_channel_vs_lumiblock_RING2->Fill(currentLS,logvalue);
            HFlumi_Occupancy_per_channel_vs_BX_RING2->Fill(bunchCrossing,logvalue);
          }

        HFlumi_ETsum_vs_BX->Fill(bunchCrossing,pow(etx*etx+ety*ety,0.5));
        int hfeta=ETA_OFFSET_HF;
        for (int i=-1*hfeta;i<=hfeta;++i)
          {
            if (i==0) continue;
            if (i>-1*ETA_BOUND_HF && i <ETA_BOUND_HF) continue;
            int index = i + ETA_OFFSET_HF;
            if (index<0 || index>= HFETASIZE) continue;
            if (HF_energy[index]==0) continue;
            double moment=pow(HF_weightedX[index],2)+pow(HF_weightedY[index],2);
            moment=pow(moment,0.5);
            moment/=HF_energy[index];
            offset = (i>0 ? 0.5: -0.5);
            if (moment!=0)
              {
                if (makeDiagnostics_) HF_CenterOfEnergyRadius[index]->Fill(moment);
                COEradiusVSeta->Fill(i+offset,moment);
              }
          } // for (int i=-1*hfeta;i<=hfeta;++i)
        float ratiom,ratiop;
          
        for(int i=0;i<13;i++){
          for(int j=0;j<36;j++){
              
            if(hitsp[i][j][0]==hitsp[i][j][1]) continue;
              
            if (hitsp[i][j][0] < 1.2 && hitsp[i][j][1] < 1.8) continue;
            //use only lumi rings
            if (((i+29) < 33) || ((i+29) > 36)) continue;
            ratiop=fabs((fabs(hitsp[i][j][0])-fabs(hitsp[i][j][1]))/(fabs(hitsp[i][j][0])+fabs(hitsp[i][j][1])));
            //cout<<ratiop<<std::endl;
            if ((hitsp_Et[i][j][0] > 5. && hitsp[i][j][1] < 1.8) || (hitsp_Et[i][j][1] > 5. &&  hitsp[i][j][0] < 1.2)){
              Etsum_ratio_p->Fill(ratiop);
              if(abs(ratiop>0.95)) Etsum_ratio_map->Fill(i,2*j+1); // i=4,5,6,7 for HFlumi rings 
            }
          }
        }
          
        for(int p=0;p<13;p++){
          for(int q=0;q<36;q++){
              
            if(hitsm[p][q][0]==hitsm[p][q][1]) continue;

            if (hitsm[p][q][0] < 1.2 && hitsm[p][q][1] < 1.8) continue;
            //use only lumi rings
            if (((p+29) < 33) || ((p+29) > 36)) continue;
            ratiom=fabs((fabs(hitsm[p][q][0])-fabs(hitsm[p][q][1]))/(fabs(hitsm[p][q][0])+fabs(hitsm[p][q][1])));         
            if ((hitsm_Et[p][q][0] > 5. && hitsm[p][q][1] < 1.8) || (hitsm_Et[p][q][1] > 5. && hitsm[p][q][0] < 1.2)){
              Etsum_ratio_m->Fill(ratiom);
              if(abs(ratiom>0.95)) Etsum_ratio_map->Fill(7-p,2*q+1); // p=4,5,6,7 for HFlumi rings
              //p=7:  ieta=-36; p=4:  ieta=-33
            }
          }
        } 
      } // if (hfHits.size()>0)
  
    totalX=HBtotalX+HEtotalX+HOtotalX+HFtotalX;
    totalY=HBtotalY+HEtotalY+HOtotalY+HFtotalY;
    totalE=HBtotalE+HEtotalE+HOtotalE+HFtotalE;

    double moment;
    if (HBtotalE>0)
      {
        moment=pow(HBtotalX*HBtotalX+HBtotalY*HBtotalY,0.5)/HBtotalE;
        HBCenterOfEnergyRadius->Fill(moment);
        HBCenterOfEnergy->Fill(HBtotalX/HBtotalE, HBtotalY/HBtotalE);
      }
    if (HEtotalE>0)
      {
        moment=pow(HEtotalX*HEtotalX+HEtotalY*HEtotalY,0.5)/HEtotalE;
        HECenterOfEnergyRadius->Fill(moment);
        HECenterOfEnergy->Fill(HEtotalX/HEtotalE, HEtotalY/HEtotalE);
      }
    if (HOtotalE>0)
      {
        moment=pow(HOtotalX*HOtotalX+HOtotalY*HOtotalY,0.5)/HOtotalE;
        HOCenterOfEnergyRadius->Fill(moment);
        HOCenterOfEnergy->Fill(HOtotalX/HOtotalE, HOtotalY/HOtotalE);
      }
    if (HFtotalE>0)
      {
        moment=pow(HFtotalX*HFtotalX+HFtotalY*HFtotalY,0.5)/HFtotalE;
        HFCenterOfEnergyRadius->Fill(moment);
        HFCenterOfEnergy->Fill(HFtotalX/HFtotalE, HFtotalY/HFtotalE);
      }
    if (totalE>0)
      {
        moment = pow(totalX*totalX+totalY*totalY,0.5)/totalE;
        CenterOfEnergyRadius->Fill(moment);
        CenterOfEnergy->Fill(totalX/totalE, totalY/totalE);
      }


    
    for (HFDigiCollection::const_iterator j=hf.begin(); j!=hf.end(); j++){
      const HFDataFrame digi = (const HFDataFrame)(*j);
      //  calibs_= cond.getHcalCalibrations(digi.id());  // Old method was made private. 
      //       float en=0;
      //       float ts =0; float bs=0;
      //       int maxi=0; float maxa=0;
      //       for(int i=sigS0_; i<=sigS1_; i++){
      //        if(digi.sample(i).adc()>maxa){maxa=digi.sample(i).adc(); maxi=i;}
      //       }
      //       for(int i=sigS0_; i<=sigS1_; i++){         
      //        float tmp1 =0;   
      //         int j1=digi.sample(i).adc();
      //         tmp1 = (LedMonAdc2fc[j1]+0.5);           
      //        en += tmp1-calibs_.pedestal(digi.sample(i).capid());
      //        if(i>=(maxi-1) && i<=maxi+1){
      //          ts += i*(tmp1-calibs_.pedestal(digi.sample(i).capid()));
      //          bs += tmp1-calibs_.pedestal(digi.sample(i).capid());
      //        }
      //       }

      //---HFlumiplots
      int theTStobeused = 6;
      // will have masking later:
      int mask=1; 
      if(mask!=1) continue;
      //if we want to sum the 10 TS instead of just taking one:
      for (int i=0; i<digi.size(); i++) {
        if (i==theTStobeused) {
          float tmpET =0;
          int jadc=digi.sample(i).adc();
          //NOW LUT used in HLX are only identy LUTs, so Et filled
          //with unlinearised adc, ie tmpET = jadc
          //      tmpET = (adc2fc[jadc]+0.5);
          tmpET = jadc;

          //-find which wedge we are in
          //  ETsum and Occupancy will be summed for both L and S
          if(digi.id().ieta()>28){
            if((digi.id().iphi()==1)||(digi.id().iphi()==71)){
              HFlumi_ETsum_perwedge->Fill(1,tmpET);
              if((digi.id().ieta()==33)||(digi.id().ieta()==34)) {
                if(jadc>100) HFlumi_Occupancy_above_thr_r1->Fill(1,1);
                if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r1->Fill(1,1);
                if(jadc<10) HFlumi_Occupancy_below_thr_r1->Fill(1,1);
              }
              else if((digi.id().ieta()==35)||(digi.id().ieta()==36)) {
                if(jadc>100) HFlumi_Occupancy_above_thr_r2->Fill(1,1);
                if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r2->Fill(1,1);
                if(jadc<10) HFlumi_Occupancy_below_thr_r2->Fill(1,1);
              }
            }
            else {
              for (int iwedge=2; iwedge<19; iwedge++) {
                int itmp=4*(iwedge-1);
                if( (digi.id().iphi()==(itmp+1)) || (digi.id().iphi()==(itmp-1))) {
                  HFlumi_ETsum_perwedge->Fill(iwedge,tmpET);
                  if((digi.id().ieta()==33)||(digi.id().ieta()==34)) {
                    if(jadc>100) HFlumi_Occupancy_above_thr_r1->Fill(iwedge,1);
                    if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r1->Fill(iwedge,1);
                    if(jadc<10) HFlumi_Occupancy_below_thr_r1->Fill(iwedge,1);
                  }
                  else if((digi.id().ieta()==35)||(digi.id().ieta()==36)) {
                    if(jadc>100) HFlumi_Occupancy_above_thr_r2->Fill(iwedge,1);
                    if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r2->Fill(iwedge,1);
                    if(jadc<10) HFlumi_Occupancy_below_thr_r2->Fill(iwedge,1);
                  }
                  iwedge=99;
                }
              }
            }
          }  //--endif ieta in HF+
          else if(digi.id().ieta()<-28){
            if((digi.id().iphi()==1)||(digi.id().iphi()==71)){
              HFlumi_ETsum_perwedge->Fill(19,tmpET);
              if((digi.id().ieta()==-33)||(digi.id().ieta()==-34)) {
                if(jadc>100) HFlumi_Occupancy_above_thr_r1->Fill(19,1);
                if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r1->Fill(19,1);
                if(jadc<10) HFlumi_Occupancy_below_thr_r1->Fill(19,1);
              }
              else if((digi.id().ieta()==-35)||(digi.id().ieta()==-36)) {
                if(jadc>100) HFlumi_Occupancy_above_thr_r2->Fill(19,1);
                if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r2->Fill(19,1);
                if(jadc<10) HFlumi_Occupancy_below_thr_r2->Fill(19,1);
              }
            }
            else {
              for (int iw=2; iw<19; iw++) {
                int itemp=4*(iw-1);
                if( (digi.id().iphi()==(itemp+1)) || (digi.id().iphi()==(itemp-1))) {
                  HFlumi_ETsum_perwedge->Fill(iw+18,tmpET);
                  if((digi.id().ieta()==-33)||(digi.id().ieta()==-34)) {
                    if(jadc>100) HFlumi_Occupancy_above_thr_r1->Fill(iw+18,1);
                    if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r1->Fill(iw+18,1);
                    if(jadc<10) HFlumi_Occupancy_below_thr_r1->Fill(iw+18,1);
                  }
                  else if((digi.id().ieta()==-35)||(digi.id().ieta()==-36)) {
                    if(jadc>100) HFlumi_Occupancy_above_thr_r2->Fill(iw+18,1);
                    if((jadc>=10)&&(jadc<=100)) HFlumi_Occupancy_between_thrs_r2->Fill(iw+18,1);
                    if(jadc<10) HFlumi_Occupancy_below_thr_r2->Fill(iw+18,1);
                  }
                  iw=99;
                }
              }
            }
          }//---endif ieta inHF-
        }//---endif TS=nr6
      } 
    }//------end loop over TS for lumi
    return;
  }
}
 // void HcalBeamMonitor::processEvent(const HBHERecHit Collection&hbheHits; ...)


void HcalBeamMonitor::beginLuminosityBlock(const edm::LuminosityBlock& lumiSeg,
                                           const edm::EventSetup& c)
  
{
  // reset histograms that get updated each luminosity section

  if (LumiInOrder(lumiSeg.luminosityBlock())==false) return;
  HcalBaseDQMonitor::beginLuminosityBlock(lumiSeg,c);

  if (lumiSeg.luminosityBlock()==lastProcessedLS_) return; // we're seeing more events from current lumi section (after some break) -- should not reset histogram
  ProblemsCurrentLB->Reset();
  HFlumi_occ_LS->Reset();
  std::stringstream title;
  title <<"HFlumi occupancy for LS # " <<currentLS;
  HFlumi_occ_LS->getTH2F()->SetTitle(title.str().c_str());
  return;
} // void HcalBeamMonitor::beginLuminosityBlock()

void HcalBeamMonitor::endLuminosityBlock(const edm::LuminosityBlock& lumiSeg,
                                         const edm::EventSetup& c)
{
  if (debug_>1) std::cout <<"<HcalBeamMonitor::endLuminosityBlock>"<<std::endl;
  if (LumiInOrder(lumiSeg.luminosityBlock())==false)
    {
      if (debug_>1)  
        std::cout <<"<HcalBeamMonitor::endLuminosityBlock>  Failed LumiInOrder test!"<<std::endl;
      return;
    }
  lastProcessedLS_=lumiSeg.luminosityBlock();
  float Nentries=HFlumi_occ_LS->getBinContent(-1,-1);
  if (debug_>3) 
    std::cout <<"Number of entries in this LB = "<<Nentries<<std::endl;

  if (Nentries<minEvents_) 
    {
      // not enough entries to determine status; fill everything with -1 and return
      HFlumi_Ring1Status_vs_LS->Fill(currentLS,-1);
      HFlumi_Ring2Status_vs_LS->Fill(currentLS,-1);
      if (Online_==false)
        return;
      // write to output file if required (Online running)
      if (lumiqualitydir_.size()==0)
        return;
      // dump out lumi quality file
      std::ofstream outStream(outfile_.str().c_str(),std::ios::app);
      outStream<<currentLS<<"\t\t-1\t\t\t-1\t\t\t-1\t\t"<<Nentries<<std::endl;
      outStream.close();
      return;
    }
  if (Nentries==0) return;


  HFlumi_total_deadcells->Fill(-1,-1,1); // counts good lumi sections in underflow bin
  HFlumi_total_hotcells->Fill(-1,-1,1);
  HFlumi_diag_deadcells->Fill(-1,-1,1); // counts good lumi sections in underflow bin
  HFlumi_diag_hotcells->Fill(-1,-1,1);

  // ADD IETA MAP
  int ietamap[8]={-36,-35,-34,-33,33,34,35,36};
  int ieta=-1, iphi = -1, depth=-1;
  int badring1=0;
  int badring2=0;
  int ndeadcells=0;
  int nhotcells=0;
  
  // Loop over cells once to count hot & dead chanels
  for (int x=1;x<=HFlumi_occ_LS->getTH2F()->GetNbinsX();++x)
    {
      for (int y=1;y<=HFlumi_occ_LS->getTH2F()->GetNbinsY();++y)
        {

          // Skip over channels that are flagged as bad
          if (x<=8)
            ieta=ietamap[x-1];
          else
            ieta=-1;
          iphi=2*y-1;
          if (abs(ieta)==33 || abs(ieta)==34)  depth=1;
          else if (abs(ieta)==35 || abs(ieta)==36) depth =2;
          else depth = -1;
          if (depth !=-1 && ieta!=1)
            {
              HcalDetId thisID(HcalForward, ieta, iphi, depth);
              if (BadCells_.find(thisID)!=BadCells_.end())
                continue;
            }
          double Ncellhits=HFlumi_occ_LS->getBinContent(x,y);
          if (Ncellhits==0)
            {
              ++ndeadcells;
              HFlumi_diag_deadcells->Fill(x-1,2*y-1,1);
            }
          // hot if present in more than 25% of events in the LS
          if (Ncellhits>hotrate_*Nentries) 
            {
              ++nhotcells;
              HFlumi_diag_hotcells->Fill(x-1,2*y-1,1);
            }
          if (Ncellhits==0 || Ncellhits>hotrate_*Nentries) // cell was either hot or dead
            {
              if (depth==1)  badring1++;
              else if (depth==2)  badring2++;
            }
        } // loop over y
    } // loop over x

  // Fill problem histogram underflow bind with number of events
  ProblemsCurrentLB->Fill(-1,-1,levt_);
  if (ndeadcells+nhotcells>=minBadCells_)
    {
      // Fill with number of error channels * events (assume bad for all events in LS)
      ProblemsCurrentLB->Fill(6,0,(ndeadcells+nhotcells)*levt_);
      for (int x=1;x<=HFlumi_occ_LS->getTH2F()->GetNbinsX();++x)
        {
          for (int y=1;y<=HFlumi_occ_LS->getTH2F()->GetNbinsY();++y)
            {
              if (x<=8)
                ieta=ietamap[x-1];
              else
                ieta=-1;
              iphi=2*y-1;
              if (abs(ieta)==33 || abs(ieta)==34)  depth=1;
              else if (abs(ieta)==35 || abs(ieta)==36) depth =2;
              else depth = -1;
              if (depth !=-1 && ieta!=1)
                {
                  // skip over channels that are flagged as bad
                  HcalDetId thisID(HcalForward, ieta, iphi, depth);
                  if (BadCells_.find(thisID)!=BadCells_.end())
                    continue;
                }
              double Ncellhits=HFlumi_occ_LS->getBinContent(x,y);
              if (Ncellhits==0)
                {
                  // One new luminosity section found with no entries for the cell in question
                  HFlumi_total_deadcells->Fill(x-1,2*y-1,1);
                } // dead cell check
              
              // hot if present in more than 25% of events in the LS
              if (Ncellhits>hotrate_*Nentries)
                {
                  HFlumi_total_hotcells->Fill(x-1,2*y-1,1);
                } // hot cell check
            } // loop over y
        } // loop over x
    } // if (ndeadcells+nhotcells>=minBadCells_)

  // Fill fraction of bad channels found in this LS
  double ring1status=0;
  double ring2status=0;
  if (ring1totalchannels_==0)
    ring1status=0;
  else
    ring1status=1-1.*badring1/ring1totalchannels_;
  HFlumi_Ring1Status_vs_LS->Fill(currentLS,ring1status);
  if (ring2totalchannels_==0)
    ring2status=0;
  else
    ring2status=1-1.*badring2/ring2totalchannels_;
  HFlumi_Ring2Status_vs_LS->Fill(currentLS,ring2status);  
  
  // Good status:  ring1 and ring2 status both > 90%
  int totalstatus=0;
  if (ring1status>0.9 && ring2status>0.9)
    totalstatus=1;
  else 
    {
      if (ring1status<=0.9)
        totalstatus-=2;
      if (ring2status<=0.9)
        totalstatus-=4;
    }

  if (lumiqualitydir_.size()==0)
    return;
  // dump out lumi quality file
  std::ofstream outStream(outfile_.str().c_str(),std::ios::app);
  outStream.precision(6);
  outStream<<currentLS<<"\t\t"<<ring1status<<"\t\t"<<ring2status<<"\t\t"<<totalstatus<<"\t\t"<<Nentries<<std::endl;
  outStream.close();
  return;
}

const float HcalBeamMonitor::area[]={0.111,0.175,0.175,0.175,0.175,0.175,0.174,0.178,0.172,0.175,0.178,0.346,0.604};
const float HcalBeamMonitor::radius[]={1300,1162,975,818,686,576,483,406,340,286,240,201,169};

void HcalBeamMonitor::SetEtaLabels(MonitorElement * h)
{
  h->getTH2F()->GetXaxis()->SetBinLabel(1,"-36S");
  h->getTH2F()->GetXaxis()->SetBinLabel(2,"-35S");
  h->getTH2F()->GetXaxis()->SetBinLabel(3,"-34L");
  h->getTH2F()->GetXaxis()->SetBinLabel(4,"-33L");
  h->getTH2F()->GetXaxis()->SetBinLabel(5,"33L");
  h->getTH2F()->GetXaxis()->SetBinLabel(6,"34L");
  h->getTH2F()->GetXaxis()->SetBinLabel(7,"35S");
  h->getTH2F()->GetXaxis()->SetBinLabel(8,"36S");
  return;
}

DEFINE_FWK_MODULE(HcalBeamMonitor);