EcalCoder.cc

Go to the documentation of this file.

#include "SimCalorimetry/EcalSimAlgos/interface/EcalCoder.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
#include "DataFormats/EcalDigi/interface/EcalMGPASample.h"
#include "DataFormats/EcalDigi/interface/EcalDataFrame.h"

#include <iostream>

//#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"


EcalCoder::EcalCoder( bool                  addNoise     , 
              bool                  PreMix1      ,
              EcalCoder::Noisifier* ebCorrNoise0 ,
              EcalCoder::Noisifier* eeCorrNoise0 ,
              EcalCoder::Noisifier* ebCorrNoise1 ,
              EcalCoder::Noisifier* eeCorrNoise1 ,
              EcalCoder::Noisifier* ebCorrNoise2 ,
              EcalCoder::Noisifier* eeCorrNoise2   ) :
   m_peds        (           nullptr ) ,
   m_gainRatios  (           nullptr ) ,
   m_intercals   (           nullptr ) ,
   m_maxEneEB    (      1668.3 ) , // 4095(MAXADC)*12(gain 2)*0.035(GeVtoADC)*0.97
   m_maxEneEE    (      2859.9 ) , // 4095(MAXADC)*12(gain 2)*0.060(GeVtoADC)*0.97
   m_addNoise    ( addNoise    ) ,
   m_PreMix1     ( PreMix1     ) 
   {
   m_ebCorrNoise[0] = ebCorrNoise0 ;
   assert( nullptr != m_ebCorrNoise[0] ) ;
   m_eeCorrNoise[0] = eeCorrNoise0 ;
   m_ebCorrNoise[1] = ebCorrNoise1 ;
   m_eeCorrNoise[1] = eeCorrNoise1 ;
   m_ebCorrNoise[2] = ebCorrNoise2 ;
   m_eeCorrNoise[2] = eeCorrNoise2 ;
}  

EcalCoder::~EcalCoder()
{
}

void 
EcalCoder::setFullScaleEnergy( double EBscale ,
                   double EEscale   )
{
   m_maxEneEB = EBscale ;
   m_maxEneEE = EEscale ;
}


void  
EcalCoder::setPedestals( const EcalPedestals* pedestals ) 
{
   m_peds = pedestals ;
}

void  
EcalCoder::setGainRatios( const EcalGainRatios* gainRatios ) 
{
   m_gainRatios = gainRatios ; 
}

void 
EcalCoder::setIntercalibConstants( const EcalIntercalibConstantsMC* ical ) 
{
   m_intercals = ical ;
}

double 
EcalCoder::fullScaleEnergy( const DetId & detId ) const 
{
   return detId.subdetId() == EcalBarrel ? m_maxEneEB : m_maxEneEE ;
}

void 
EcalCoder::analogToDigital( CLHEP::HepRandomEngine* engine,
                            const EcalSamples& clf ,
                EcalDataFrame&     df    ) const 
{
   df.setSize( clf.size() ) ;
   encode( clf, df, engine );
/*   std::cout<<"   **Id=" ;
   if( CaloGenericDetId( clf.id() ).isEB() )
   {
      std::cout<<EBDetId( clf.id() ) ;
   }
   else
   {
      std::cout<<EEDetId( clf.id() ) ;
   }
   std::cout<<", size="<<df.size();
   for( int i ( 0 ) ; i != df.size() ; ++i )
   {
      std::cout<<", "<<df[i];
   }
   std::cout<<std::endl ;*/
}

void 
EcalCoder::encode( const EcalSamples& ecalSamples , 
           EcalDataFrame&     df,
                   CLHEP::HepRandomEngine* engine ) const
{
   assert( nullptr != m_peds ) ;

   const unsigned int csize ( ecalSamples.size() ) ;

  
   DetId detId = ecalSamples.id();             
   double Emax = fullScaleEnergy(detId);       

   //....initialisation
   if ( ecalSamples[5] > 0. ) LogDebug("EcalCoder") << "Input caloSample" << "\n" << ecalSamples;
  
   double LSB[NGAINS+1];
   double pedestals[NGAINS+1];
   double widths[NGAINS+1];
   double gains[NGAINS+1];
   int    maxADC[NGAINS+1];
   double trueRMS[NGAINS+1];

   double icalconst = 1. ;
   findIntercalibConstant( detId, icalconst );

   for( unsigned int igain ( 0 ); igain <= NGAINS ; ++igain ) 
   {
      // fill in the pedestal and width
      findPedestal( detId ,
            igain , 
            pedestals[igain] , 
            widths[igain]      ) ;

      if( 0 < igain ) 
     trueRMS[igain] = std::sqrt( widths[igain]*widths[igain] - 1./12. ) ;

      // set nominal value first
      findGains( detId , 
         gains  );               

      LSB[igain] = 0.;
      if ( igain > 0 ) LSB[igain]= Emax/(MAXADC*gains[igain]);
      maxADC[igain] = ADCGAINSWITCH;                   // saturation at 4080 for middle and high gains x6 & x12
      if ( igain == NGAINS ) maxADC[igain] = MAXADC;   // saturation at 4095 for low gain x1 
   }

   CaloSamples noiseframe[] = { CaloSamples( detId , csize ) ,
                CaloSamples( detId , csize ) ,
                CaloSamples( detId , csize )   } ;

   const Noisifier* noisy[3] = { ( nullptr == m_eeCorrNoise[0]          ||
                   EcalBarrel == detId.subdetId()    ?
                   m_ebCorrNoise[0] :
                   m_eeCorrNoise[0]                  ) ,
                 ( EcalBarrel == detId.subdetId() ?
                   m_ebCorrNoise[1] :
                   m_eeCorrNoise[1]                  ) ,
                 ( EcalBarrel == detId.subdetId() ?
                   m_ebCorrNoise[2] :
                   m_eeCorrNoise[2]                  )   } ;

   if( m_addNoise )
   {
     noisy[0]->noisify( noiseframe[0], engine ) ; // high gain
      if( nullptr == noisy[1] ) noisy[0]->noisify( noiseframe[1] ,
                                             engine,
                         &noisy[0]->vecgau() ) ; // med 
      if( nullptr == noisy[2] ) noisy[0]->noisify( noiseframe[2] ,
                                             engine,
                         &noisy[0]->vecgau() ) ; // low
   }


   //   std::cout << " intercal, LSBs, gains " << icalconst << " " << LSB[0] << " " << LSB[1] << " " << gains[0] << " " << gains[1] << " " << Emax <<  std::endl;

   int wait = 0 ;
   int gainId = 0 ;
   bool isSaturated = false;

   for( unsigned int i ( 0 ) ; i != csize ; ++i )
   {    
      bool done ( false ) ;
      int adc = MAXADC ;
      if( 0 == wait ) gainId = 1 ;

      // see which gain bin it fits in
      int igain ( gainId - 1 ) ;
      do
      {
     ++igain ;

//   if( igain != gainId ) std::cout <<"$$$$ Gain switch from " << gainId
//                   <<" to "<< igain << std::endl ;

     if( 1 != igain                    &&   // not high gain
         m_addNoise                    &&   // want to add noise
         nullptr != noisy[igain-1]           &&   // exists
         noiseframe[igain-1].isBlank()    ) // not already done
     {
        noisy[igain-1]->noisify( noiseframe[igain-1] ,
                                     engine,
                     &noisy[0]->vecgau()   ) ;
        //std::cout<<"....noisifying gain level = "<<igain<<std::endl ;
     }
    
     double signal;

     if(!m_PreMix1) {

       // noiseframe filled with zeros if !m_addNoise
       const double asignal ( pedestals[igain] +
                   ecalSamples[i] /( LSB[igain]*icalconst ) +
                   trueRMS[igain]*noiseframe[igain-1][i]      ) ;
       signal = asignal;
     }
     else {  // Any changes made here must be reverse-engineered in EcalSignalGenerator!

           if( igain == 1) {
             const double asignal ( ecalSamples[i]*1000. );  // save low level info                   
             signal = asignal;
           }
           else if( igain == 2) {
             const double asignal ( ecalSamples[i]/( LSB[1]*icalconst ));
             signal = asignal;
           }
           else if( igain == 3) {   // bet that no pileup hit has an energy over Emax/2             
             const double asignal ( ecalSamples[i]/( LSB[2]*icalconst ) );
             signal = asignal;
           }
       else { //not sure we ever get here at gain=0, but hit wil be saturated anyway
         const double asignal ( ecalSamples[i]/( LSB[3]*icalconst ) ); // just calculate something
             signal = asignal;
       }
       // old version
       //const double asignal ( // no pedestals for pre-mixing
       //            ecalSamples[i] /( LSB[igain]*icalconst ) );
       //signal = asignal;
     }

     //  std::cout << " " << ecalSamples[i] << " " << noiseframe[igain-1][i] << std::endl;


     const int isignal ( signal ) ;
     const int tmpadc ( signal - (double)isignal < 0.5 ?
                isignal : isignal + 1 ) ;
     // LogDebug("EcalCoder") << "DetId " << detId.rawId() << " gain " << igain << " caloSample " 
     //                       <<  ecalSamples[i] << " pededstal " << pedestals[igain] 
     //                       << " noise " << widths[igain] << " conversion factor " << LSB[igain] 
     //                       << " result (ped,tmpadc)= " << ped << " " << tmpadc;
         
     if( tmpadc <= maxADC[igain] ) 
     {
        adc = tmpadc;
        done = true ;
     }
      }
      while( !done       &&
         igain < 3    ) ;

      if (igain == 1 ) 
      {
         wait = 0 ;
         gainId = igain ;
      }
      else 
      {
         if (igain == gainId) --wait ;
         else 
     {
        wait = 5 ;
        gainId = igain ;
     }
      }


      // change the gain for saturation
      int storeGainId = gainId;
      if ( gainId == 3 && adc == MAXADC ) 
      {
     storeGainId = 0;
     isSaturated = true;
      }
      // LogDebug("EcalCoder") << " Writing out frame " << i << " ADC = " << adc << " gainId = " << gainId << " storeGainId = " << storeGainId ; 
     
      df.setSample( i, EcalMGPASample( adc, storeGainId ) );   
   }
   // handle saturation properly according to IN-2007/056
   // N.B. - FIXME 
   // still missing the possibility for a signal to pass the saturation threshold
   // again within the 5 subsequent samples (higher order effect).

   if ( isSaturated ) 
   {
      for (unsigned int i = 0 ; i < ecalSamples.size() ; ++i) 
      {
     if ( df.sample(i).gainId() == 0 ) 
     {
        unsigned int hyst = i+1+2;
        for ( unsigned int j = i+1; j < hyst && j < ecalSamples.size(); ++j ) 
        {
           df.setSample(j, EcalMGPASample(MAXADC, 0));   
        }
     }
      }
   }
}

void 
EcalCoder::findPedestal( const DetId & detId  , 
             int           gainId , 
             double&       ped    , 
             double&       width     ) const
{
   /*
     EcalPedestalsMapIterator mapItr 
     = m_peds->m_pedestals.find(detId.rawId());
     // should I care if it doesn't get found?
     if(mapItr == m_peds->m_pedestals.end()) {
     edm::LogError("EcalCoder") << "Could not find pedestal for " << detId.rawId() << " among the " << m_peds->m_pedestals.size();
     } else {
     EcalPedestals::Item item = mapItr->second;
   */
   
   /*   
    EcalPedestals::Item const & item = (*m_peds)(detId);
    ped = item.mean(gainId);
    width = item.rms(gainId);
   */

   EcalPedestalsMap::const_iterator itped = m_peds->getMap().find( detId );
   ped   = (*itped).mean(gainId);
   width = (*itped).rms(gainId);
  
   if ( (detId.subdetId() != EcalBarrel) && (detId.subdetId() != EcalEndcap) ) 
   { 
      edm::LogError("EcalCoder") << "Could not find pedestal for " << detId.rawId() << " among the " << m_peds->getMap().size();
   } 

  /*
    switch(gainId) {
    case 0:
      ped = 0.;
      width = 0.;
    case 1:
      ped = item.mean_x12;
      width = item.rms_x12;
      break;
    case 2:
      ped = item.mean_x6;
      width = item.rms_x6;
      break;
    case 3:
      ped = item.mean_x1;
      width = item.rms_x1;
      break;
    default:
      edm::LogError("EcalCoder") << "Bad Pedestal " << gainId;
      break;
    }
  */

   LogDebug("EcalCoder") << "Pedestals for " << detId.rawId() << " gain range " << gainId << " : \n" << "Mean = " << ped << " rms = " << width;
}


void 
EcalCoder::findGains( const DetId & detId , 
              double Gains[]        ) const
{
  /*
    EcalGainRatios::EcalGainRatioMap::const_iterator grit=m_gainRatios->getMap().find(detId.rawId());
    EcalMGPAGainRatio mgpa;
    if( grit!=m_gainRatios->getMap().end() ){
    mgpa = grit->second;
    Gains[0] = 0.;
    Gains[3] = 1.;
    Gains[2] = mgpa.gain6Over1() ;
    Gains[1] = Gains[2]*(mgpa.gain12Over6()) ;
    LogDebug("EcalCoder") << "Gains for " << detId.rawId() << "\n" << " 1 = " << Gains[1] << "\n" << " 2 = " << Gains[2] << "\n" << " 3 = " << Gains[3];
    } else {
    edm::LogError("EcalCoder") << "Could not find gain ratios for " << detId.rawId() << " among the " << m_gainRatios->getMap().size();
    }
  */

   EcalGainRatioMap::const_iterator grit = m_gainRatios->getMap().find( detId );
   Gains[0] = 0.;
   Gains[3] = 1.;
   Gains[2] = (*grit).gain6Over1();
   Gains[1] = Gains[2]*( (*grit).gain12Over6() );   
   
  
   if ( (detId.subdetId() != EcalBarrel) && (detId.subdetId() != EcalEndcap) ) 
   { 
      edm::LogError("EcalCoder") << "Could not find gain ratios for " << detId.rawId() << " among the " << m_gainRatios->getMap().size();
   }   
  
}

void 
EcalCoder::findIntercalibConstant( const DetId& detId, 
                   double&      icalconst ) const
{
   EcalIntercalibConstantMC thisconst = 1.;
   // find intercalib constant for this xtal
   const EcalIntercalibConstantMCMap &icalMap = m_intercals->getMap();
   EcalIntercalibConstantMCMap::const_iterator icalit = icalMap.find(detId);
   if( icalit!=icalMap.end() )
   {
      thisconst = (*icalit);
      if ( icalconst == 0. ) { thisconst = 1.; }
   } 
   else
   {
      edm::LogError("EcalCoder") << "No intercalib const found for xtal " << detId.rawId() << "! something wrong with EcalIntercalibConstants in your DB? ";
   }
   icalconst = thisconst;
}