PedestalsTask.cc

Go to the documentation of this file.

#include "DQM/SiStripCommissioningSources/interface/PedestalsTask.h"
#include "DataFormats/SiStripCommon/interface/SiStripConstants.h"
#include "DataFormats/SiStripCommon/interface/SiStripHistoTitle.h"
#include "DQMServices/Core/interface/DQMStore.h"
#include "DQM/SiStripCommon/interface/ExtractTObject.h"
#include "DQM/SiStripCommon/interface/UpdateTProfile.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <algorithm>
#include <cmath>

using namespace sistrip;

// -----------------------------------------------------------------------------
//
PedestalsTask::PedestalsTask( DQMStore* dqm,
                    const FedChannelConnection& conn ) :
  CommissioningTask( dqm, conn, "PedestalsTask" ),
  peds_(),
  cm_()
{
  LogTrace(mlDqmSource_)
    << "[PedestalsTask::" << __func__ << "]"
    << " Constructing object...";
}

// -----------------------------------------------------------------------------
//
PedestalsTask::~PedestalsTask() {
  LogTrace(mlDqmSource_)
    << "[PedestalsTask::" << __func__ << "]"
    << " Destructing object...";
}

// -----------------------------------------------------------------------------
//
void PedestalsTask::book() {
  LogTrace(mlDqmSource_) << "[PedestalsTask::" << __func__ << "]";
  
  uint16_t nbins;
  std::string title;
  std::string extra_info;
  peds_.resize(2);
  nbins = 256;
  
  // Pedestals histogram
  extra_info = sistrip::extrainfo::pedestals_; 
  peds_[0].isProfile_ = true;
  
  title = SiStripHistoTitle( sistrip::EXPERT_HISTO, 
                 sistrip::PEDESTALS, 
                 sistrip::FED_KEY, 
                 fedKey(),
                 sistrip::LLD_CHAN, 
                 connection().lldChannel(),
                 extra_info ).title();
  
  peds_[0].histo( dqm()->bookProfile( title, title, 
                      nbins, -0.5, nbins*1.-0.5,
                      1025, 0., 1025. ) );
  
  peds_[0].vNumOfEntries_.resize(nbins,0);
  peds_[0].vSumOfContents_.resize(nbins,0);
  peds_[0].vSumOfSquares_.resize(nbins,0);

  // Noise histogram
  extra_info = sistrip::extrainfo::noise_; 
  peds_[1].isProfile_ = true;
  
  title = SiStripHistoTitle( sistrip::EXPERT_HISTO, 
                 sistrip::PEDESTALS, 
                 sistrip::FED_KEY, 
                 fedKey(),
                 sistrip::LLD_CHAN, 
                 connection().lldChannel(),
                 extra_info ).title();
  
  peds_[1].histo( dqm()->bookProfile( title, title, 
                      nbins, -0.5, nbins*1.-0.5,
                      1025, 0., 1025. ) );
  
  peds_[1].vNumOfEntries_.resize(nbins,0);
  peds_[1].vSumOfContents_.resize(nbins,0);
  peds_[1].vSumOfSquares_.resize(nbins,0);

  // Common mode histograms
  cm_.resize(2);
  nbins = 1024;
  for ( uint16_t iapv = 0; iapv < 2; iapv++ ) { 
    
    title = SiStripHistoTitle( sistrip::EXPERT_HISTO, 
                   sistrip::PEDESTALS, 
                   sistrip::FED_KEY, 
                   fedKey(),
                   sistrip::APV, 
                   connection().i2cAddr(iapv),
                   sistrip::extrainfo::commonMode_ ).title();
    
    cm_[iapv].histo( dqm()->book1D( title, title, nbins, -0.5, nbins*1.-0.5 ) );
    cm_[iapv].isProfile_ = false;
    
    cm_[iapv].vNumOfEntries_.resize(nbins,0);
    cm_[iapv].vNumOfEntries_.resize(nbins,0);
    
  }
  
}

// -----------------------------------------------------------------------------
//
void PedestalsTask::fill( const SiStripEventSummary& summary,
              const edm::DetSet<SiStripRawDigi>& digis ) {
  
  if ( digis.data.size() != peds_[0].vNumOfEntries_.size() ) {
    edm::LogWarning(mlDqmSource_)
      << "[PedestalsTask::" << __func__ << "]"
      << " Unexpected number of digis: " 
      << digis.data.size(); 
    return;
  }
  
  // Check number of digis
  uint16_t nbins = peds_[0].vNumOfEntries_.size();
  if ( digis.data.size() < nbins ) { nbins = digis.data.size(); }

  //@@ Inefficient!!!
  uint16_t napvs = nbins / 128;
  std::vector<uint32_t> cm; cm.resize(napvs,0);
  
  // Calc common mode for both APVs
  std::vector<uint16_t> adc;
  for ( uint16_t iapv = 0; iapv < napvs; iapv++ ) { 
    adc.clear(); adc.reserve(128);
    for ( uint16_t ibin = 0; ibin < 128; ibin++ ) { 
      if ( (iapv*128)+ibin < nbins ) { 
    adc.push_back( digis.data[(iapv*128)+ibin].adc() );  //@@ VIRGIN RAW DATA (MUX, APV READOUT)
      }
    }
    sort( adc.begin(), adc.end() ); 
    uint16_t index = adc.size()%2 ? adc.size()/2 : adc.size()/2-1;
    if ( !adc.empty() ) { cm[iapv] = static_cast<uint32_t>( adc[index] ); }
  }
  
  for ( uint16_t ibin = 0; ibin < nbins; ibin++ ) {
    updateHistoSet( peds_[0], ibin, digis.data[ibin].adc() ); // peds and raw noise
    float diff = static_cast<float>( digis.data[ibin].adc() ) - static_cast<float>( cm[ibin/128] );
    updateHistoSet( peds_[1], ibin, diff ); // residuals and real noise
  }
  
  if ( cm.size() < cm_.size() ) {
    edm::LogWarning(mlDqmSource_)
      << "[PedestalsTask::" << __func__ << "]"
      << " Fewer CM values than expected: " << cm.size();
  }
  
  updateHistoSet( cm_[0], cm[0] );
  updateHistoSet( cm_[1], cm[1] );
  
}

// -----------------------------------------------------------------------------
//
void PedestalsTask::update() {
  
  // Pedestals 
  updateHistoSet( peds_[0] );
  
  // Noise (cannot use HistoSet directly, as want to plot noise as "contents", not "error")
  TProfile* histo = ExtractTObject<TProfile>().extract( peds_[1].histo() );
  for ( uint16_t ii = 0; ii < peds_[1].vNumOfEntries_.size(); ++ii ) {
    
    float mean    = 0.;
    float spread  = 0.;
    float entries = peds_[1].vNumOfEntries_[ii];
    if ( entries > 0. ) {
      mean = peds_[1].vSumOfContents_[ii] / entries;
      spread = sqrt( peds_[1].vSumOfSquares_[ii] / entries - mean * mean ); 
    }
    
    float noise = spread;
    float error = 0; // sqrt(entries) / entries;
    
    UpdateTProfile::setBinContent( histo, ii+1, entries, noise, error );
    
  }
  
  // Common mode
  updateHistoSet( cm_[0] );
  updateHistoSet( cm_[1] );
  
}