CMSSW_4_4_3_patch1/src/SimCalorimetry/EcalSimAlgos/src/ESDigitizer.cc

Go to the documentation of this file.

#include "SimCalorimetry/EcalSimAlgos/interface/ESDigitizer.h"
#include "SimCalorimetry/EcalSimAlgos/interface/ESElectronicsSimFast.h"
#include "SimCalorimetry/EcalSimAlgos/interface/ESHitResponse.h"
#include "SimDataFormats/CaloHit/interface/PCaloHit.h"

#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
#include "CLHEP/Random/RandGeneral.h"
#include "CLHEP/Random/RandPoissonQ.h"
#include "CLHEP/Random/RandFlat.h"

#include <gsl/gsl_sf_erf.h>
#include <gsl/gsl_sf_result.h>

ESDigitizer::ESDigitizer( EcalHitResponse*      hitResponse    , 
                          ESElectronicsSimFast* electronicsSim ,
                          bool                  addNoise         ) :
   EcalTDigitizer< ESDigitizerTraits >( hitResponse, electronicsSim, addNoise ) ,
   m_detIds      ( 0              ) ,
   m_engine      ( 0              ) ,
   m_ranGeneral  ( 0              ) ,
   m_ranPois     ( 0              ) ,
   m_ranFlat     ( 0              ) ,
   m_ESGain      ( 0              ) ,
   m_histoBin    ( 0              ) ,
   m_histoInf    ( 0              ) ,
   m_histoWid    ( 0              ) ,
   m_meanNoisy   ( 0              ) ,
   m_trip        (                )
{
   m_trip.reserve( 2500 ) ; 

   edm::Service<edm::RandomNumberGenerator> rng ;
   if( !rng.isAvailable() )
   {
      throw cms::Exception( "Configuration" )
         << "ESDigitizer requires the RandomNumberGeneratorService\n"
         "which is not present in the configuration file.  You must add the service\n"
         "in the configuration file or remove the modules that require it.";
   }
   m_engine = &rng->getEngine() ;
}

ESDigitizer::~ESDigitizer() 
{
   delete m_ranGeneral ;
   delete m_ranPois    ;
   delete m_ranFlat    ;
}

void 
ESDigitizer::setDetIds( const std::vector<DetId>& detIds )
{
   assert( 0       == m_detIds ||
           &detIds == m_detIds    ) ; // sanity check; don't allow to change midstream
   m_detIds = &detIds ;
}

void 
ESDigitizer::setGain( const int gain ) 
{
   if( 0 != m_ESGain )
   {
      assert( gain == m_ESGain ) ; // only allow one value
   }
   else
   {
      assert( 0 != m_detIds &&
              0 != m_detIds->size() ) ; // detIds must already be set as we need size

      assert( 1 == gain ||
              2 == gain    ) ; // legal values
      
      m_ESGain = gain ;
      
      if( addNoise() ) 
      {
         assert( 0 != m_engine ) ; // sanity check

         double zsThresh ( 0. ) ;
         std::string refFile ;

         if( 1 == m_ESGain ) 
         {
            zsThresh = 3 ;
            refFile = "SimCalorimetry/EcalSimProducers/data/esRefHistosFile_LG.txt";
         }
         else
         {
            zsThresh = 4 ;
            refFile = "SimCalorimetry/EcalSimProducers/data/esRefHistosFile_HG.txt";
         }

         gsl_sf_result result ;
         int status  = gsl_sf_erf_Q_e( zsThresh, &result ) ;
         if( status != 0 ) std::cerr << "ESDigitizer::could not compute gaussian tail probability for the threshold chosen" << std::endl ;

         const double probabilityLeft ( result.val ) ;
         m_meanNoisy = probabilityLeft * m_detIds->size() ;

         m_ranPois = new CLHEP::RandPoissonQ( *m_engine, m_meanNoisy      ) ;
         m_ranFlat = new CLHEP::RandFlat(     *m_engine, m_detIds->size() ) ;

         std::ifstream histofile ( edm::FileInPath( refFile ).fullPath().c_str() ) ;
         if( !histofile.good() )
         { 
            throw edm::Exception(edm::errors::InvalidReference,"NullPointer")
               << "Reference histos file not opened" ;
         }
         else
         {
            // number of bins
            char buffer[200] ;
            int thisLine = 0 ;
            while( 0 == thisLine ) 
            {
               histofile.getline( buffer, 200 ) ;
               if( !strstr(buffer,"#")  && 
                   !(strspn(buffer," ") == strlen(buffer) ) )
               {        
                  float histoBin ; 
                  sscanf( buffer, "%f" , &histoBin ) ; 
                  m_histoBin = (double) histoBin ;
                  ++thisLine ;
               }
            }
            const uint32_t histoBin1 ( (int) m_histoBin    ) ;
            const uint32_t histoBin2 ( histoBin1*histoBin1 ) ;

            double t_histoSup ( 0 ) ;

            std::vector<double> t_refHistos ;
            t_refHistos.reserve( 2500 ) ;

            int thisBin = -2 ;
            while( !( histofile.eof() ) )
            {
               histofile.getline( buffer, 200 ) ;
               if( !strstr( buffer, "#" ) &&
                   !( strspn( buffer, " " ) == strlen( buffer ) ) )
               {
                  if( -2 == thisBin )
                  {
                     float histoInf ;
                     sscanf( buffer, "%f" , &histoInf ) ;
                     m_histoInf = (double) histoInf ;
                  }
                  if( -1 == thisBin  )
                  {
                     float histoSup ;
                     sscanf( buffer, "%f" , &histoSup ) ;
                     t_histoSup = (double) histoSup ;
                  }
                  if( 0 <= thisBin )
                  { 
                     float refBin ; 
                     sscanf( buffer, "%f", &refBin ) ;
                     if( 0.5 < refBin ) 
                     {
                        t_refHistos.push_back( (double) refBin ) ;
                        const uint32_t i2 ( thisBin/histoBin2 ) ;
                        const uint32_t off ( i2*histoBin2 ) ;
                        const uint32_t i1 ( ( thisBin - off )/histoBin1 ) ;
                        const uint32_t i0 ( thisBin - off - i1*histoBin1 ) ;
                        m_trip.push_back( Triplet( i0, i1, i2 ) ) ;
                     }
                  }
                  ++thisBin ;
               }
            }
            m_histoWid = ( t_histoSup - m_histoInf )/m_histoBin ;

            m_histoInf -= 1000. ;

            // creating the reference distribution to extract random numbers
            m_ranGeneral = new CLHEP::RandGeneral( *m_engine            ,
                                                   &t_refHistos.front() ,
                                                   t_refHistos.size()   ,
                                                   0             ) ;
            histofile.close();
         }
      }
   }
}

void 
ESDigitizer::run( MixCollection<PCaloHit>& input  ,
                  ESDigiCollection&        output   )
{
    assert( 0 != m_detIds         &&
            0 != m_detIds->size() &&
            ( !addNoise()         ||
              ( 0 != m_engine     &&
                0 != m_ranPois    &&
                0 != m_ranFlat    &&
                0 != m_ranGeneral        ) ) ) ; // sanity check

    // reserve space for how many digis we expect, with some cushion
    output.reserve( 2*( (int) m_meanNoisy ) + hitResponse()->samplesSize() ) ;

    EcalTDigitizer< ESDigitizerTraits >::run( input, output ) ;

    // random generation of channel above threshold
    std::vector<DetId> abThreshCh ;
    if( addNoise() ) createNoisyList( abThreshCh ) ;

    // first make a raw digi for every cell where we have noise
    for( std::vector<DetId>::const_iterator idItr ( abThreshCh.begin() ) ;
         idItr != abThreshCh.end(); ++idItr ) 
    {
       if( hitResponse()->findDetId( *idItr )->zero() ) // only if no true hit!
       {
          ESHitResponse::ESSamples analogSignal ( *idItr, 3 ) ; // space for the noise hit
          uint32_t myBin ( (uint32_t) m_trip.size()*m_ranGeneral->fire() ) ;
          if( myBin == m_trip.size() ) --myBin ; // guard against roundup
          assert( myBin < m_trip.size() ) ;
          const Triplet& trip ( m_trip[ myBin ] ) ;
          analogSignal[ 0 ] = m_histoInf + m_histoWid*trip.first  ;
          analogSignal[ 1 ] = m_histoInf + m_histoWid*trip.second ;
          analogSignal[ 2 ] = m_histoInf + m_histoWid*trip.third  ;
          ESDataFrame digi( *idItr ) ;
          const_cast<ESElectronicsSimFast*>(elecSim())->
             analogToDigital( analogSignal ,
                              digi         ,
                              true           ) ;        
          output.push_back( digi ) ;  
       }
    }
}

// preparing the list of channels where the noise has to be generated
void 
ESDigitizer::createNoisyList( std::vector<DetId>& abThreshCh )
{
   const unsigned int nChan ( m_ranPois->fire() ) ;
   abThreshCh.reserve( nChan ) ;

   for( unsigned int i ( 0 ) ; i != nChan ; ++i )
   {
      std::vector<DetId>::const_iterator idItr ( abThreshCh.end() ) ;
      uint32_t iChan ( 0 ) ;
      DetId id ;
      do 
      {
         iChan = (uint32_t) m_ranFlat->fire() ;
         if( iChan == m_detIds->size() ) --iChan ; //protect against roundup at end
         assert( m_detIds->size() > iChan ) ;      // sanity check
         id = (*m_detIds)[ iChan ] ;
         idItr = find( abThreshCh.begin() ,
                       abThreshCh.end()   ,
                       id                  ) ;
      }
      while( idItr != abThreshCh.end() ) ;

      abThreshCh.push_back( id ) ;
   }
}