EcalTrivialObjectAnalyzer.cc

Go to the documentation of this file.

//
// Created: 2 Mar 2006
//          Shahram Rahatlou, University of Rome & INFN
//
#include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialObjectAnalyzer.h"

#include <stdexcept>
#include <string>
#include <iostream>
#include <iomanip>
#include <map>
#include <vector>
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/EcalDetId/interface/EBDetId.h"

#include "CondFormats/EcalObjects/interface/EcalPedestals.h"
#include "CondFormats/DataRecord/interface/EcalPedestalsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalXtalGroupId.h"
#include "CondFormats/EcalObjects/interface/EcalWeightXtalGroups.h"
#include "CondFormats/DataRecord/interface/EcalWeightXtalGroupsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
#include "CondFormats/EcalObjects/interface/EcalTBWeights.h"
#include "CondFormats/DataRecord/interface/EcalTBWeightsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
#include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalIntercalibErrors.h"
#include "CondFormats/DataRecord/interface/EcalIntercalibErrorsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
#include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalTimeCalibErrors.h"
#include "CondFormats/DataRecord/interface/EcalTimeCalibErrorsRcd.h"

#include "CondFormats/EcalObjects/interface/EcalTimeOffsetConstant.h"
#include "CondFormats/DataRecord/interface/EcalTimeOffsetConstantRcd.h"

#include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
#include "CondFormats/EcalObjects/interface/EcalGainRatios.h"
#include "CondFormats/DataRecord/interface/EcalGainRatiosRcd.h"

#include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
#include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"

#include "CondFormats/EcalObjects/interface/EcalLaserAlphas.h"
#include "CondFormats/DataRecord/interface/EcalLaserAlphasRcd.h"
 
#include "CondFormats/EcalObjects/interface/EcalLaserAPDPNRatiosRef.h"
#include "CondFormats/DataRecord/interface/EcalLaserAPDPNRatiosRefRcd.h"
 
#include "CondFormats/EcalObjects/interface/EcalLaserAPDPNRatios.h"
#include "CondFormats/DataRecord/interface/EcalLaserAPDPNRatiosRcd.h"

#include "CondFormats/EcalObjects/interface/EcalClusterLocalContCorrParameters.h"
#include "CondFormats/DataRecord/interface/EcalClusterLocalContCorrParametersRcd.h"
#include "CondFormats/EcalObjects/interface/EcalClusterCrackCorrParameters.h"
#include "CondFormats/DataRecord/interface/EcalClusterCrackCorrParametersRcd.h"
#include "CondFormats/EcalObjects/interface/EcalClusterEnergyCorrectionParameters.h"
#include "CondFormats/DataRecord/interface/EcalClusterEnergyCorrectionParametersRcd.h"
#include "CondFormats/EcalObjects/interface/EcalClusterEnergyUncertaintyParameters.h"
#include "CondFormats/DataRecord/interface/EcalClusterEnergyUncertaintyParametersRcd.h"
#include "CondFormats/EcalObjects/interface/EcalClusterEnergyCorrectionObjectSpecificParameters.h"
#include "CondFormats/DataRecord/interface/EcalClusterEnergyCorrectionObjectSpecificParametersRcd.h"

#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"

#include "CondFormats/EcalObjects/interface/EcalSampleMask.h"
#include "CondFormats/DataRecord/interface/EcalSampleMaskRcd.h"

#include "CLHEP/Matrix/Matrix.h"

using namespace std;

void EcalTrivialObjectAnalyzer::analyze(const edm::Event& e, const edm::EventSetup& context)
{
    using namespace edm::eventsetup;
    // Context is not used.
    edm::LogInfo(">>> EcalTrivialObjectAnalyzer: processing run ")<<e.id().run() << " event: " << e.id().event();

    edm::ESHandle<EcalPedestals> pPeds;
    context.get<EcalPedestalsRcd>().get(pPeds);

    // ADC -> GeV Scale
    edm::ESHandle<EcalADCToGeVConstant> pAgc;
    context.get<EcalADCToGeVConstantRcd>().get(pAgc);
    const EcalADCToGeVConstant* agc = pAgc.product();
    edm::LogInfo("Global ADC->GeV scale: EB ") << std::setprecision(6) << agc->getEBValue() << " GeV/ADC count" 
          " EE " << std::setprecision(6) << agc->getEEValue() << " GeV/ADC count";

    // use a channel to fetch values from DB
    double r1 = (double)std::rand()/( double(RAND_MAX)+double(1) );
    int ieta =  int( 1 + r1*85 );
    r1 = (double)std::rand()/( double(RAND_MAX)+double(1) );
    int iphi =  int( 1 + r1*20 );

    EBDetId ebid(ieta,iphi); //eta,phi
    edm::LogInfo("EcalTrivialObjectAnalyzer: using EBDetId: ") << ebid;

    const EcalPedestals* myped=pPeds.product();
    EcalPedestals::const_iterator it = myped->find(ebid.rawId());
    if( it!=myped->end() ){
      edm::LogInfo("EcalPedestal: ")
                << "  mean_x1:  " << std::setprecision(8) << (*it).mean_x1 << " rms_x1: " << (*it).rms_x1
                << "  mean_x6:  " <<(*it).mean_x6 << " rms_x6: " << (*it).rms_x6
                << "  mean_x12: " <<(*it).mean_x12 << " rms_x12: " << (*it).rms_x12
                ;
    } else {
      edm::LogInfo("No pedestal found for this xtal! something wrong with EcalPedestals in your DB? ")
               ;
    }

    // fetch map of groups of xtals
    edm::ESHandle<EcalWeightXtalGroups> pGrp;
    context.get<EcalWeightXtalGroupsRcd>().get(pGrp);
    const EcalWeightXtalGroups* grp = pGrp.product();

    EcalXtalGroupsMap::const_iterator git = grp->getMap().find( ebid.rawId() );
    EcalXtalGroupId gid;
    if( git != grp->getMap().end() ) {
      edm::LogInfo("XtalGroupId.id() = ") << std::setprecision(3) << (*git).id() ;
      gid = (*git);
    } else {
      edm::LogInfo("No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?")
                ;
   }

    // Gain Ratios
    edm::ESHandle<EcalGainRatios> pRatio;
    context.get<EcalGainRatiosRcd>().get(pRatio);
    const EcalGainRatios* gr = pRatio.product();

    EcalGainRatioMap::const_iterator grit=gr->getMap().find(ebid.rawId());
    EcalMGPAGainRatio mgpa;
    if( grit!=gr->getMap().end() ){
      mgpa = (*grit);

      edm::LogInfo("EcalMGPAGainRatio: ")
                << "gain 12/6 :  " << std::setprecision(4) << mgpa.gain12Over6() << " gain 6/1: " << mgpa.gain6Over1()
                ;
    } else {
      edm::LogInfo("No MGPA Gain Ratio found for this xtal! something wrong with EcalGainRatios in your DB? ")
               ;
    }

    // Intercalib constants
    edm::ESHandle<EcalIntercalibConstants> pIcal;
    context.get<EcalIntercalibConstantsRcd>().get(pIcal);
    const EcalIntercalibConstants* ical = pIcal.product();

    EcalIntercalibConstantMap::const_iterator icalit=ical->getMap().find(ebid.rawId());
    EcalIntercalibConstant icalconst;
    if( icalit!=ical->getMap().end() ){
      icalconst = (*icalit);

      edm::LogInfo("EcalIntercalibConstant: ")
                <<std::setprecision(6)
                << icalconst
                ;
    } else {
      edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalIntercalibConstants in your DB? ")
               ;
    }

    // Intercalib errors
    edm::ESHandle<EcalIntercalibErrors> pIcalErr;
    context.get<EcalIntercalibErrorsRcd>().get(pIcalErr);
    const EcalIntercalibErrors* icalErr = pIcalErr.product();

    EcalIntercalibErrorMap::const_iterator icalitErr=icalErr->getMap().find(ebid.rawId());
    EcalIntercalibError icalconstErr;
    if( icalitErr!=icalErr->getMap().end() ){
      icalconstErr = (*icalitErr);

      edm::LogInfo("EcalIntercalibError: ")
                <<std::setprecision(6)
                << icalconstErr
                ;
    } else {
      edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalIntercalibErrors in your DB? ")
               ;
    }



    { // quick and dirty for cut and paste ;) it is a test program...
            // TimeCalib constants
            edm::ESHandle<EcalTimeCalibConstants> pIcal;
            context.get<EcalTimeCalibConstantsRcd>().get(pIcal);
            const EcalTimeCalibConstants* ical = pIcal.product();

            EcalTimeCalibConstantMap::const_iterator icalit=ical->getMap().find(ebid.rawId());
            EcalTimeCalibConstant icalconst;
            if( icalit!=ical->getMap().end() ){
                    icalconst = (*icalit);

                    edm::LogInfo("EcalTimeCalibConstant: ")
                            <<std::setprecision(6)
                            << icalconst
                            ;
            } else {
          edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalTimeCalibConstants in your DB? ")
                            ;
            }

            // TimeCalib errors
            edm::ESHandle<EcalTimeCalibErrors> pIcalErr;
            context.get<EcalTimeCalibErrorsRcd>().get(pIcalErr);
            const EcalTimeCalibErrors* icalErr = pIcalErr.product();

            EcalTimeCalibErrorMap::const_iterator icalitErr=icalErr->getMap().find(ebid.rawId());
            EcalTimeCalibError icalconstErr;
            if( icalitErr!=icalErr->getMap().end() ){
                    icalconstErr = (*icalitErr);

                    edm::LogInfo("EcalTimeCalibError: ")
                            <<std::setprecision(6)
                            << icalconstErr
                            ;
            } else {
          edm::LogInfo("No intercalib const found for this xtal! something wrong with EcalTimeCalibErrors in your DB? ")
                            ;
            }
    }

   // fetch Time Offset
   //std::cout <<"Fetching TimeOffsetConstant from DB " << std::endl;

   // Time Offset constants
   edm::ESHandle<EcalTimeOffsetConstant> pTOff;
   context.get<EcalTimeOffsetConstantRcd>().get(pTOff);
   const EcalTimeOffsetConstant* TOff = pTOff.product();

   edm::LogInfo("EcalTimeOffsetConstant: ")
         << " EB " << TOff->getEBValue()
         << " EE " << TOff->getEEValue()
         ;

   // fetch TB weights
   edm::LogInfo("Fetching EcalTBWeights from DB ");
   edm::ESHandle<EcalTBWeights> pWgts;
   context.get<EcalTBWeightsRcd>().get(pWgts);
   const EcalTBWeights* wgts = pWgts.product();
   edm::LogInfo("EcalTBWeightMap.size(): ") << std::setprecision(3) << wgts->getMap().size();


   // look up the correct weights for this  xtal
   //EcalXtalGroupId gid( (*git) );
   EcalTBWeights::EcalTDCId tdcid(1);

   edm::LogInfo("Lookup EcalWeightSet for groupid: ") << std::setprecision(3) 
             << gid.id() << " and TDC id " << tdcid ;
   EcalTBWeights::EcalTBWeightMap::const_iterator wit = wgts->getMap().find( std::make_pair(gid,tdcid) );
   EcalWeightSet  wset;
   if( wit != wgts->getMap().end() ) {
      wset = wit->second;
      edm::LogInfo("check size of data members in EcalWeightSet") ;
      //wit->second.print(std::cout);


      //typedef std::vector< std::vector<EcalWeight> > EcalWeightMatrix;
      const EcalWeightSet::EcalWeightMatrix& mat1 = wit->second.getWeightsBeforeGainSwitch();
      const EcalWeightSet::EcalWeightMatrix& mat2 = wit->second.getWeightsAfterGainSwitch();

//       std::cout << "WeightsBeforeGainSwitch.size: " << mat1.size()
//                 << ", WeightsAfterGainSwitch.size: " << mat2.size() << std::endl;


      CLHEP::HepMatrix clmat1(3,10,0);
      CLHEP::HepMatrix clmat2(3,10,0);
      for(int irow=0; irow<3; irow++) {
       for(int icol=0; icol<10; icol++) {
         clmat1[irow][icol] = mat1(irow,icol);
         clmat2[irow][icol] = mat2(irow,icol);
       }
     }
      edm::LogInfo("weight matrix before gain switch:") << clmat1;
      edm::LogInfo("weight matrix after gain switch:") << clmat2;

   } else {
     edm::LogInfo("No weights found for EcalGroupId: ") << gid.id() << " and  EcalTDCId: " << tdcid ;
  }

   // cluster functions/corrections
   edm::ESHandle<EcalClusterLocalContCorrParameters> pLocalCont;
   context.get<EcalClusterLocalContCorrParametersRcd>().get(pLocalCont);
   const EcalClusterLocalContCorrParameters* paramLocalCont = pLocalCont.product();
   edm::LogInfo("LocalContCorrParameters:");
   for ( EcalFunctionParameters::const_iterator it = paramLocalCont->params().begin(); it != paramLocalCont->params().end(); ++it ) {
     edm::LogInfo(" ") << *it;
   }
   // std::cout << "\n";
   edm::ESHandle<EcalClusterCrackCorrParameters> pCrack;
   context.get<EcalClusterCrackCorrParametersRcd>().get(pCrack);
   const EcalClusterCrackCorrParameters* paramCrack = pCrack.product();
   edm::LogInfo("CrackCorrParameters:");
   for ( EcalFunctionParameters::const_iterator it = paramCrack->params().begin(); it != paramCrack->params().end(); ++it ) {
     edm::LogInfo(" ") << *it;
   }
   // std::cout << "\n";
   edm::ESHandle<EcalClusterEnergyCorrectionParameters> pEnergyCorrection;
   context.get<EcalClusterEnergyCorrectionParametersRcd>().get(pEnergyCorrection);
   const EcalClusterEnergyCorrectionParameters* paramEnergyCorrection = pEnergyCorrection.product();
   edm::LogInfo("EnergyCorrectionParameters:");
   for ( EcalFunctionParameters::const_iterator it = paramEnergyCorrection->params().begin(); it != paramEnergyCorrection->params().end(); ++it ) {
     edm::LogInfo(" ") << *it;
   }
   // std::cout << "\n";
   edm::ESHandle<EcalClusterEnergyUncertaintyParameters> pEnergyUncertainty;
   context.get<EcalClusterEnergyUncertaintyParametersRcd>().get(pEnergyUncertainty);
   const EcalClusterEnergyUncertaintyParameters* paramEnergyUncertainty = pEnergyUncertainty.product();
   edm::LogInfo("EnergyCorrectionParameters:");
   for ( EcalFunctionParameters::const_iterator it = paramEnergyUncertainty->params().begin(); it != paramEnergyUncertainty->params().end(); ++it ) {
     edm::LogInfo(" ") << *it;
   }
   // std::cout << "\n";
   edm::ESHandle<EcalClusterEnergyCorrectionObjectSpecificParameters> pEnergyCorrectionObjectSpecific;
   context.get<EcalClusterEnergyCorrectionObjectSpecificParametersRcd>().get(pEnergyCorrectionObjectSpecific);
   const EcalClusterEnergyCorrectionObjectSpecificParameters* paramEnergyCorrectionObjectSpecific = pEnergyCorrectionObjectSpecific.product();
   edm::LogInfo("EnergyCorrectionObjectSpecificParameters:");
   for ( EcalFunctionParameters::const_iterator it = paramEnergyCorrectionObjectSpecific->params().begin(); it != paramEnergyCorrectionObjectSpecific->params().end(); ++it ) {
     edm::LogInfo(" ") << *it;
   }
   // std::cout << "\n";


   // laser correction
   
   // laser alphas
   edm::ESHandle<EcalLaserAlphas> pAlpha; 
   context.get<EcalLaserAlphasRcd>().get(pAlpha);
   const EcalLaserAlphas* lalpha = pAlpha.product();
   
   EcalLaserAlphaMap::const_iterator lalphait;
   lalphait = lalpha->getMap().find(ebid.rawId());
   if( lalphait!=lalpha->getMap().end() ){
     edm::LogInfo("EcalLaserAlpha: ")
           <<std::setprecision(6)
           << (*lalphait)
           ;
   } else {
     edm::LogInfo("No laser alpha found for this xtal! something wrong with EcalLaserAlphas in your DB? ")
               ;
   }

   // laser apdpnref
   edm::ESHandle<EcalLaserAPDPNRatiosRef> pAPDPNRatiosRef;          
   context.get<EcalLaserAPDPNRatiosRefRcd>().get(pAPDPNRatiosRef);  
   const EcalLaserAPDPNRatiosRef* lref = pAPDPNRatiosRef.product(); 
   
   EcalLaserAPDPNRatiosRef::const_iterator lrefit;
   lrefit = lref->getMap().find(ebid.rawId());
   if( lrefit!=lref->getMap().end() ){
     edm::LogInfo("EcalLaserAPDPNRatiosRef: ")
           <<std::setprecision(6)
           << (*lrefit)
           ;
   } else {
     edm::LogInfo("No laser apd/pn ref found for this xtal! something wrong with EcalLaserAPDPNRatiosRef in your DB? ")
           ;
   }

   // laser apdpn ratios 
   edm::ESHandle<EcalLaserAPDPNRatios> pAPDPNRatios;          
   context.get<EcalLaserAPDPNRatiosRcd>().get(pAPDPNRatios);  
   const EcalLaserAPDPNRatios* lratio = pAPDPNRatios.product();

   EcalLaserAPDPNRatios::EcalLaserAPDPNRatiosMap::const_iterator lratioit; 
   lratioit = lratio->getLaserMap().find(ebid.rawId());
   EcalLaserAPDPNRatios::EcalLaserAPDPNpair lratioconst;
   
   if( lratioit!=lratio->getLaserMap().end() ){
     lratioconst = (*lratioit);
     
     edm::LogInfo("EcalLaserAPDPNRatios: ")
       //          <<e.id().run() << " " << e.id().event() << " " 
               << std::setprecision(6)
               << lratioconst.p1 << " " << lratioconst.p2
               ;
   } else {
     edm::LogInfo("No laser apd/pn ratio found for this xtal! something wrong with EcalLaserAPDPNRatios in your DB? ")
           ;
   }

   // laser timestamps
   EcalLaserAPDPNRatios::EcalLaserTimeStamp ltimestamp;
   // EcalLaserAPDPNRatios::EcalLaserTimeStampMap::const_iterator ltimeit;
   for (int i=1; i<=92; i++) {     
     ltimestamp = lratio->getTimeMap()[i];
     edm::LogInfo("i = ") << std::setprecision(6) << i
             << ltimestamp.t1.value() << " " << ltimestamp.t2.value() << " : " ;
   }
   edm::LogInfo("Tests finished.");

   // channel status
   edm::ESHandle<EcalChannelStatus> pChannelStatus;
   context.get<EcalChannelStatusRcd>().get(pChannelStatus);
   const EcalChannelStatus *ch_status = pChannelStatus.product();

   EcalChannelStatusMap::const_iterator chit;
   chit = ch_status->getMap().find(ebid.rawId());
   if( chit != ch_status->getMap().end() ){
           EcalChannelStatusCode ch_code = (*chit);
           edm::LogInfo("EcalChannelStatus: ")
                   <<std::setprecision(6)
                   << ch_code.getStatusCode()
                   ;
   } else {
     edm::LogInfo("No channel status found for this xtal! something wrong with EcalChannelStatus in your DB? ")
                   ;
   }


   // laser transparency correction

   // Mask to ignore sample
   edm::ESHandle<EcalSampleMask> pSMask;
   context.get<EcalSampleMaskRcd>().get(pSMask);
   const EcalSampleMask* smask = pSMask.product();
   edm::LogInfo("Sample Mask EB ") << std::hex << smask->getEcalSampleMaskRecordEB() 
         << " EE " << std::hex << smask->getEcalSampleMaskRecordEE();
   

/*
    std::cout << "make CLHEP matrices from vector<vector<Ecalweight>>" << std::endl;
    CLHEP::HepMatrix clmat1(3,8,0);
    CLHEP::HepMatrix clmat2(3,8,0);
    for(int irow=0; irow<3; irow++) {
     for(int icol=0; icol<8; icol++) {
       clmat1[irow][icol] = (mat1[irow])[icol]();
       clmat2[irow][icol] = (mat2[irow])[icol]();
     }
   }
   std::cout << clmat1 << std::endl;
   std::cout << clmat2 << std::endl;
*/

} //end of ::Analyze()
  //DEFINE_FWK_MODULE(EcalTrivialObjectAnalyzer);