DQMHcalPhiSymAlCaReco.cc

Go to the documentation of this file.

/*
 * \file DQMHcalPhiSymAlCaReco.cc
 *
 * \author Olga Kodolova
 * 
 *
 *
 * Description: Monitoring of Phi Symmetry Calibration Stream  
*/

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

// DQM include files

#include "DQMServices/Core/interface/MonitorElement.h"

// work on collections

#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "DataFormats/FEDRawData/interface/FEDRawData.h"
#include "DataFormats/FEDRawData/interface/FEDNumbering.h"
#include "DataFormats/FEDRawData/interface/FEDHeader.h"
#include "EventFilter/HcalRawToDigi/interface/HcalHTRData.h"
#include "EventFilter/HcalRawToDigi/interface/HcalDCCHeader.h"

#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMOffline/CalibCalo/src/DQMHcalPhiSymAlCaReco.h"

using namespace std;
using namespace edm;

// ******************************************
// constructors
// *****************************************

DQMHcalPhiSymAlCaReco::DQMHcalPhiSymAlCaReco( const edm::ParameterSet& ps ) :
eventCounter_(0)
{
  //
  // Input from configurator file 
  //
  folderName_   = ps.getUntrackedParameter<string>("FolderName","ALCAStreamHcalPhiSym");
  
  hbherecoMB    = consumes<HBHERecHitCollection>(ps.getParameter<edm::InputTag>("hbheInputMB"));
  horecoMB      = ps.getParameter<edm::InputTag>("hoInputMB");
  hfrecoMB      = consumes<HFRecHitCollection>(ps.getParameter<edm::InputTag>("hfInputMB"));
  
  hbherecoNoise = consumes<HBHERecHitCollection>(ps.getParameter<edm::InputTag>("hbheInputNoise"));
  horecoNoise   = ps.getParameter<edm::InputTag>("hoInputNoise");
  hfrecoNoise   = consumes<HFRecHitCollection>(ps.getParameter<edm::InputTag>("hfInputNoise"));

  rawInLabel_   = consumes<FEDRawDataCollection>(ps.getParameter<edm::InputTag>("rawInputLabel"));

  period_       = ps.getParameter<unsigned int>("period") ;

  saveToFile_   = ps.getUntrackedParameter<bool>("SaveToFile",false);
  fileName_     = ps.getUntrackedParameter<string>("FileName","MonitorAlCaHcalPhiSym.root");

  // histogram parameters

  // Distribution of rechits in iPhi, iEta 
  hiDistr_y_nbin_ = ps.getUntrackedParameter<int>("hiDistr_y_nbin",72);
  hiDistr_y_min_  =  ps.getUntrackedParameter<double>("hiDistr_y_min",0.5);
  hiDistr_y_max_  =  ps.getUntrackedParameter<double>("hiDistr_y_max",72.5);
  hiDistr_x_nbin_ = ps.getUntrackedParameter<int>("hiDistr_x_nbin",41);
  hiDistr_x_min_  =  ps.getUntrackedParameter<double>("hiDistr_x_min",0.5);
  hiDistr_x_max_  =  ps.getUntrackedParameter<double>("hiDistr_x_max",41.5);
  // Check for NZS
  hiDistr_r_nbin_ = ps.getUntrackedParameter<int>("hiDistr_r_nbin",100);
  ihbhe_size_     = ps.getUntrackedParameter<double>("ihbhe_size_",5184.);
  ihf_size_       = ps.getUntrackedParameter<double>("ihf_size_",1728.);
    
}

DQMHcalPhiSymAlCaReco::~DQMHcalPhiSymAlCaReco()
{}

//--------------------------------------------------------
void DQMHcalPhiSymAlCaReco::bookHistograms(DQMStore::IBooker & ibooker, edm::Run const & irun, edm::EventSetup const & isetup) {
 
  // create and cd into new folder
  ibooker.setCurrentFolder(folderName_);

  eventCounter_ = 0;

  hFEDsize = ibooker.book1D("hFEDsize","HCAL FED size (kB)",200,-0.5,20.5);
  hFEDsize->setAxisTitle("kB",1);

  hHcalIsZS = ibooker.book1D("hHcalIsZS", "Hcal Is ZS", 4, -1.5, 2.5);
  hHcalIsZS->setBinLabel(2, "NZS");
  hHcalIsZS->setBinLabel(3, "ZS");

  char hname[50];
  sprintf(hname, "L1 Event Number %% %i", period_);
  hL1Id = ibooker.book1D("hL1Id", hname,4200,-99.5,4099.5);
  hL1Id->setAxisTitle(hname);
  

  // book some histograms 1D
  double xmin = 0.1;
  double xmax = 1.1;
  hiDistrHBHEsize1D_ =
    ibooker.book1D("DistrHBHEsize","Size of HBHE Collection",
                  hiDistr_r_nbin_,
                  xmin,
                  xmax
                );
  hiDistrHFsize1D_ =
    ibooker.book1D("DistrHFsize","Size of HF Collection",
                  hiDistr_r_nbin_,
                  xmin,
                  xmax
                );


  // First moment
  hiDistrMBPl2D_ = 
    ibooker.book2D("MBdepthPl1", "iphi- +ieta signal distribution at depth1",
         hiDistr_x_nbin_, 
         hiDistr_x_min_,
         hiDistr_x_max_,
         hiDistr_y_nbin_, 
         hiDistr_y_min_,
         hiDistr_y_max_
         );

  hiDistrMBPl2D_->setAxisTitle("i#phi ", 2);
  hiDistrMBPl2D_->setAxisTitle("i#eta ", 1);


  hiDistrNoisePl2D_ = 
    ibooker.book2D("NoisedepthPl1", "iphi-ieta noise distribution at depth1",
         hiDistr_x_nbin_+1, 
         hiDistr_x_min_-1.,
         hiDistr_x_max_,
         hiDistr_y_nbin_+1, 
         hiDistr_y_min_-1.,
         hiDistr_y_max_
         );

  hiDistrNoisePl2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoisePl2D_->setAxisTitle("i#eta ", 1);
// Second moment
  hiDistrMB2Pl2D_ =
    ibooker.book2D("MB2depthPl1", "iphi- +ieta signal distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrMB2Pl2D_->setAxisTitle("i#phi ", 2);
  hiDistrMB2Pl2D_->setAxisTitle("i#eta ", 1);


  hiDistrNoise2Pl2D_ =
    ibooker.book2D("Noise2depthPl1", "iphi-ieta noise distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrNoise2Pl2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoise2Pl2D_->setAxisTitle("i#eta ", 1);

// Variance
  hiDistrVarMBPl2D_ =
    ibooker.book2D("VarMBdepthPl1", "iphi- +ieta signal distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrVarMBPl2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarMBPl2D_->setAxisTitle("i#eta ", 1);


  hiDistrVarNoisePl2D_ =
    ibooker.book2D("VarNoisedepthPl1", "iphi-ieta noise distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrVarNoisePl2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarNoisePl2D_->setAxisTitle("i#eta ", 1);

//==================================================================================
// First moment
  hiDistrMBMin2D_ = 
    ibooker.book2D("MBdepthMin1", "iphi- +ieta signal distribution at depth1",
         hiDistr_x_nbin_, 
         hiDistr_x_min_,
         hiDistr_x_max_,
         hiDistr_y_nbin_, 
         hiDistr_y_min_,
         hiDistr_y_max_
         );

  hiDistrMBMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrMBMin2D_->setAxisTitle("i#eta ", 1);


  hiDistrNoiseMin2D_ = 
    ibooker.book2D("NoisedepthMin1", "iphi-ieta noise distribution at depth1",
         hiDistr_x_nbin_, 
         hiDistr_x_min_,
         hiDistr_x_max_,
         hiDistr_y_nbin_, 
         hiDistr_y_min_,
         hiDistr_y_max_
         );

  hiDistrNoiseMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoiseMin2D_->setAxisTitle("i#eta ", 1);
// Second moment
  hiDistrMB2Min2D_ =
    ibooker.book2D("MB2depthMin1", "iphi- +ieta signal distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrMB2Min2D_->setAxisTitle("i#phi ", 2);
  hiDistrMB2Min2D_->setAxisTitle("i#eta ", 1);


  hiDistrNoise2Min2D_ =
    ibooker.book2D("Noise2depthMin1", "iphi-ieta noise distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrNoise2Min2D_->setAxisTitle("i#phi ", 2);
  hiDistrNoise2Min2D_->setAxisTitle("i#eta ", 1);

// Variance
  hiDistrVarMBMin2D_ =
    ibooker.book2D("VarMBdepthMin1", "iphi- +ieta signal distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrVarMBMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarMBMin2D_->setAxisTitle("i#eta ", 1);


  hiDistrVarNoiseMin2D_ =
    ibooker.book2D("VarNoisedepthMin1", "iphi-ieta noise distribution at depth1",
                 hiDistr_x_nbin_,
                 hiDistr_x_min_,
                 hiDistr_x_max_,
                 hiDistr_y_nbin_,
                 hiDistr_y_min_,
                 hiDistr_y_max_
                 );

  hiDistrVarNoiseMin2D_->setAxisTitle("i#phi ", 2);
  hiDistrVarNoiseMin2D_->setAxisTitle("i#eta ", 1);


}

//--------------------------------------------------------
//void DQMHcalPhiSymAlCaReco::beginRun(const edm::Run& r, const EventSetup& context) {
//}

//--------------------------------------------------------

//-------------------------------------------------------------

void DQMHcalPhiSymAlCaReco::analyze(const Event& iEvent, 
                   const EventSetup& iSetup ){  
 

  eventCounter_++;
 
  edm::Handle<FEDRawDataCollection> rawIn;
  iEvent.getByToken(rawInLabel_,rawIn);

  if(!rawIn.isValid()){
     LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
      return ;
  }

  //get HCAL FEDs:
  std::vector<int> selFEDs;
  for (int i=FEDNumbering::MINHCALFEDID; i<=FEDNumbering::MAXHCALFEDID; i++)
    {
      selFEDs.push_back(i);
    }

  //  std::cout<<" Size of FED "<<selFEDs.size()<<std::endl;

  const FEDRawDataCollection *rdc=rawIn.product(); 

  bool hcalIsZS = false;
  int lvl1ID = 0;
  bool lvl1IDFound = false;
  for (unsigned int k=0; k<selFEDs.size(); k++) {
    const FEDRawData & fedData = rdc->FEDData(selFEDs[k]);
    // std::cout<<fedData.size()*std::pow(1024.,-1)<<std::endl;
    hFEDsize->Fill(fedData.size()*std::pow(1024.,-1),1);
    
    // get HCAL DCC Header for each FEDRawData
    const HcalDCCHeader* dccHeader=(const HcalDCCHeader*)(fedData.data());
    if (dccHeader) {
      // walk through the HTR data...
      HcalHTRData htr;

      int nspigot =0; 
      for (int spigot=0; spigot<HcalDCCHeader::SPIGOT_COUNT; spigot++) {    
       nspigot++;

        if (!dccHeader->getSpigotPresent(spigot)) continue;

        // Load the given decoder with the pointer and length from this spigot.
        dccHeader->getSpigotData(spigot,htr, fedData.size()); 

        if(k != 20 && nspigot !=14 ) {      
          if ( !htr.isUnsuppressed() ) { hcalIsZS = true; }
        }
      } 
    }
 
    // try to get the lvl1ID from the HCAL fed
    if (!lvl1IDFound && (fedData.size() > 0)) { 
      // get FED Header for FEDRawData
      FEDHeader fedHeader(fedData.data());
      lvl1ID = fedHeader.lvl1ID();
      lvl1IDFound = true;
    }
  } // loop over HcalFEDs
  
  hHcalIsZS->Fill( hcalIsZS );
  hL1Id->Fill( lvl1ID%period_ );
 
  edm::Handle<HBHERecHitCollection> hbheNS;
  iEvent.getByToken(hbherecoNoise, hbheNS);
  
  if(!hbheNS.isValid()){
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return ;
  }
  
  edm::Handle<HBHERecHitCollection> hbheMB;
  iEvent.getByToken(hbherecoMB, hbheMB);
  
  if(!hbheMB.isValid()){
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return ;
  }
  
  edm::Handle<HFRecHitCollection> hfNS;
  iEvent.getByToken(hfrecoNoise, hfNS);
  
  if(!hfNS.isValid()){
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return ;
  }
  
  edm::Handle<HFRecHitCollection> hfMB;
  iEvent.getByToken(hfrecoMB, hfMB);
  
  if(!hfMB.isValid()){
    LogDebug("") << "HcalCalibAlgos: Error! can't get hbhe product!" << std::endl;
    return ;
  }
  
  const HBHERecHitCollection HithbheNS = *(hbheNS.product());
  
  hiDistrHBHEsize1D_->Fill(HithbheNS.size()/ihbhe_size_);
  
     
  for(HBHERecHitCollection::const_iterator hbheItr=HithbheNS.begin(); hbheItr!=HithbheNS.end(); hbheItr++)
        {
             DetId id = (*hbheItr).detid(); 
             HcalDetId hid=HcalDetId(id);
                 
         if(hid.depth() == 1) {
                 if( hid.ieta() > 0 ) {
         hiDistrNoisePl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy());
                 hiDistrNoise2Pl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy()*hbheItr->energy());
                 } else {
         hiDistrNoiseMin2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy());
                 hiDistrNoise2Min2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy()*hbheItr->energy());
         }
         }
        }

  const HBHERecHitCollection HithbheMB = *(hbheMB.product());

  for(HBHERecHitCollection::const_iterator hbheItr=HithbheMB.begin(); hbheItr!=HithbheMB.end(); hbheItr++)
        {
             DetId id = (*hbheItr).detid(); 
             HcalDetId hid=HcalDetId(id);
                 
         if(hid.depth() == 1) {
                 if( hid.ieta() > 0 ) {
         hiDistrMBPl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy());
                 hiDistrMB2Pl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy()*hbheItr->energy());
                 } else {
         hiDistrMBMin2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy());
                 hiDistrMB2Min2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy()*hbheItr->energy());
         }
         }

        }

  const HFRecHitCollection HithfNS = *(hfNS.product());

  hiDistrHFsize1D_->Fill(HithfNS.size()/ihf_size_);
  
  for(HFRecHitCollection::const_iterator hbheItr=HithfNS.begin(); hbheItr!=HithfNS.end(); hbheItr++)
        {
    
             DetId id = (*hbheItr).detid(); 
             HcalDetId hid=HcalDetId(id);
                 
         if(hid.depth() == 1) {
                 if( hid.ieta() > 0 ) {
         hiDistrNoisePl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy());
                 hiDistrNoise2Pl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy()*hbheItr->energy());
                 } else {
         hiDistrNoiseMin2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy());
                 hiDistrNoise2Min2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy()*hbheItr->energy());
         }
         }
    
        }

  const HFRecHitCollection HithfMB = *(hfMB.product());
 
  for(HFRecHitCollection::const_iterator hbheItr=HithfMB.begin(); hbheItr!=HithfMB.end(); hbheItr++)
        {
             DetId id = (*hbheItr).detid(); 
             HcalDetId hid=HcalDetId(id);
                 
         if(hid.depth() == 1) {
                 if( hid.ieta() > 0 ) {
         hiDistrMBPl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy());
                 hiDistrMB2Pl2D_->Fill(hid.ieta(),hid.iphi(),hbheItr->energy()*hbheItr->energy());
                 } else {
         hiDistrMBMin2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy());
                 hiDistrMB2Min2D_->Fill(fabs(hid.ieta()),hid.iphi(),hbheItr->energy()*hbheItr->energy());
         }
         }
        }   
    
    
} //analyze




//--------------------------------------------------------
//--------------------------------------------------------
void DQMHcalPhiSymAlCaReco::endRun(const Run& r, const EventSetup& context){
// Keep Variances
  if(eventCounter_ > 0) {
  for(int k=0; k<=hiDistr_x_nbin_;k++)
  {
    for(int j=0; j<=hiDistr_y_nbin_;j++)
    {
// First moment
       float cc1=hiDistrMBPl2D_->getBinContent(k,j);
       cc1 = cc1 * 1./eventCounter_;
       float cc2=hiDistrNoisePl2D_->getBinContent(k,j);
       cc2 = cc2 * 1./eventCounter_;
       float cc3=hiDistrMBMin2D_->getBinContent(k,j);
       cc3 = cc3 * 1./eventCounter_;
       float cc4=hiDistrNoiseMin2D_->getBinContent(k,j);
       cc4 = cc4 * 1./eventCounter_;
// Second moment
       float cc11=hiDistrMB2Pl2D_->getBinContent(k,j);
       cc11 = cc11 * 1./eventCounter_;
       hiDistrVarMBPl2D_->setBinContent(k,j,cc11-cc1*cc1);
       float cc22=hiDistrNoise2Pl2D_->getBinContent(k,j);
       cc22 = cc22 * 1./eventCounter_;
       hiDistrVarNoisePl2D_->setBinContent(k,j,cc22-cc2*cc2);
       float cc33=hiDistrMB2Min2D_->getBinContent(k,j);
       cc33 = cc33 * 1./eventCounter_;
       hiDistrVarMBMin2D_->setBinContent(k,j,cc33-cc3*cc3);
       float cc44=hiDistrNoise2Min2D_->getBinContent(k,j);
       cc44 = cc44 * 1./eventCounter_;
       hiDistrVarNoiseMin2D_->setBinContent(k,j,cc44-cc4*cc4);
    }     
  }
  }
}