#include <ESDigitizer.h>

Inheritance diagram for ESDigitizer:

Classes
class	Triplet
Public Types
typedef ESDigitizerTraits::ElectronicsSim	ElectronicsSim
Public Member Functions
	ESDigitizer (EcalHitResponse hitResponse, ElectronicsSim electronicsSim, bool addNoise)
virtual void	run (MixCollection< PCaloHit > &input, DigiCollection &output)
	turns hits into digis
void	setDetIds (const std::vector< DetId > &detIds)
	tell the digitizer which cells exist; cannot change during a run
void	setGain (const int gain)
virtual	~ESDigitizer ()
Private Member Functions
void	createNoisyList (std::vector< DetId > &abThreshCh)
Private Attributes
const std::vector< DetId > *	m_detIds
CLHEP::HepRandomEngine *	m_engine
int	m_ESGain
double	m_histoBin
double	m_histoInf
double	m_histoWid
double	m_meanNoisy
CLHEP::RandFlat *	m_ranFlat
CLHEP::RandGeneral *	m_ranGeneral
CLHEP::RandPoissonQ *	m_ranPois
std::vector< Triplet >	m_trip

Detailed Description

Definition at line 15 of file ESDigitizer.h.

Member Typedef Documentation

typedef ESDigitizerTraits::ElectronicsSim ESDigitizer::ElectronicsSim

Reimplemented from EcalTDigitizer< ESDigitizerTraits >.

Definition at line 19 of file ESDigitizer.h.

Constructor & Destructor Documentation

ESDigitizer::ESDigitizer	(	EcalHitResponse *	hitResponse,
		ElectronicsSim *	electronicsSim,
		bool	addNoise
	)

Definition at line 15 of file ESDigitizer.cc.

References Exception, edm::Service< T >::isAvailable(), m_engine, and m_trip.

                                                                   :
   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 ( ) [virtual]

Definition at line 44 of file ESDigitizer.cc.

References m_ranFlat, m_ranGeneral, and m_ranPois.

{
   delete m_ranGeneral ;
   delete m_ranPois    ;
   delete m_ranFlat    ;
}

Member Function Documentation

void ESDigitizer::createNoisyList ( std::vector< DetId > & abThreshCh ) [private]

Definition at line 235 of file ESDigitizer.cc.

References spr::find(), i, m_detIds, m_ranFlat, and m_ranPois.

Referenced by run().

{
   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 ) ;
   }
}

void ESDigitizer::run	(	MixCollection< PCaloHit > &	input,
		DigiCollection &	output
	)		`[virtual]`

turns hits into digis

Reimplemented from EcalTDigitizer< ESDigitizerTraits >.

Definition at line 189 of file ESDigitizer.cc.

References EcalTDigitizer< ESDigitizerTraits >::addNoise(), createNoisyList(), EcalTDigitizer< ESDigitizerTraits >::elecSim(), EcalHitResponse::findDetId(), EcalTDigitizer< ESDigitizerTraits >::hitResponse(), m_detIds, m_engine, m_histoInf, m_histoWid, m_meanNoisy, m_ranFlat, m_ranGeneral, m_ranPois, m_trip, ESDigiCollection::push_back(), edm::DataFrameContainer::reserve(), and CaloTSamplesBase< Ttype >::zero().

Referenced by EcalDigiProducer::produce().

{
    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 ) ;  
       }
    }
}

void ESDigitizer::setDetIds ( const std::vector< DetId > & detIds )

tell the digitizer which cells exist; cannot change during a run

Definition at line 53 of file ESDigitizer.cc.

References m_detIds.

Referenced by EcalDigiProducer::updateGeometry().

{
   assert( 0       == m_detIds ||
           &detIds == m_detIds    ) ; // sanity check; don't allow to change midstream
   m_detIds = &detIds ;
}

void ESDigitizer::setGain ( const int gain )

Definition at line 61 of file ESDigitizer.cc.

References EcalTDigitizer< ESDigitizerTraits >::addNoise(), dtNoiseDBValidation_cfg::cerr, Exception, jetmet_cfg::histofile, edm::errors::InvalidReference, m_detIds, m_engine, m_ESGain, m_histoBin, m_histoInf, m_histoWid, m_meanNoisy, m_ranFlat, m_ranGeneral, m_ranPois, m_trip, query::result, and ntuplemaker::status.

Referenced by EcalDigiProducer::checkCalibrations().

{
   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();
         }
      }
   }
}