LArHitEMapToDigitAlg.cxx

Go to the documentation of this file.

/*
 *   Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
 */
#include "LArHitEMapToDigitAlg.h"
#include "AthenaKernel/ITriggerTime.h"
#include "CLHEP/Random/RandGaussZiggurat.h"
#include "CaloDetDescr/CaloDetDescrManager.h"
#include "CaloIdentifier/CaloIdManager.h"
#include "CaloIdentifier/LArID.h"
#include "CaloIdentifier/LArID_Exception.h"
#include "CaloIdentifier/LArNeighbours.h"
#include "Identifier/IdentifierHash.h"
#include "LArIdentifier/LArOnlineID.h"
 
#include "AthenaKernel/RNGWrapper.h"
#include "CLHEP/Random/RandomEngine.h"
#include <CLHEP/Random/Randomize.h>
#include <GaudiKernel/StatusCode.h>
 
#include "StoreGate/WriteHandle.h"
 
using CLHEP::RandFlat;
using CLHEP::RandGaussZiggurat;
 
StatusCode LArHitEMapToDigitAlg::initialize()
{
 
  if (m_NSamples > s_MaxNSamples) {
    ATH_MSG_ERROR("Requested Nsamples " << m_NSamples << " larger than max "
                                        << s_MaxNSamples);
    return StatusCode::FAILURE;
  }
  ATH_CHECK(m_noiseKey.initialize((!m_RndmEvtOverlay || m_isMcOverlay) &&
                                  !m_pedestalNoise && m_NoiseOnOff));
 
  ATH_CHECK(m_shapeKey.initialize());
  ATH_CHECK(m_fSamplKey.initialize());
  ATH_CHECK(m_OFCKey.initialize());
  ATH_CHECK(m_pedestalKey.initialize());
  ATH_CHECK(m_autoCorrNoiseKey.initialize(!m_RndmEvtOverlay && m_NoiseOnOff) );
  ATH_CHECK(m_bcContKey.initialize());
  ATH_CHECK(m_badFebKey.initialize());
  ATH_CHECK(m_adc2mevKey.initialize());
  ATH_CHECK(m_caloMgrKey.initialize());
 
  ATH_CHECK(m_cablingKey.initialize());
  
  // helpers
  //retrieve ID helpers
  ATH_CHECK(detStore()->retrieve(m_calocell_id,"CaloCell_ID"));
 
 
  const CaloIdManager* caloIdMgr = nullptr;
  StatusCode sc = detStore()->retrieve(caloIdMgr);
  if (sc.isFailure()) {
    ATH_MSG_ERROR(" Unable to retrieve CaloIdManager from DetectoreStore");
    return StatusCode::FAILURE;
  }
  m_larem_id   = caloIdMgr->getEM_ID();
  m_larhec_id  = caloIdMgr->getHEC_ID();
  m_larfcal_id = caloIdMgr->getFCAL_ID();
 
  sc = detStore()->retrieve(m_laronline_id);
  if (sc.isFailure()) {
    ATH_MSG_ERROR(" Unable to retrieve LArOnlineId from DetectoreStore");
    return StatusCode::FAILURE;
  }
  ATH_CHECK(m_bcMask.buildBitMask(m_problemsToMask,msg()));
 
  // Services
  ATH_CHECK(m_rndmGenSvc.retrieve());
  ATH_CHECK(m_hitMapKey.initialize());
  ATH_CHECK(m_hitMapKey_DigiHSTruth.initialize(m_doDigiTruth));
 
  ATH_CHECK(m_DigitContainerName.initialize());
  ATH_CHECK(m_DigitContainerName_DigiHSTruth.initialize(m_doDigiTruth));
 
  // Check consistency of gain-ranges and gain-switching thresholds:
  std::array<std::pair<int,std::string>,4> iCaloToStr{{{EM,"EM"},{HEC,"HEC"},{FCAL,"FCAL"},{EMIW,"EMIW"}}};
 
  for (int iCalo = EM; iCalo <= EMIW; ++iCalo) {
    if ((m_gainRange[iCalo].value().first == CaloGain::LARHIGHGAIN) == (m_HighGainThresh[iCalo] != 0))
      ATH_MSG_INFO("jobO consistency check: " << iCaloToStr[iCalo] << " has " << ((m_gainRange[iCalo].value().first == CaloGain::LARHIGHGAIN) ? "" : "no ")
                                              << " HIGH gain and high Gain threshold=" << m_HighGainThresh[iCalo]);
    else {
      ATH_MSG_ERROR("jobO inconsistency! " << iCaloToStr[iCalo] << " has " << ((m_gainRange[iCalo].value().first == CaloGain::LARHIGHGAIN) ? "" : "no ")
                                           << " HIGH gain but high Gain threshold=" << m_HighGainThresh[iCalo]);
      return StatusCode::FAILURE;
    }
    if ((m_gainRange[iCalo].value().second == CaloGain::LARLOWGAIN) == (m_LowGainThresh[iCalo] <= m_maxADC))
      ATH_MSG_INFO("jobO consistency check: Calo " << iCaloToStr[iCalo] << " has " << ((m_gainRange[iCalo].value().first == CaloGain::LARHIGHGAIN) ? "" : "no ")
                                                   << "LOW gain and high Gain threshold=" << m_LowGainThresh[iCalo] << " (maxADC=" << m_maxADC << ")");
    else {
      ATH_MSG_ERROR("jobO inconsistency! Calo " << iCaloToStr[iCalo] << " has " << ((m_gainRange[iCalo].value().first == CaloGain::LARHIGHGAIN) ? "" : "no ")
                                                << "LOW gain but high Gain threshold=" << m_LowGainThresh[iCalo] << " (maxADC=" << m_maxADC << ")");
      return StatusCode::FAILURE;
    }
 
    if (m_gainRange[iCalo].value().first == m_gainRange[iCalo].value().second) {
      ATH_MSG_ERROR(" Calo " << iCaloToStr[iCalo] << " configured to have only one gain. This is not supported.");
      return Status::FAILURE;
    }
    if (m_HighGainThresh[iCalo] >= m_LowGainThresh[iCalo] ) {
      ATH_MSG_ERROR(" Calo " << iCaloToStr[iCalo] << " High gain threshold > low gain threshold! " << m_HighGainThresh[iCalo] << " >= " << m_LowGainThresh[iCalo]);
      return Status::FAILURE;  
    }
 
 
  } //end  iCalo loop
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode LArHitEMapToDigitAlg::execute(const EventContext& context) const {
 
   // load many conditions
   const LArBadFebCont* badFebs = pointerFromKey<LArBadFebCont>(context,m_badFebKey);
   SG::ReadCondHandle<LArOnOffIdMapping> cablingHdl{m_cablingKey, context};
   const LArOnOffIdMapping* cabling=*cablingHdl;
   if(!cabling) {
     ATH_MSG_ERROR("Failed to retrieve LAr Cabling map with key " << m_cablingKey.key() );
     return StatusCode::FAILURE;
   }
 
  SG::ReadCondHandle<LArADC2MeV> adc2mevHdl(m_adc2mevKey, context);
  const LArADC2MeV* adc2MeVs=*adc2mevHdl;
 
  SG::ReadCondHandle<ILArfSampl> fSamplHdl(m_fSamplKey, context);
  const ILArfSampl* fSampl=*fSamplHdl;
 
  SG::ReadCondHandle<ILArPedestal> pedHdl(m_pedestalKey, context);
  const ILArPedestal* pedestal=*pedHdl;
 
  const ILArNoise* noise=nullptr;
  if ( (!m_RndmEvtOverlay || m_isMcOverlay)  && !m_pedestalNoise && m_NoiseOnOff ){
    SG::ReadCondHandle<ILArNoise> noiseHdl(m_noiseKey, context);
    noise=*noiseHdl;
  }
 
  const LArAutoCorrNoise* autoCorrNoise=nullptr;
  if ( !m_RndmEvtOverlay  &&  m_NoiseOnOff) {
    SG::ReadCondHandle<LArAutoCorrNoise>  autoCorrNoiseHdl(m_autoCorrNoiseKey, context);
    autoCorrNoise=*autoCorrNoiseHdl;
  }

  SG::ReadCondHandle<LArBadChannelCont> bch{m_bcContKey,context};
  const LArBadChannelCont* bcCont{*bch};
 
  SG::ReadCondHandle<ILArShape> shapeHdl(m_shapeKey, context);
  const ILArShape* shape=*shapeHdl;
 
   // Inputs
   SG::ReadHandle<LArHitEMap> hitmap(m_hitMapKey,context);
   const LArHitEMap* hitmapPtr = hitmap.cptr();
   const LArHitEMap* hitmapPtr_DigiHSTruth = nullptr;
   if ( m_doDigiTruth ) {
     SG::ReadHandle<LArHitEMap> hitmap_DigitHSTruth(m_hitMapKey_DigiHSTruth,context);
     hitmapPtr_DigiHSTruth = hitmap_DigitHSTruth.cptr();
   }
 
   const size_t nCells=hitmapPtr->GetNbCells();
   // Prepare Output
   //
   // For the standard one lets use a DataPool
   auto DigitContainer = std::make_unique<LArDigitContainer>(SG::VIEW_ELEMENTS);
   DigitContainer->reserve(nCells);
   DataPool<LArDigit> dataItemsPool(context);
   dataItemsPool.reserve(nCells);
   //
   // HSTruth output might not be needed so avoid doing anything
   // in that case
   std::unique_ptr<LArDigitContainer> DigitContainer_DigiHSTruth = nullptr;
   if (m_doDigiTruth){
     DigitContainer_DigiHSTruth = std::make_unique<LArDigitContainer>();
     DigitContainer_DigiHSTruth->reserve(nCells);
   }
   //
   const std::vector<std::pair<float,float> >* TimeE;
   const std::vector<std::pair<float,float> >* TimeE_DigiHSTruth = nullptr;
 
   ATHRNG::RNGWrapper* rngWrapper = m_rndmGenSvc->getEngine(this, m_randomStreamName);
   CLHEP::HepRandomEngine * engine = rngWrapper->getEngine(context);
   ATHRNG::RNGWrapper::SeedingOptionType seedingmode=m_useLegacyRandomSeeds ? ATHRNG::RNGWrapper::MC16Seeding : ATHRNG::RNGWrapper::SeedingDefault;
   rngWrapper->setSeedLegacy( m_randomStreamName, context, m_randomSeedOffset, seedingmode );
 
   for (size_t it=0;it<nCells;++it) 
   {
      const LArHitList& hitlist = hitmapPtr->GetCell(it);
 
      if (!m_Windows || hitlist.inWindows()) {
        TimeE = &(hitlist.getData());
        if(m_doDigiTruth) {
          const auto& hitlist_DigiHSTruth=hitmapPtr_DigiHSTruth->GetCell(it);
          TimeE_DigiHSTruth = &(hitlist_DigiHSTruth.getData());
        }
 
        if (!TimeE->empty() || m_NoiseOnOff || m_RndmEvtOverlay) {
           const Identifier cellID=m_calocell_id->cell_id(IdentifierHash(it));
           HWIdentifier ch_id = cabling->createSignalChannelIDFromHash(IdentifierHash(it));
           HWIdentifier febId = m_laronline_id->feb_Id(ch_id);
           bool missing=!(badFebs->status(febId).good());
           if (!missing) {
             const LArDigit * digit = nullptr ;
             if(m_RndmEvtOverlay) digit = hitmapPtr->GetDigit(it);
             // MakeDigit called if in no overlay mode or
             // if in overlay mode and random digit exists
             if ((!m_RndmEvtOverlay) || (m_RndmEvtOverlay && digit)) {
               LArDigit* Digit = nullptr;
               LArDigit* Digit_DigiHSTruth = nullptr;
               auto sc = MakeDigit(context, cellID, ch_id, Digit,
                                   dataItemsPool,
                                   Digit_DigiHSTruth, TimeE, digit, engine,
                                   TimeE_DigiHSTruth,
                                   adc2MeVs,
                                   fSampl,
                                   pedestal,
                                   noise,
                                   autoCorrNoise,
                                   bcCont,
                                   shape);
 
               if (sc.isFailure()){
                 ATH_MSG_ERROR("LArHitEMapToDigitAlg::execute failed in MakeDigit");
                 delete Digit_DigiHSTruth;
                 return sc;
               }
               DigitContainer->push_back(Digit);
               if (DigitContainer_DigiHSTruth){
                 DigitContainer_DigiHSTruth->push_back(Digit_DigiHSTruth);
               }
             }
           }
        }
      }     // check window
   }        // end of loop over the cells
 
 
  ATH_MSG_DEBUG(" number of created digits  = " << DigitContainer->size());
 
  SG::WriteHandle<LArDigitContainer> DigitContainerHandle( m_DigitContainerName, context);
  ATH_CHECK(DigitContainerHandle.record( std::move(DigitContainer) ) );
  if ( DigitContainer_DigiHSTruth ){
    SG::WriteHandle<LArDigitContainer> DigitContainer_DigiHSTruthHandle( m_DigitContainerName_DigiHSTruth, context);
    ATH_CHECK(DigitContainer_DigiHSTruthHandle.record( std::move(DigitContainer_DigiHSTruth) ) );
  }
 
 
  return StatusCode::SUCCESS;
}
 
StatusCode LArHitEMapToDigitAlg::MakeDigit(
    const EventContext& ctx,
    const Identifier& cellId,
    const HWIdentifier& ch_id,
    LArDigit*& Digit,
    DataPool<LArDigit>& dataItemsPool,
    LArDigit*& Digit_DigiHSTruth,
    const std::vector<std::pair<float, float>>* TimeE,
    const LArDigit* rndmEvtDigit,
    CLHEP::HepRandomEngine* engine,
    const std::vector<std::pair<float, float>>* TimeE_DigiHSTruth,
    const LArADC2MeV* adc2MeVs,
    const ILArfSampl* fSampl,
    const ILArPedestal* pedestal,
    const ILArNoise* noise,
    const LArAutoCorrNoise* autoCorrNoise,
    const LArBadChannelCont* bcCont,
    const ILArShape* shape) const {
 
  bool createDigit_DigiHSTruth = true;
 
 
  int i;
  short Adc;
  short Adc_DigiHSTruth;
 
  std::vector<short> AdcSample(m_NSamples);
  std::vector<short> AdcSample_DigiHSTruth(m_NSamples);
 
  float SF=1.;
  float SigmaNoise;
  staticVecFloat_t rndm_energy_samples(m_NSamples) ;
 
  int iCalo = EM;
  if (m_larem_id->is_lar_hec(cellId))
    iCalo = HEC;
  else if (m_larem_id->is_lar_fcal(cellId))
    iCalo = FCAL;
  else if (m_larem_id->is_em_endcap_inner(cellId))
    iCalo = EMIW;
 
  CaloGain::CaloGain  initialGain=static_cast<CaloGain::CaloGain>(m_gainRange[iCalo].value().first); 
  CaloGain::CaloGain rndmGain = CaloGain::LARHIGHGAIN;
 
// ........ retrieve data (1/2) ................................
//
  SF=fSampl->FSAMPL(ch_id);
 
//
// ....... dump info ................................
//
#ifndef NDEBUG
  ATH_MSG_DEBUG("    Cellid " << m_larem_id->show_to_string(cellId));
  ATH_MSG_DEBUG("    SF: " << SF);
#endif
 
  staticVecDouble_t Samples;
  staticVecDouble_t Samples_DigiHSTruth;
  staticVecDouble_t Noise;
  Samples.resize(m_NSamples,0);
  if(m_doDigiTruth) Samples_DigiHSTruth.resize(m_NSamples,0);
  Noise.resize(m_NSamples,0);
 
//
// ....... make the five samples
//
 
#ifndef NDEBUG
  ATH_MSG_DEBUG(" number of hit for this cell " << TimeE->size());
#endif
 
//
// convert Hits into energy samples and add result to Samples assuming LARHIGHGAIN for pulse shape
//
  bool isDead = m_bcMask.cellShouldBeMasked(bcCont,ch_id);
 
 
  if (!isDead) {
    if( this->ConvertHits2Samples(cellId,ch_id,initialGain,TimeE, Samples, shape).isFailure() ) {
      return StatusCode::SUCCESS;
    }
    if(m_doDigiTruth && TimeE_DigiHSTruth){
      if( this->ConvertHits2Samples(cellId,ch_id,initialGain,TimeE_DigiHSTruth, Samples_DigiHSTruth, shape).isFailure() ) {
        return StatusCode::SUCCESS;
      }
    }
  }
 
//
// .... add random event digit if needed
//
 float energy2adc ;
 float rAdc ;
 if(m_RndmEvtOverlay && rndmEvtDigit ) // no overlay if missing random digit
 {
  rndmGain= rndmEvtDigit->gain();
  auto polynom_adc2mev =adc2MeVs->ADC2MEV(ch_id,rndmEvtDigit->gain());
  if (polynom_adc2mev.size() > 1) {
     float adc2energy = SF * polynom_adc2mev[1];
     const std::vector<short> & rndm_digit_samples = rndmEvtDigit->samples() ;
     float Pedestal = pedestal->pedestal(ch_id,rndmEvtDigit->gain());
     if (Pedestal <= (1.0+LArElecCalib::ERRORCODE)) {
       ATH_MSG_WARNING("  Pedestal not found in database for this channel offID " << cellId << " Use sample 0 for random");
       Pedestal = rndm_digit_samples[0];
     }
     ATH_MSG_DEBUG(" Params for inverting LAr Digitization: pedestal " << Pedestal << "  adc2energy " << adc2energy);
 
// in case Medium or low gain, take into account ramp intercept in ADC->"energy" computation
//   this requires to take into account the sum of the optimal filter coefficients, as they don't compute with ADC shift
     float adc0=0.;
     if (!m_isMcOverlay && rndmEvtDigit->gain()>0) {
        SG::ReadCondHandle<ILArOFC> larOFC(m_OFCKey, ctx);
        if (larOFC.cptr() != nullptr) {
          ILArOFC::OFCRef_t ofc_a = larOFC->OFC_a(ch_id,rndmEvtDigit->gain(),0);
          float sumOfc=0.;
          if (ofc_a.size()>0) {
            for (unsigned int j=0;j<ofc_a.size();j++) sumOfc += ofc_a.at(j);
          }
          if (sumOfc>0) adc0 =  polynom_adc2mev[0] * SF /sumOfc;
        }
     }
 
     int nmax=m_NSamples;
     if ((int)(rndm_digit_samples.size()) < m_NSamples) {
        ATH_MSG_WARNING(
            "Less digit Samples than requested in digitization for cell "
            << ch_id.get_compact() << " Digit has " << rndm_digit_samples.size()
            << " samples.  Digitization request " << m_NSamples);
        nmax = rndm_digit_samples.size();
     }
     for(i=0 ; i<nmax ; i++)
     {
       rAdc = (rndm_digit_samples[i] - Pedestal ) * adc2energy + adc0;
       rndm_energy_samples[i] = rAdc ;
       Samples[i] += rAdc ;
     }
  }
  else {
    ATH_MSG_ERROR(" No ramp found for this random cell " << m_larem_id->show_to_string(cellId) << " for gain " << rndmEvtDigit->gain());
    if (adc2MeVs->nGains()<=(int)rndmEvtDigit->gain()) {
      ATH_MSG_ERROR("Found ramp for only " << adc2MeVs->nGains() << " gains. Are you trying to overlay of Run 3 random with run 4 signal?");
    }
    return  StatusCode::FAILURE;
  }
 }
 
 
  CaloGain::CaloGain igain=chooseGain(Samples,ch_id,static_cast<CaloNum>(iCalo),pedestal,adc2MeVs,SF);
  if (igain==CaloGain::INVALIDGAIN) {
        return StatusCode::FAILURE;
  }
 
 // check that select gain is never lower (higher index number) than random gain in case of overlay
  igain=std::max(rndmGain,igain);
 
//
// recompute Samples if igain != HIGHGAIN
//
   if (igain != initialGain ){
     for (i=0;i<m_NSamples;i++) {
       if(m_doDigiTruth) Samples_DigiHSTruth[i] = 0.;
       if (m_RndmEvtOverlay) Samples[i]= rndm_energy_samples[i] ;
       else Samples[i] = 0.;
     }
 
     if (!isDead) {
       if( this->ConvertHits2Samples(cellId,ch_id,igain,TimeE, Samples, shape) == StatusCode::FAILURE ) {
         return StatusCode::SUCCESS;
       }
       if(m_doDigiTruth){
         if( this->ConvertHits2Samples(cellId,ch_id,igain,TimeE_DigiHSTruth, Samples_DigiHSTruth, shape) == StatusCode::FAILURE ) {
           return StatusCode::SUCCESS;
         }
       }
     }
   }
 
//
// ........ add the noise ................................
//
 
  double Rndm[32]{};
  int BvsEC=0;
  if(iCalo==EM || iCalo==EMIW) BvsEC=std::abs(m_larem_id->barrel_ec(cellId));
 
  bool addedNoise=false;
  if (m_NoiseOnOff &&
      ((BvsEC == 1 && m_NoiseInEMB) || (BvsEC > 1 && m_NoiseInEMEC) ||
       (iCalo == HEC && m_NoiseInHEC) || (iCalo == FCAL && m_NoiseInFCAL)))
  // add the noise only in the wanted sub-detectors
  {
     if (!m_RndmEvtOverlay) {
       if (!m_pedestalNoise) {
         SigmaNoise = noise->noise(ch_id, igain);
       } else {
         float thisNoise = pedestal->pedestalRMS(ch_id, igain);
         if (thisNoise >= (1.0 + LArElecCalib::ERRORCODE))
          SigmaNoise = thisNoise;
         else
          SigmaNoise = 0.;
       }
       // Sqrt of noise covariance matrix
       const std::vector<float>& CorGen =
           autoCorrNoise->autoCorrSqrt(cellId, igain);
       if (CorGen.size() < (unsigned)m_NSamples * m_NSamples) {
         ATH_MSG_ERROR("Noise AutoCorr too small, need "
                       << m_NSamples * m_NSamples << " points for "
                       << m_NSamples << " samples.");
         return StatusCode::FAILURE;
       }
 
       RandGaussZiggurat::shootArray(engine, m_NSamples, Rndm, 0., 1.);
 
       int index;
       for (int i = 0; i < m_NSamples; i++) {
         Noise[i] = 0.;
         for (int j = 0; j <= i; j++) {
          index = i * m_NSamples + j;
          Noise[i] += Rndm[j] * CorGen[index];
         }
         Noise[i] = Noise[i] * SigmaNoise;
       }
       addedNoise = true;
     } else {
       // overlay case a priori don't add any noise
       for (int i = 0; i < m_NSamples; i++)
         Noise[i] = 0.;
       // if gain from zerobias events is < gain from mixed events => add extra
       // noise to account for gain vs noise dependance
       //   done in a simple way without taking into account the time
       //   correlation of this extra noise properly
       if (rndmEvtDigit) {
         // if gain of cell is different from ZB event gain
         if (igain > rndmEvtDigit->gain()) {
          double SigmaNoiseZB = 0.;  // noise in ZB event for gain of ZB event
          double SigmaNoise = 0.;    // noise expected for new gain value
          double SigmaExtraNoise = 0.;  // quadratic difference of noise values
          if (!m_pedestalNoise) {
            SigmaNoiseZB = noise->noise(ch_id, rndmEvtDigit->gain());
            SigmaNoise = noise->noise(ch_id, igain);
          } else {
            float thisNoise = pedestal->pedestalRMS(ch_id, rndmEvtDigit->gain());
            if (thisNoise >= (1.0 + LArElecCalib::ERRORCODE))
                 SigmaNoiseZB = thisNoise;
            else
                 SigmaNoiseZB = 0.;
            thisNoise = pedestal->pedestalRMS(ch_id, igain);
            if (thisNoise >= (1.0 + LArElecCalib::ERRORCODE))
                 SigmaNoise = thisNoise;
            else
                 SigmaNoise = 0.;
          }
          // Convert SigmaNoiseZB in noise in ADC counts for igain conversion
          auto polynom_adc2mevZB =
              adc2MeVs->ADC2MEV(cellId, rndmEvtDigit->gain());
          auto polynom_adc2mev = adc2MeVs->ADC2MEV(cellId, igain);
          if (polynom_adc2mevZB.size() > 1 && polynom_adc2mev.size() > 1) {
            if (polynom_adc2mev[1] > 0.) {
                 SigmaNoiseZB = SigmaNoiseZB * (polynom_adc2mevZB[1]) /
                                (polynom_adc2mev[1]);
                 if (SigmaNoise > SigmaNoiseZB)
                   SigmaExtraNoise = sqrt(SigmaNoise * SigmaNoise -
                                          SigmaNoiseZB * SigmaNoiseZB);
            }
          }  // check that AC2MeV factors are there
          RandGaussZiggurat::shootArray(engine, m_NSamples, Rndm, 0.,
                                        1.);  // generate noise
          for (int i = 0; i < m_NSamples; i++)
            Noise[i] = SigmaExtraNoise * Rndm[i];
          addedNoise = true;
         }  // different gains
       }    // rndm Digit is there
     }      // rndm Overlay test
  }         // add noise ?
            //
// ......... convert into adc counts  ................................
//
  float Pedestal = pedestal->pedestal(ch_id,igain);
  if (Pedestal <= (1.0+LArElecCalib::ERRORCODE)) {
     ATH_MSG_WARNING(" pedestal not found for cellId " << cellId << " assume 1000" );
     Pedestal=1000.;
  }
  const auto polynom_adc2mev = adc2MeVs->ADC2MEV(cellId,igain);
  if (polynom_adc2mev.size() < 2) {
    ATH_MSG_WARNING(" No ramp found for requested gain " << igain << " for cell " << m_larem_id->show_to_string(cellId) << " no digit made...");
    return StatusCode::SUCCESS;
  }
 
  energy2adc=1./(polynom_adc2mev[1])/SF;
 
// in case Medium or low gain, take into account ramp intercept in energy->ADC computation
//   this requires to take into account the sum of the optimal filter coefficients, as they don't compute with ADC shift
  if(!m_isMcOverlay && m_RndmEvtOverlay  && igain>0)
  {
    SG::ReadCondHandle<ILArOFC> larOFC(m_OFCKey, ctx);
    if (larOFC.cptr() != nullptr) {
      float sumOfc=0.;
      ILArOFC::OFCRef_t ofc_a = larOFC->OFC_a(ch_id,igain,0);
      if (ofc_a.size()>0) {
        for (unsigned int j=0;j<ofc_a.size();j++) sumOfc+= ofc_a.at(j);
      }
      if ((polynom_adc2mev[1])>0 && sumOfc>0) Pedestal = Pedestal - (polynom_adc2mev[0])/(polynom_adc2mev[1])/sumOfc;
      ATH_MSG_DEBUG("  Params for final LAr Digitization  gain: " << igain << "    pedestal: " << Pedestal <<  "   energy2adc: " << energy2adc);
    }
  }
  for(i=0;i<m_NSamples;i++)
  {
    double xAdc;
    double xAdc_DigiHSTruth = 0;
 
    if ( addedNoise ){
      xAdc =  Samples[i]*energy2adc + Noise[i] + Pedestal + 0.5;
      if(m_doDigiTruth) {
        xAdc_DigiHSTruth =  Samples_DigiHSTruth[i]*energy2adc + Noise[i] + Pedestal + 0.5;
      }
    }
 
    else {
      if (m_roundingNoNoise) {
        float flatRndm = RandFlat::shoot(engine);
        xAdc =  Samples[i]*energy2adc + Pedestal + flatRndm;
        if(m_doDigiTruth) {
          xAdc_DigiHSTruth =  Samples_DigiHSTruth[i]*energy2adc + Pedestal + flatRndm;
        }
 
      }
      else{
         xAdc =  Samples[i]*energy2adc + Pedestal + 0.5;
         if(m_doDigiTruth) {
           xAdc_DigiHSTruth =  Samples_DigiHSTruth[i]*energy2adc + Pedestal + 0.5;
         }
      }
 
    }
 
//
// ........ truncate at maximum value + 1
//          add possibility to saturate at 0 for negative signals
//
    if (xAdc <0) {
        ATH_MSG_WARNING("Found below-zero digit " << xAdc);
        Adc=0;
    }
    else if (xAdc >= m_maxADC) {
      ATH_MSG_WARNING("Found saturating digit " << i << ": " << xAdc << " chid="<< ch_id.get_identifier32().get_compact() << ", gain=" << igain << 
                      ", energy2adc=" << energy2adc << ", E[0]= " << (TimeE->size()>0 ? std::format("{}",TimeE->at(0).first) : "No value"));
        Adc=m_maxADC;
        }
    else 
       Adc = (short) xAdc;
 
    AdcSample[i]=Adc;
 
    if(m_doDigiTruth){
      if (xAdc_DigiHSTruth <0)  Adc_DigiHSTruth=0;
      else if (xAdc_DigiHSTruth >= m_maxADC) Adc_DigiHSTruth=m_maxADC;
      else Adc_DigiHSTruth = (short) xAdc_DigiHSTruth;
      AdcSample_DigiHSTruth[i] = Adc_DigiHSTruth;
    }
 
#ifndef NDEBUG
    ATH_MSG_DEBUG(" Sample " << i << "  Energy= " << Samples[i] << "  Adc=" << Adc);
#endif
 
  }
 
//
// ...... create the LArDigit .............
//
  Digit = dataItemsPool.nextElementPtr();
  (*Digit)=LArDigit(ch_id,igain,std::move(AdcSample));
 
  if (m_doDigiTruth && createDigit_DigiHSTruth) {
    createDigit_DigiHSTruth = false;
    Digit_DigiHSTruth = nullptr;
 
    for (int i = 0; i < m_NSamples; i++) {
      if (Samples_DigiHSTruth[i] != 0)
         createDigit_DigiHSTruth = true;
    }
 
    Digit_DigiHSTruth =
        new LArDigit(ch_id, igain, std::move(AdcSample_DigiHSTruth));
  }
 
  return StatusCode::SUCCESS;
}
 
// ---------------------------------------------------------------------------------------
 
StatusCode LArHitEMapToDigitAlg::ConvertHits2Samples(const Identifier & cellId, const HWIdentifier ch_id, CaloGain::CaloGain igain,
                                                     const std::vector<std::pair<float,float> >  *TimeE, staticVecDouble_t &sampleList,
                                                     const ILArShape* shape) const
 
{
// Converts  hits of a particular LAr cell into energy samples
// declarations
   int nsamples ;
   int nsamples_der ;
   int i ;
   int j ;
 
// ........ retrieve data (1/2) ................................
//
   ILArShape::ShapeRef_t Shape = shape->Shape(ch_id,igain);
   ILArShape::ShapeRef_t ShapeDer = shape->ShapeDer(ch_id,igain);
 
  nsamples = Shape.size();
  nsamples_der = ShapeDer.size();
 
  if (nsamples==0) {
    ATH_MSG_INFO(" No samples for cell = " << cellId );
    return StatusCode::FAILURE;
  }
 
#ifndef NDEBUG
  ATH_MSG_DEBUG(" Cellid " << m_larem_id->show_to_string(cellId));
  for (i=0;i<nsamples;i++)
  {
    msg() << MSG::DEBUG << Shape[i] << " ";
  }
  msg() << endmsg;
  ATH_MSG_DEBUG("m_NSamples, m_usePhase " << m_NSamples << " " << m_usePhase);
#endif
 
 
for (const auto& [energy, time] : *TimeE) {
  // fix the shift +1 if HEC  and nSamples 4 and firstSample 0
   // in case of data overlay this should NOT  be done as the pulse shape read from the database is already shifted
   //   but this should still be done in case of MC overlay
   int ihecshift=0;
   if((!m_RndmEvtOverlay || m_isMcOverlay) && m_larem_id->is_lar_hec(cellId) && m_NSamples.value() == 4 && m_firstSample.value() == 0) ihecshift=1;
 
 
   if (!m_usePhase) {
 
 // Atlas like mode where we use 25ns binned pulse shape and derivative to deal with time offsets
 
// shift between reference shape and this time
      int ishift=(int)(rint(time*(1./25.)));
      double dtime=time-25.*((double)(ishift));
      for (i=0;i<m_NSamples.value();i++)
      {
       j = i - ishift + m_firstSample + ihecshift;
#ifndef NDEBUG
       msg() << MSG::DEBUG << " time/i/j " << time << " "<< i << " " << j;
#endif
       if (j >=0 && j < nsamples ) {
         if (j<nsamples_der && std::abs(ShapeDer[j])<10. )
              sampleList[i] += (Shape[j]- ShapeDer[j]*dtime)*energy ;
         else sampleList[i] += Shape[j]*energy ;
       }
      }
      msg() << MSG::DEBUG << endmsg;
   }
// Mode to use phase (tbin) to get pulse shape ( pulse shape with fine time binning should be available)
 
   else {
 
     // FIXME hardcode 8phases3ns configuration (cannot access parameters from ILArShape interface now)
     int nTimeBins = 8;
     float timeBinWidth = 25./24.*3.;
 
//    -50<t<-25 phase=-t-25, shift by one peak time (for s2 uses shape(3) with tbin)
// for -25<t<0  phase = -t, no shift of peak time
// for 0<t<25   phase=25-t, shift by one peak time (for s2 uses shape(1) with tbin)
//    25<t<50   phase=50-t, shift by two
//  etc...
 
      int ishift = (int)(time*(1./25.));
      int tbin;
      if (time>0) {
         tbin = (int)(fmod(time,25)/timeBinWidth);
          if (tbin>0) {
             tbin=nTimeBins-tbin;
             ishift +=1;
          }
      } else {
         tbin = (int)(fmod(-time,25)/timeBinWidth);
      }
 
      double dtime = time - ( 25.*((float)(ishift)) - timeBinWidth*tbin);
 
      Shape = shape->Shape(ch_id,igain,tbin);
      ShapeDer = shape->ShapeDer(ch_id,igain,tbin);
 
      nsamples = Shape.size();
      nsamples_der = ShapeDer.size();
 
 
      for (i=0;i<m_NSamples.value();i++)
      {
       j = i - ishift+m_firstSample + ihecshift;
#ifndef NDEBUG
       ATH_MSG_DEBUG(" time/i/j " << time << " "<< i << " " << j);
#endif
       if (j >=0 && j < nsamples ) {
         if (j<nsamples_der && std::abs(ShapeDer[j])<10. )
              sampleList[i] += (Shape[j]- ShapeDer[j]*dtime)*energy ;
         else sampleList[i] += Shape[j]*energy ;
       }
      }
 
   }     // else if of m_usePhase
  } // loop over hits
 
   return StatusCode::SUCCESS;
 
}
 
CaloGain::CaloGain LArHitEMapToDigitAlg::chooseGain(const staticVecDouble_t& samples, const HWIdentifier ch_id, const CaloNum iCalo,
                                                    const ILArPedestal* pedestal, const LArADC2MeV* adc2MeVs, const float SF) const {
 
 
  const CaloGain::CaloGain gainChoosingGain=CaloGain::LARMEDIUMGAIN;                                                       
  int sampleGainChoice{2};
  if (m_firstSample < 0)
    sampleGainChoice -= m_firstSample;
 
  // fix the shift +1 if HEC  and nSamples 4 and firstSample 0
  if (iCalo == HEC && m_NSamples.value() == 4 && m_firstSample.value() == 0)
    sampleGainChoice -= 1;  // ihecshift
 
  // We choose the gain in applying thresholds on the 3rd Sample (index "2")
  // converted in ADC counts in the second-lowest gain (eg MEDIUM gain for run 1,2,3, HIGH gain for run 4)
  // Indeed, thresholds in ADC counts are defined with respect to the MediumGain.
  //
  //               1300              3900
  //  ---------------|----------------|--------------> ADC counts in MediumGain
  //     HighGain  <---  MediumGain  --->  LowGain
 
  float Pedestal = pedestal->pedestal(ch_id, gainChoosingGain);
  if (Pedestal <= (1.0 + LArElecCalib::ERRORCODE)) {
    ATH_MSG_DEBUG(" Pedestal not found for channel " << m_laronline_id->channel_name(ch_id) << " assume 1000 ");
    Pedestal = 1000.;
  }
  const auto& polynom_adc2mev = adc2MeVs->ADC2MEV(ch_id, gainChoosingGain);
  if (polynom_adc2mev.size() < 2) {
    ATH_MSG_WARNING(" No ramp found for channel  " << m_laronline_id->channel_name(ch_id) << ", gain " << gainChoosingGain << ",  gain choice failed...");
    return CaloGain::INVALIDGAIN;
  }
  const float pseudoADC3 = samples[sampleGainChoice] / (polynom_adc2mev[1]) / SF + Pedestal;
 
  CaloGain::CaloGain igain = gainChoosingGain;
  // if we are not yet already in the highest gain and we are below high-gain theshold, switch to high gain
  if (gainChoosingGain > m_gainRange[iCalo].value().first && pseudoADC3 < m_HighGainThresh[iCalo]) {
    igain=CaloGain::LARHIGHGAIN;
  }
 
  else if (gainChoosingGain < m_gainRange[iCalo].value().second && pseudoADC3 > m_LowGainThresh[iCalo]) {
    igain = CaloGain::LARLOWGAIN;
  }
 
  return igain;
}