EcalBarrelDigisValidation.cc

Go to the documentation of this file.

/*
 * \file EcalBarrelDigisValidation.cc
 *
 * \author F. Cossutti
 *
*/

#include <Validation/EcalDigis/interface/EcalBarrelDigisValidation.h>
#include "CalibCalorimetry/EcalTrivialCondModules/interface/EcalTrivialConditionRetriever.h"
#include "DQMServices/Core/interface/DQMStore.h"

using namespace cms;
using namespace edm;
using namespace std;

EcalBarrelDigisValidation::EcalBarrelDigisValidation(const ParameterSet& ps):
  EBdigiCollection_(consumes<EBDigiCollection>(ps.getParameter<edm::InputTag>("EBdigiCollection")))
{
  
  // verbosity switch
  verbose_ = ps.getUntrackedParameter<bool>("verbose", false);
    
                                                                                                                                          
  // get hold of back-end interface
                                                                                                                                          

  gainConv_[1] = 1.;
  gainConv_[2] = 2.;
  gainConv_[3] = 12.;
  gainConv_[0] = 12.;   // saturated channels
  barrelADCtoGeV_ = 0.035;
  endcapADCtoGeV_ = 0.06;
 
  meEBDigiOccupancy_ = nullptr;

  meEBDigiMultiplicity_ = nullptr;

  meEBDigiADCGlobal_ = nullptr;

  for (int i = 0; i < 10 ; i++ ) {
    meEBDigiADCAnalog_[i] = nullptr;
    meEBDigiADCgS_[i]  = nullptr;
    meEBDigiADCg1_[i]  = nullptr;
    meEBDigiADCg6_[i]  = nullptr;
    meEBDigiADCg12_[i] = nullptr;
    meEBDigiGain_[i] = nullptr;
  }

  meEBPedestal_ = nullptr;
                                 
  meEBMaximumgt100ADC_ = nullptr; 
                                 
  meEBMaximumgt10ADC_ = nullptr; 

  meEBnADCafterSwitch_ = nullptr;
 
}

EcalBarrelDigisValidation::~EcalBarrelDigisValidation(){
 
}

void EcalBarrelDigisValidation::bookHistograms(DQMStore::IBooker &ibooker, edm::Run const&, edm::EventSetup const&){

    Char_t histo[200];

    ibooker.setCurrentFolder("EcalDigisV/EcalDigiTask");
  
    sprintf (histo, "EcalDigiTask Barrel occupancy" ) ;
    meEBDigiOccupancy_ = ibooker.book2D(histo, histo, 360, 0., 360., 170, -85., 85.);

    sprintf (histo, "EcalDigiTask Barrel digis multiplicity" ) ;
    meEBDigiMultiplicity_ = ibooker.book1D(histo, histo, 612, 0., 61200);
  
    sprintf (histo, "EcalDigiTask Barrel global pulse shape" ) ;
    meEBDigiADCGlobal_ = ibooker.bookProfile(histo, histo, 10, 0, 10, 10000, 0., 1000.) ;
    
    for (int i = 0; i < 10 ; i++ ) {

      sprintf (histo, "EcalDigiTask Barrel analog pulse %02d", i+1) ;
      meEBDigiADCAnalog_[i] = ibooker.book1D(histo, histo, 4000, 0., 400.);

      sprintf (histo, "EcalDigiTask Barrel ADC pulse %02d Gain 0 - Saturated", i+1) ;
      meEBDigiADCgS_[i] = ibooker.book1D(histo, histo, 4096, -0.5, 4095.5);

      sprintf (histo, "EcalDigiTask Barrel ADC pulse %02d Gain 1", i+1) ;
      meEBDigiADCg1_[i] = ibooker.book1D(histo, histo, 4096, -0.5, 4095.5);

      sprintf (histo, "EcalDigiTask Barrel ADC pulse %02d Gain 6", i+1) ;
      meEBDigiADCg6_[i] = ibooker.book1D(histo, histo, 4096, -0.5, 4095.5);

      sprintf (histo, "EcalDigiTask Barrel ADC pulse %02d Gain 12", i+1) ;
      meEBDigiADCg12_[i] = ibooker.book1D(histo, histo, 4096, -0.5, 4095.5);

      sprintf (histo, "EcalDigiTask Barrel gain pulse %02d", i+1) ;
      meEBDigiGain_[i] = ibooker.book1D(histo, histo, 4, 0, 4);

    }
    
    sprintf (histo, "EcalDigiTask Barrel pedestal for pre-sample" ) ;
    meEBPedestal_ = ibooker.book1D(histo, histo, 4096, -0.5, 4095.5) ;

    sprintf (histo, "EcalDigiTask Barrel maximum position gt 100 ADC" ) ;
    meEBMaximumgt100ADC_ = ibooker.book1D(histo, histo, 10, 0., 10.) ;

    sprintf (histo, "EcalDigiTask Barrel maximum position gt 10 ADC" ) ;
    meEBMaximumgt10ADC_ = ibooker.book1D(histo, histo, 10, 0., 10.) ;

    sprintf (histo, "EcalDigiTask Barrel ADC counts after gain switch" ) ;
    meEBnADCafterSwitch_ = ibooker.book1D(histo, histo, 10, 0., 10.) ;

}

void EcalBarrelDigisValidation::dqmBeginRun(edm::Run const&, edm::EventSetup const& c){

  checkCalibrations(c);

}


void EcalBarrelDigisValidation::analyze(Event const & e, EventSetup const & c){

  //LogInfo("EventInfo") << " Run = " << e.id().run() << " Event = " << e.id().event();

  Handle<EBDigiCollection> EcalDigiEB;

  e.getByToken( EBdigiCollection_ , EcalDigiEB );

  //Return if no Barrel data 
  if( !EcalDigiEB.isValid() ) return;

  // BARREL

  // loop over Digis

  const EBDigiCollection * barrelDigi = EcalDigiEB.product () ;

  std::vector<double> ebAnalogSignal ;
  std::vector<double> ebADCCounts ;
  std::vector<double> ebADCGains ;
  ebAnalogSignal.reserve(EBDataFrame::MAXSAMPLES);
  ebADCCounts.reserve(EBDataFrame::MAXSAMPLES);
  ebADCGains.reserve(EBDataFrame::MAXSAMPLES);

  int nDigis = 0;
  
  for (unsigned int digis=0; digis<EcalDigiEB->size(); ++digis) 
    {
      
      EBDataFrame ebdf = (*barrelDigi)[digis];
      int nrSamples = ebdf.size();
      
      EBDetId ebid = ebdf.id () ;

      nDigis++;
      if (meEBDigiOccupancy_) meEBDigiOccupancy_->Fill( ebid.iphi(), ebid.ieta() );
      
      double Emax = 0. ;
      int Pmax = 0 ;
      double pedestalPreSample = 0.;
      double pedestalPreSampleAnalog = 0.;
      int countsAfterGainSwitch = -1;
      double higherGain = 1.;
      int higherGainSample = 0;
      
      for (int sample = 0 ; sample < nrSamples; ++sample) {
        ebAnalogSignal[sample] = 0.;
        ebADCCounts[sample] = 0.;
        ebADCGains[sample] = 0.;
      }
      
      for (int sample = 0 ; sample < nrSamples; ++sample) 
        {
      EcalMGPASample thisSample = ebdf[sample];
      
          ebADCCounts[sample] = (thisSample.adc());
          ebADCGains[sample]  = (thisSample.gainId());
          ebAnalogSignal[sample] = (ebADCCounts[sample]*gainConv_[(int)ebADCGains[sample]]*barrelADCtoGeV_);
      
          if (Emax < ebAnalogSignal[sample] ) {
            Emax = ebAnalogSignal[sample] ;
            Pmax = sample ;
          }
      
          if ( sample < 3 ) {
            pedestalPreSample += ebADCCounts[sample] ;
            pedestalPreSampleAnalog += ebADCCounts[sample]*gainConv_[(int)ebADCGains[sample]]*barrelADCtoGeV_ ;
          }
      
          if ( sample > 0 && ( ((ebADCGains[sample] > ebADCGains[sample-1]) && (ebADCGains[sample-1]!=0)) || (countsAfterGainSwitch<0 && ebADCGains[sample]==0)) ) {
            higherGain = ebADCGains[sample];
            higherGainSample = sample;
            countsAfterGainSwitch = 1;
          }
      
          if ( (higherGain > 1 && (higherGainSample != sample) && (ebADCGains[sample] == higherGain)) || (higherGain==3 && (higherGainSample != sample) && (ebADCGains[sample]==0)) || (higherGain==0 && (higherGainSample != sample) && ((ebADCGains[sample] == 3) || (ebADCGains[sample]==0))) ) { countsAfterGainSwitch++ ; }
        }
      
      pedestalPreSample /= 3. ; 
      pedestalPreSampleAnalog /= 3. ; 
      
      LogDebug("DigiInfo") << "Barrel Digi for EBDetId = " << ebid.rawId() << " eta,phi " << ebid.ieta() << " " << ebid.iphi() ;
      for ( int i = 0; i < 10 ; i++ ) {
    LogDebug("DigiInfo") << "sample " << i << " ADC = " << ebADCCounts[i] << " gain = " << ebADCGains[i] << " Analog = " << ebAnalogSignal[i];  
      }
      LogDebug("DigiInfo") << "Maximum energy = " << Emax << " in sample " << Pmax << " Pedestal from pre-sample = " << pedestalPreSampleAnalog;
      if ( countsAfterGainSwitch > 0 ) LogDebug("DigiInfo") << "Counts after switch " << countsAfterGainSwitch;
      
      if ( countsAfterGainSwitch > 0 && countsAfterGainSwitch < 5 ) {
    edm::LogWarning("DigiWarning") << "Wrong number of counts after gain switch before next switch! " << countsAfterGainSwitch ;
    for ( int i = 0; i < 10 ; i++ ) {
      edm::LogWarning("DigiWarning") << "sample " << i << " ADC = " << ebADCCounts[i] << " gain = " << ebADCGains[i] << " Analog = " << ebAnalogSignal[i];
    } 
      }
      
      for ( int i = 0 ; i < 10 ; i++ ) {
    if (meEBDigiADCGlobal_ && (Emax-pedestalPreSampleAnalog*gainConv_[(int)ebADCGains[Pmax]]) > 100.*barrelADCtoGeV_) meEBDigiADCGlobal_->Fill( i , ebAnalogSignal[i] ) ;
    if (meEBDigiADCAnalog_[i]) meEBDigiADCAnalog_[i]->Fill( ebAnalogSignal[i] ) ;
    
    if ( ebADCGains[i] == 0) {
      if (meEBDigiADCgS_[i]) meEBDigiADCgS_[i]->Fill( ebADCCounts[i] ) ;
    }
    else if ( ebADCGains[i] == 3 ) {
      if (meEBDigiADCg1_[i]) meEBDigiADCg1_[i]->Fill( ebADCCounts[i] ) ;
    }
    else if ( ebADCGains[i] == 2 ) {
      if (meEBDigiADCg6_[i]) meEBDigiADCg6_[i]->Fill( ebADCCounts[i] ) ;
    }
    else if ( ebADCGains[i] == 1 ) {
      if (meEBDigiADCg12_[i]) meEBDigiADCg12_[i]->Fill( ebADCCounts[i] ) ;
    }
    if (meEBDigiGain_[i]) meEBDigiGain_[i]->Fill( ebADCGains[i] ) ;
      }
      
      if (meEBPedestal_) meEBPedestal_->Fill ( pedestalPreSample ) ;
      if (meEBMaximumgt10ADC_ && (Emax-pedestalPreSampleAnalog*gainConv_[(int)ebADCGains[Pmax]]) > 10.*barrelADCtoGeV_) meEBMaximumgt10ADC_->Fill( Pmax ) ;
      if (meEBMaximumgt100ADC_ && (Emax-pedestalPreSampleAnalog*gainConv_[(int)ebADCGains[Pmax]]) > 100.*barrelADCtoGeV_) meEBMaximumgt100ADC_->Fill( Pmax ) ;
      if (meEBnADCafterSwitch_) meEBnADCafterSwitch_->Fill( countsAfterGainSwitch ) ;
      
    }
  
  if ( meEBDigiMultiplicity_ ) meEBDigiMultiplicity_->Fill(nDigis);
    
}

void  EcalBarrelDigisValidation::checkCalibrations(edm::EventSetup const & eventSetup) 
  {
    
  // ADC -> GeV Scale
  edm::ESHandle<EcalADCToGeVConstant> pAgc;
  eventSetup.get<EcalADCToGeVConstantRcd>().get(pAgc);
  const EcalADCToGeVConstant* agc = pAgc.product();
  
  EcalMGPAGainRatio * defaultRatios = new EcalMGPAGainRatio();

  gainConv_[1] = 1.;
  gainConv_[2] = defaultRatios->gain12Over6() ;
  gainConv_[3] = gainConv_[2]*(defaultRatios->gain6Over1()) ;
  gainConv_[0] = gainConv_[2]*(defaultRatios->gain6Over1()) ;

  LogDebug("EcalDigi") << " Gains conversions: " << "\n" << " g0 = " << gainConv_[0] << "\n" << " g1 = " << gainConv_[1] << "\n" << " g2 = " << gainConv_[2]  << "\n" << " g3 = " << gainConv_[3]; 

  delete defaultRatios;

  const double barrelADCtoGeV_  = agc->getEBValue();
  LogDebug("EcalDigi") << " Barrel GeV/ADC = " << barrelADCtoGeV_;
  const double endcapADCtoGeV_ = agc->getEEValue();
  LogDebug("EcalDigi") << " Endcap GeV/ADC = " << endcapADCtoGeV_;

}