LArShapeDumper.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "LArCafJobs/LArShapeDumper.h"
#include "StoreGate/StoreGateSvc.h"
#include "LArRawEvent/LArDigit.h"
#include "GaudiKernel/INTupleSvc.h"
#include "LArRawEvent/LArOFIterResultsContainer.h"
#include "LArRawEvent/LArFebErrorSummary.h"
#include "LArElecCalib/ILArShape.h"
 
#include "LArRawConditions/LArPhysWave.h"
#include "LArRawConditions/LArShapeComplete.h"
#include "LArRawConditions/LArAutoCorrComplete.h"
 
#include "CaloDetDescr/CaloDetDescriptor.h"
#include "CaloDetDescr/CaloDetDescrElement.h"
#include "xAODEventInfo/EventInfo.h"
#include "LArRecConditions/LArBadChannel.h"
#include "LArCafJobs/DataContainer.h"
#include "LArCafJobs/ShapeInfo.h"
#include "LArCafJobs/CellInfo.h"
#include "LArCafJobs/RunData.h"
#include "LArCafJobs/EventData.h"
#include "TrigT1Result/RoIBResult.h"
#include "TrigSteeringEvent/TrigRoiDescriptor.h"
#include "TFile.h"
#include "TMath.h"
 
#include <regex>
 
 
#include <vector>
#include <iostream>
using std::cout;
using std::endl;
 
using namespace LArSamples;
 
 
LArShapeDumper::LArShapeDumper(const std::string & name, ISvcLocator * pSvcLocator) : 
  AthAlgorithm(name, pSvcLocator),
  m_count(0),
  m_nWrongBunchGroup(0),
  m_nPrescaledAway(0),
  m_nLArError(0),
  m_nNoDigits(0),
  m_nNoDigitsSC(0),
  m_onlineHelper(nullptr),
  m_onlineHelperSC(nullptr),
  m_doEM(false),
  m_doHEC(false),
  m_doFCAL(false),
  m_doSC(false),
  m_samples(nullptr)
{
  declareProperty("FileName", m_fileName = "samples.root");
  declareProperty("MaxChannels", m_maxChannels = 200000);
  declareProperty("Prescale", m_prescale = 1);
  declareProperty("CaloType", m_caloType = "EMHECFCALSC");
  declareProperty("EnergyCut", m_energyCut = -1);
  declareProperty("EnergyCutSC", m_energyCutSC = -1);
  declareProperty("NoiseSignifCut", m_noiseSignifCut = 3);
  declareProperty("MinADCMax", m_minADCMax = -1);
  declareProperty("MinADCMaxSC", m_minADCMaxSC = -1);
  declareProperty("Gains", m_gainSpec = "HIGH,MEDIUM,LOW");
  declareProperty("DumpDisconnected", m_dumpDisc = false);
  declareProperty("DoStream", m_doStream = false);
  declareProperty("DoTrigger", m_doTrigger = true);
  declareProperty("DoOFCIter", m_doOFCIter = true);
  declareProperty("DoAllEvents", m_doAllEvents = true);
  declareProperty("DumpChannelInfos", m_dumpChannelInfos = false);
  declareProperty("DoRoIs", m_doRoIs = true);
  declareProperty("TriggerNames", m_triggerNames);
  declareProperty("DoAllLvl1", m_doAllLvl1 = true);
  declareProperty("onlyEmptyBC",m_onlyEmptyBC=false);
}
 
 
LArShapeDumper::~LArShapeDumper() 
{
}
 
 
StatusCode LArShapeDumper::initialize()
{
  ATH_MSG_DEBUG ("in initialize()");
 
  m_samples = new DataStore();
 
  std::transform(m_caloType.begin(), m_caloType.end(), m_caloType.begin(), toupper);
  m_doEM   = (m_caloType.find("EM")   != std::string::npos);
  m_doHEC  = (m_caloType.find("HEC")  != std::string::npos);
  m_doFCAL = (m_caloType.find("FCAL") != std::string::npos);
  m_doSC   = (m_caloType.find("SC") != std::string::npos);
 
  ATH_CHECK( m_cablingKey.initialize() );
  ATH_CHECK( m_BCKey.initialize() );
  ATH_CHECK( m_BCKeySC.initialize(m_doSC) );
  ATH_CHECK( m_noiseCDOKey.initialize() );
  ATH_CHECK( m_adc2mevKey.initialize() );
  ATH_CHECK( m_pedestalKey.initialize() );
  ATH_CHECK( m_bcDataKey.initialize() );
 
  ATH_CHECK( m_digitsKey.initialize() );
  ATH_CHECK( m_channelsKey.initialize() );
 
  ATH_CHECK( m_digitsKeySC.initialize(m_doSC) );
  ATH_CHECK( m_rawscKey.initialize(m_doSC) );
  ATH_CHECK( m_cablingKeySC.initialize(m_doSC) );
  ATH_CHECK( m_pedestalKeySC.initialize(m_doSC) );
 
  ATH_CHECK( detStore()->retrieve(m_onlineHelper, "LArOnlineID") );
  ATH_CHECK(m_caloMgrKey.initialize());
  ATH_CHECK(m_caloSuperCellMgrKey.initialize(m_doSC));
  if(m_doSC) ATH_CHECK( detStore()->retrieve(m_onlineHelperSC, "LArOnline_SuperCellID") );
 
  ATH_CHECK(m_bcMask.buildBitMask(m_problemsToMask,msg()));
  if(m_doSC) ATH_CHECK(m_bcMaskSC.buildBitMask(m_problemsToMaskSC,msg()));
 
 
  if (m_doTrigger) {
    CHECK( m_trigDec.retrieve() );
  }
 
  ATH_CHECK( m_dumperTool.retrieve() );
 
  ATH_CHECK(m_acorrKey.initialize(m_dumperTool->doShape()) );
  
  std::transform(m_gainSpec.begin(), m_gainSpec.end(), m_gainSpec.begin(), toupper);
  m_gains[CaloGain::LARHIGHGAIN]   = (m_gainSpec.find("HIGH")   != std::string::npos);
  m_gains[CaloGain::LARMEDIUMGAIN] = (m_gainSpec.find("MEDIUM") != std::string::npos);
  m_gains[CaloGain::LARLOWGAIN]    = (m_gainSpec.find("LOW")    != std::string::npos);
  
  //if (m_onlyEmptyBC)
  //ATH_CHECK(m_bcidTool.retrieve());
 
  return StatusCode::SUCCESS; 
}
 
 
StatusCode LArShapeDumper::start()
{
  m_runData = std::make_unique<RunData>(0);
 
  if (m_doTrigger) {
    std::vector<std::regex> regexs;
    regexs.reserve(m_triggerNames.size());
for (const std::string& name : m_triggerNames) {
      regexs.push_back(std::regex(name));
    }
 
    std::vector<std::string> chains = m_trigDec->getListOfTriggers();
    std::vector<std::string> myChains;
    std::cmatch match;
 
    for (const std::string& chain : chains) {
      ATH_MSG_INFO ( "Configured chain : " << chain );
      for (const std::regex& regex : regexs) {
        if (std::regex_match(chain.c_str(), match, regex)) myChains.push_back(chain);
      }
    }
    for (const std::string& group : m_trigDec->getListOfGroups())
      ATH_MSG_INFO ( "Configured group : " << group );
    const Trig::ChainGroup* calibStreamGroup = m_trigDec->getChainGroup("Calibration"); 
    if (calibStreamGroup) {
      std::vector<std::string> chains = calibStreamGroup->getListOfTriggers();
      ATH_MSG_INFO ( "Chains for Calibration group:" );
      for (const std::string& chain : chains)
        ATH_MSG_INFO ( "Calib chain : " << chain );
    }
 
    unsigned int idx = 0;
  
    if (m_doAllLvl1) {
      const Trig::ChainGroup* group = m_trigDec->getChainGroup("L1_.*");
      for (const std::string& l1Item : group->getListOfTriggers()) {
        const TrigConf::TriggerItem* confItem = m_trigDec->ExperimentalAndExpertMethods().getItemConfigurationDetails(l1Item);
        if (!confItem) {
          ATH_MSG_WARNING ( "LVL1 item " << l1Item << ", obtained from TrigConfig, cannot be retrieved!" );
          continue;
        }
        int pos = confItem->ctpId();
        if (pos < 0 || pos >= 256) {
          ATH_MSG_WARNING ( "LVL1 item " << l1Item << "has out-of-range ctpId " << pos );
          continue;
        }
        m_runData->addBit(l1Item.c_str(), pos);
        ATH_MSG_INFO ( "Adding LVL1 trigger bit for " << l1Item << " at position " << pos );
      }
      idx = 256;
    }
 
    for (const std::string& name : myChains) {
      if (m_trigDec->getListOfTriggers(name).empty()) {
        ATH_MSG_WARNING ( "Requested trigger name " << name << " is not configured in this run" );
        continue;
      }
      const Trig::ChainGroup* group = m_trigDec->getChainGroup(name);
      if (!group) {
        ATH_MSG_WARNING ( "Could not retrieve chain group for trigger " << name );
        continue;
      }
      m_runData->addBit(name.c_str(), idx++);
      m_triggerGroups.push_back(group);
      ATH_MSG_INFO ( "Adding trigger bit for " << name << " at position " << idx-1 );
    }
  }
  return StatusCode::SUCCESS;
}
 
 
StatusCode LArShapeDumper::execute()
{    
  m_count++;
  const EventContext& ctx = Gaudi::Hive::currentContext();
  if ((m_prescale > 1 && m_random.Rndm() > 1.0/m_prescale) || m_prescale <= 0) {
    ATH_MSG_VERBOSE ( "======== prescaling event "<< m_count << " ========" );
    m_nPrescaledAway++;
    return StatusCode::SUCCESS;
  }
 
  ATH_MSG_VERBOSE ( "======== executing event "<< m_count << " ========" );
 
  const xAOD::EventInfo* eventInfo = nullptr;
  ATH_CHECK( evtStore()->retrieve(eventInfo) );
  
  int event     = eventInfo->eventNumber();
  int run       = eventInfo->runNumber();
  int lumiBlock = eventInfo->lumiBlock();
  int bunchId   = eventInfo->bcid();
 
  
  if (eventInfo->errorState(xAOD::EventInfo::LAr)==xAOD::EventInfo::Error) {
    ATH_MSG_DEBUG("Ignoring Event b/c of LAr ERROR");
    m_nLArError++;
    return StatusCode::SUCCESS;
  }
 
 
  SG::ReadCondHandle<BunchCrossingCondData> bccd (m_bcDataKey,ctx);
  const BunchCrossingCondData* bunchCrossing=*bccd;
  if (!bunchCrossing) {
    ATH_MSG_ERROR("Failed to retrieve Bunch Crossing obj");
    return StatusCode::FAILURE;
  }
 
  SG::ReadCondHandle<CaloDetDescrManager> caloMgrHandle{m_caloMgrKey};
  ATH_CHECK(caloMgrHandle.isValid());
  const CaloDetDescrManager* caloMgr = *caloMgrHandle;  
 
  if (m_onlyEmptyBC) {
    if (!bccd->isFilled(bunchId)) {
      ATH_MSG_DEBUG("Ignoring Event with bunch crossing type ");
      m_nWrongBunchGroup++;
      return StatusCode::SUCCESS;
    }
   }
 
 
  EventData* eventData = nullptr;
  int eventIndex = -1;
  if (m_doAllEvents) {
    eventIndex = makeEvent(eventData, run, event, lumiBlock, bunchId);
    if (eventIndex < 0) return StatusCode::FAILURE;
  }
 
  SG::ReadHandle<LArDigitContainer> hdlDigit(m_digitsKey, ctx);
  if(!hdlDigit.isValid()) {
     ATH_MSG_WARNING( "Unable to retrieve LArDigitContainer with key " << m_digitsKey << " from DetectorStore. " );
     return StatusCode::SUCCESS;
  } else
      ATH_MSG_DEBUG( "Got LArDigitContainer with key " << m_digitsKey.key() );
  const LArDigitContainer* larDigitContainer = &(*hdlDigit);
 
  if (larDigitContainer->empty()) {
    ATH_MSG_WARNING ( "LArDigitContainer with key=" << m_digitsKey.key() << " is empty!" );
    m_nNoDigits++;
    return StatusCode::SUCCESS;
  }
 
  SG::ReadHandle<LArRawChannelContainer> hdlRaw(m_channelsKey, ctx);
  if(!hdlRaw.isValid()) {
     ATH_MSG_WARNING( "Unable to retrieve LArRawChannelsContainer with key " << m_channelsKey.key() << " from DetectorStore. " );
     return StatusCode::SUCCESS;
  } else
      ATH_MSG_DEBUG( "Got LArRawChannelsContainer with key " << m_channelsKey.key() );
  const LArRawChannelContainer* rawChannelContainer = &(*hdlRaw);;
 
  if (rawChannelContainer->empty()) {
    ATH_MSG_WARNING ( "LArRawChannelContainer with key=" << m_channelsKey << " is empty!" );
    m_nNoDigits++;
    return StatusCode::SUCCESS;
  }
 
  
  
  SG::ReadCondHandle<LArOnOffIdMapping> cablingHdl{m_cablingKey};
  const LArOnOffIdMapping* cabling=*cablingHdl;
  if(!cabling) {
     ATH_MSG_ERROR( "Do not have cabling object LArOnOffIdMapping" );
     return StatusCode::FAILURE;
  }
 
  SG::ReadCondHandle<LArADC2MeV> adc2mevHdl(m_adc2mevKey, ctx);
  const LArADC2MeV* adc2MeV=*adc2mevHdl;
  if(!adc2MeV) {
     ATH_MSG_ERROR( "Failed to retreive ADC2MeV cond obj" );
     return StatusCode::FAILURE;
  }
 
  SG::ReadCondHandle<ILArPedestal> pedHdl(m_pedestalKey, ctx);
  const ILArPedestal* pedestals=*pedHdl;
  if (!pedestals) {
    ATH_MSG_ERROR("Failed to retrieve pedestal cond obj");
     return StatusCode::FAILURE;
  }
 
  const ILArAutoCorr* aCorr=nullptr;
  if(m_dumperTool->doShape()) {
     SG::ReadCondHandle<ILArAutoCorr> acorrHdl(m_acorrKey, ctx);
     aCorr=*acorrHdl;
     if (!aCorr) {
        ATH_MSG_ERROR("Failed to retrieve AutoCorr cond obj");
        return StatusCode::FAILURE;
     }
  }
 
  SG::ReadCondHandle<LArBadChannelCont> readHandle{m_BCKey};
  const LArBadChannelCont *bcCont {*readHandle};
  if(!bcCont) {
     ATH_MSG_ERROR( "Do not have Bad chan container " << m_BCKey.key() );
     return StatusCode::FAILURE;
  }
  
  const LArOFIterResultsContainer* ofIterResult = nullptr;
  if (m_doOFCIter) {
    if (evtStore()->contains<LArOFIterResultsContainer> ("LArOFIterResult")) {
       ATH_CHECK( evtStore()->retrieve(ofIterResult, "LArOFIterResult") );
    } else {
       ATH_MSG_WARNING("Do not have LArOFIterResult in this event");
    }
  }
 
  const LArFebErrorSummary* larFebErrorSummary = nullptr;
  ATH_CHECK( evtStore()->retrieve(larFebErrorSummary, "LArFebErrorSummary") );
  const std::map<unsigned int,uint16_t>& febErrorMap = larFebErrorSummary->get_all_febs();
  std::map<unsigned int, const LArRawChannel*> channelsToKeep;
  
  for (LArRawChannelContainer::const_iterator channel = rawChannelContainer->begin(); 
       channel != rawChannelContainer->end(); ++channel) 
  {
    if (m_energyCut > 0 && TMath::Abs(channel->energy()) < m_energyCut) continue;
    if (m_bcMask.cellShouldBeMasked(bcCont,channel->channelID())) continue;
 
    IdentifierHash hash = m_onlineHelper->channel_Hash(channel->channelID());
    
    if (!hash.is_valid()) {
      ATH_MSG_FATAL ( "Found a LArRawChannel whose HWIdentifier (" << channel->channelID()
                      << ") does not correspond to a valid hash -- returning StatusCode::FAILURE." );
      return StatusCode::FAILURE;
    }    
    channelsToKeep[hash] = &*channel;
    if (m_dumpChannelInfos) {
      HistoryContainer* histCont = m_samples->hist_cont(hash);
      CellInfo* info = nullptr;
      if (!histCont) {
        HWIdentifier channelID = channel->hardwareID();
        const Identifier id = cabling->cnvToIdentifier(channelID);
        const CaloDetDescrElement* caloDetElement = caloMgr->get_element(id);
        info = m_dumperTool->makeCellInfo(channelID, id, caloDetElement);
        if (!info) continue;
        m_samples->makeNewHistory(hash, info);
      }      
    }
  }
 
  std::map<HWIdentifier, LArOFIterResultsContainer::const_iterator> ofcResultPosition;
  if (m_doOFCIter && ofIterResult) { 
    for (LArOFIterResultsContainer::const_iterator ofResult = ofIterResult->begin();
         ofResult != ofIterResult->end(); ++ofResult) 
      ofcResultPosition[ofResult->getChannelID()] = ofResult;
  
    ATH_MSG_DEBUG ( "njpbSizes : " << larDigitContainer->size()
                 << " " << (ofIterResult ? ofIterResult->size() : 0) << " " 
                 << rawChannelContainer->size() << " " << channelsToKeep.size() );
  }
  SG::ReadCondHandle<CaloNoise> noiseHdl{m_noiseCDOKey};
  const CaloNoise* noiseCDO=*noiseHdl;
  
  for (LArDigitContainer::const_iterator digit = larDigitContainer->begin();
       digit != larDigitContainer->end(); ++digit) 
  {    
    //Check Energy selection
    IdentifierHash hash = m_onlineHelper->channel_Hash((*digit)->channelID());
    
    std::map<unsigned int, const LArRawChannel*>::const_iterator findChannel = channelsToKeep.find(hash);
    if (findChannel == channelsToKeep.end()) continue;
    const LArRawChannel* rawChannel = findChannel->second;
 
    //Check detector part
    HWIdentifier channelID = (*digit)->hardwareID();
    if ((m_onlineHelper->isEMBchannel(channelID) || m_onlineHelper->isEMECchannel(channelID)) && !m_doEM) continue;
    if (m_onlineHelper->isHECchannel(channelID) && !m_doHEC) continue;
    if (m_onlineHelper->isFCALchannel(channelID) && !m_doFCAL) continue;
 
    //Check gain
    CaloGain::CaloGain gain=(*digit)->gain();
    if (gain >= CaloGain::LARNGAIN || m_gains[gain] == false) continue;
 
    //Check if connected
    const bool connected = cabling->isOnlineConnected(channelID);
    if (!connected && !m_dumpDisc) continue;
   
    // Check ADCMax selection
    float pedestal = pedestals->pedestal(channelID, gain);
    float pedestalRMS = pedestals->pedestalRMS(channelID, gain);
    if (m_minADCMax > 0 || m_noiseSignifCut > 0) {
      const std::vector<short>& samples = (*digit)->samples();
      double maxValue = -1;
      for (short sample : samples)
        if (sample - pedestal > maxValue) maxValue = sample - pedestal;
      if (m_minADCMax > 0 && fabs(maxValue) < m_minADCMax) continue;
      if (m_noiseSignifCut > 0 && fabs(maxValue) < pedestalRMS*m_noiseSignifCut) continue;
    }
   
    const Identifier id = cabling->cnvToIdentifier(channelID);
    const CaloDetDescrElement* caloDetElement = nullptr;
  
    HistoryContainer* histCont = m_samples->hist_cont(hash);
    CellInfo* info = nullptr;
    if (!histCont) {
      if (!caloDetElement) caloDetElement = caloMgr->get_element(id);
      info = m_dumperTool->makeCellInfo(channelID, id, caloDetElement);
      if (!info) continue;
      histCont = m_samples->makeNewHistory(hash, info);
    }
    else 
      info = histCont->cell_info();
 
    float noise = -1;
    unsigned int status = 0xFFFFFFFF;
    if (connected) {
      if (!caloDetElement) caloDetElement = caloMgr->get_element(id);
      noise = noiseCDO->getNoise(id,gain);
      status = bcCont->status(channelID).packedData();
      HWIdentifier febId = m_onlineHelper->feb_Id(m_onlineHelper->feedthrough_Id(channelID), m_onlineHelper->slot(channelID));
      std::map<unsigned int,uint16_t>::const_iterator findError = febErrorMap.find(febId.get_identifier32().get_compact());
      if (findError != febErrorMap.end()) {
    unsigned int febErrWord = findError->second;
    status = status & (febErrWord >> 17);
      }
    }
 
    //std::vector<float> autoCorr;
    ILArAutoCorr::AutoCorrRef_t autoCorr;
    if (m_dumperTool->doShape()) {
      const LArAutoCorrComplete* autoCorrObj = dynamic_cast<const LArAutoCorrComplete*>(aCorr);
      if (!autoCorrObj)
        ATH_MSG_WARNING ( "AutoCorr object is not of type LArAutoCorrComplete!" );
      else
        autoCorr = autoCorrObj->autoCorr(channelID, gain);
    }
 
    if (!info->shape((*digit)->gain())) // this happens if doAllShapes is off
      info->setShape((*digit)->gain(), m_dumperTool->retrieveShape(channelID, gain));
    
    if (!eventData) {
      eventIndex = makeEvent(eventData, run, event, lumiBlock, bunchId); // this happens if doAllEvents is off
      if (eventIndex < 0) return StatusCode::FAILURE;
    }
    
    DataContainer* data = 
        new DataContainer((*digit)->gain(), (*digit)->samples(),
                             rawChannel->energy(),
                             rawChannel->time()/double(1000),
                             rawChannel->quality(),    
                             eventIndex,
                             autoCorr, 
                 noise, pedestal, pedestalRMS, status);
    
   //  std::map<HWIdentifier, LArOFIterResultsContainer::const_iterator>::const_iterator findResult = ofcResultPosition.find(channelID);    
//     if (findResult != ofcResultPosition.end()) {
//       LArOFIterResultsContainer::const_iterator ofResult = findResult->second;
//       if (ofResult->getValid() && ofResult->getConverged())
//         data->setADCMax(ofResult->getAmplitude());
//     }
//     else
//       msg() << MSG::INFO << "OFResult for channel 0x" << MSG::hex << channelID << MSG::dec 
//           << " not found. (size was " << ofcResultPosition.size() << ")" << endmsg;
    
   
    const auto ramp=adc2MeV->ADC2MEV(channelID,gain); //dudu
    data->setADCMax(rawChannel->energy()/ramp[1]); //pow(ADCPeak,i); //dudu
        
 
    histCont->add(data);
  }
  
  if(m_doSC) {
    std::map<unsigned int, const LArRawSC*> scToKeep;
    
    SG::ReadHandle<LArRawSCContainer> hdlSC(m_rawscKey, ctx);
    if(!hdlSC.isValid()) {
       ATH_MSG_WARNING( "Unable to retrieve LArRawSCContainer with key " << m_rawscKey << " from EventStore. " );
       return StatusCode::SUCCESS;
    } else
        ATH_MSG_DEBUG( "Got LArRawSCContainer with key " << m_rawscKey.key() );
    const LArRawSCContainer*  etcontainer = &(*hdlSC);
 
    if (etcontainer->empty()) {
      ATH_MSG_WARNING ( "LArRawSCContainer with key=" << m_rawscKey.key() << " is empty!" );
      return StatusCode::SUCCESS;
    }
 
    SG::ReadCondHandle<LArOnOffIdMapping> cablingHdl{m_cablingKeySC, ctx};
    const LArOnOffIdMapping* cablingSC=*cablingHdl;
    if(!cablingSC) {
       ATH_MSG_ERROR( "Do not have cabling object LArOnOffIdMapping" );
       return StatusCode::FAILURE;
    }
 
    SG::ReadCondHandle<CaloSuperCellDetDescrManager> caloMgrHandle{m_caloSuperCellMgrKey};
    ATH_CHECK(caloMgrHandle.isValid());
    const CaloSuperCellDetDescrManager* caloMgrSC = *caloMgrHandle;  
 
    SG::ReadCondHandle<ILArPedestal> pedHdl(m_pedestalKeySC, ctx);
    const ILArPedestal* pedestalsSC=*pedHdl;
    if (!pedestals) {
      ATH_MSG_ERROR("Failed to retrieve pedestal cond obj for SC");
       return StatusCode::FAILURE;
    }
 
    std::map<unsigned int, std::pair<float,float> > channelsToKeepSC;
 
    for (const LArRawSC* rawSC : *hdlSC) {
 
      const std::vector<unsigned short>& bcids = rawSC->bcids();
      const std::vector<int>& energies = rawSC->energies();
      const std::vector<int>& tauenergies = rawSC->tauEnergies();
      const std::vector<bool>& satur = rawSC->satur();
  
      // Look for bcid:
      float scEne = 0;
      float scTim = -99999999.;
  
      const size_t nBCIDs = bcids.size();
      size_t i = 0;
      for (i = 0; i < nBCIDs && bcids[i] != bunchId; i++)
        ;
      if(i==nBCIDs) continue;
      if (satur[i]) continue;
  
      scEne = energies[i]; 
      if (m_energyCutSC > 0 && TMath::Abs(scEne) < m_energyCut) continue;
      if (m_bcMaskSC.cellShouldBeMasked(bcCont,rawSC->hardwareID())) continue;
 
      IdentifierHash hash = m_onlineHelperSC->channel_Hash(rawSC->hardwareID());
      
      if (!hash.is_valid()) {
        ATH_MSG_FATAL ( "Found a LArRawSC whose HWIdentifier (" << rawSC->hardwareID()
                        << ") does not correspond to a valid hash -- returning StatusCode::FAILURE." );
        return StatusCode::FAILURE;
      }    
 
      if(tauenergies.size() && scEne != 0) scTim = tauenergies[i] / scEne;
      channelsToKeepSC[hash] = std::make_pair(scEne, scTim);
 
      if (m_dumpChannelInfos) {
        HistoryContainer* histCont = m_samples->hist_cont_sc(hash);
        CellInfo* info = nullptr;
        if (!histCont) {
          HWIdentifier channelID = rawSC->hardwareID();
          const Identifier id = cablingSC->cnvToIdentifier(channelID);
          const CaloDetDescrElement* caloDetElement = caloMgrSC->get_element(id);
          info = m_dumperToolSC->makeCellInfo(channelID, id, caloDetElement);
          if (!info) continue;
          m_samples->makeNewHistorySC(hash, info);
        }      
      }
    }
 
    
    SG::ReadHandle<LArDigitContainer> hdlSCDigit(m_digitsKeySC, ctx);
    if(!hdlSCDigit.isValid()) {
       ATH_MSG_WARNING( "Unable to retrieve LArDigitContainer with key " << m_digitsKeySC << " from DetectorStore. " );
       return StatusCode::SUCCESS;
    } else
        ATH_MSG_DEBUG( "Got LArDigitContainer with key " << m_digitsKeySC.key() );
    const LArDigitContainer* larSCDigitContainer = &(*hdlSCDigit);
 
    if (larSCDigitContainer->empty()) {
      ATH_MSG_WARNING ( "LArDigitContainer with key=" << m_digitsKeySC.key() << " is empty!" );
      m_nNoDigitsSC++;
      return StatusCode::SUCCESS;
    }
 
    for (LArDigitContainer::const_iterator digit = larSCDigitContainer->begin();
         digit != larSCDigitContainer->end(); ++digit) 
    {    
      //Check Energy selection
      IdentifierHash hash = m_onlineHelperSC->channel_Hash((*digit)->channelID());
      
      std::map<unsigned int, std::pair<float,float> >::const_iterator findChannel = channelsToKeepSC.find(hash);
      if (findChannel == channelsToKeepSC.end()) continue;
 
      // Check ADCMax selection
      float pedestal = pedestalsSC->pedestal((*digit)->channelID(), 0);
      if (m_minADCMaxSC > 0 ) {
        const std::vector<short>& samples = (*digit)->samples();
        double maxValue = -1;
        for (short sample : samples)
          if (sample - pedestal > maxValue) maxValue = sample - pedestal;
        if (fabs(maxValue) < m_minADCMax) continue;
      }
     
      const Identifier id = cablingSC->cnvToIdentifier((*digit)->channelID());
      const CaloDetDescrElement* caloDetElement = nullptr;
    
      HistoryContainer* histCont = m_samples->hist_cont_sc(hash);
      CellInfo* info = nullptr;
      if (!histCont) {
        if (!caloDetElement) caloDetElement = caloMgrSC->get_element(id);
        info = m_dumperToolSC->makeCellInfo((*digit)->channelID(), id, caloDetElement);
        if (!info) continue;
        histCont = m_samples->makeNewHistorySC(hash, info);
      }
      else 
        info = histCont->cell_info();
 
      if (!eventData) {
        eventIndex = makeEvent(eventData, run, event, lumiBlock, bunchId); // this happens if doAllEvents is off
        if (eventIndex < 0) return StatusCode::FAILURE;
      }
      
      DataContainer* data = 
          new DataContainer((*digit)->gain(), (*digit)->samples(),
                               findChannel->second.first,
                               findChannel->second.second/double(1000),
                               -1,    
                               eventIndex,
                               LArVectorProxy(), 
                     -1, pedestal);
      
     
      histCont->add(data);
    }
    
  } //m_doSC
  
  return StatusCode::SUCCESS;
}
 
 
StatusCode LArShapeDumper::stop()
{
  m_samples->addRun(m_runData.release());
  return StatusCode::SUCCESS;
}
 
 
StatusCode LArShapeDumper::finalize()
{
  ATH_MSG_DEBUG ("in finalize() ");
 
  if (m_prescale>1) ATH_MSG_INFO("Prescale dropped " <<  m_nPrescaledAway << "/" << m_count << " events"); 
  if (m_onlyEmptyBC) ATH_MSG_INFO("Dropped " << m_nWrongBunchGroup << "/" << m_count << " events b/c of wrong bunch group"); 
  ATH_MSG_INFO("Dropped " << m_nLArError << "/" <<  m_count << " Events b/c of LAr Veto (Noise burst or corruption)");
 
 
  int n = 0;
  for (unsigned int i = 0; i < m_samples->nChannels(); i++)
    if (m_samples->historyContainer(i)) {
      if (m_samples->historyContainer(i)->cellInfo() == nullptr)
    ATH_MSG_INFO ( "Cell with no cellInfo at index " << i << " !!" );
      //else if (m_samples->historyContainer(i)->cellInfo()->shape() == 0)
    //msg() << MSG::INFO << "Cell with no ShapeInfo at index " << i << " !!" << endmsg;
      //msg() << MSG::INFO << "Non-zero cell at index " << i << " " << m_samples->shape(i)->size() << endmsg;
      n++;
    }
 
  //for (unsigned int i = 0; i < m_samples->nEvents(); i++) {
  //   msg() << MSG::INFO << "Event " << i << " = " 
  //            << m_samples->eventData(i)->run() << " " << m_samples->eventData(i)->event()
  //            << "trigger = " << m_samples->eventData(i)->triggers() << ", nRoIs = " << m_samples->eventData(i)->nRoIs() << endmsg;
  // }
  ATH_MSG_INFO ( "Non-zero cells = " << n << ", footprint = " << m_samples->footprint() );
  ATH_MSG_INFO ( "Writing..." );
 
  if (!m_doStream) {
    m_samples->writeTrees(m_fileName.c_str());
/*    TFile* f = TFile::Open(m_fileName.c_str(), "RECREATE");
    msg() << MSG::INFO << "Writing (2)..." << endmsg;
    f->WriteObjectAny(m_samples, "Container", "LArSamples");
    msg() << MSG::INFO << "Closing..." << endmsg;
    f->Close();
    msg() << MSG::INFO << "Deleting..." << endmsg;
    delete m_samples;*/
    msg() << MSG::INFO << "Done!" << endmsg;
  }
 
  return StatusCode::SUCCESS;
}
 
 
int LArShapeDumper::makeEvent(EventData*& eventData,
                          int run, int event, 
                  int lumiBlock, int bunchXing) const
{
  std::vector<unsigned int> triggerWords;
  if (m_doTrigger) {
    const ROIB::RoIBResult* l1Result = nullptr;
    if (evtStore()->retrieve(l1Result).isFailure() || !l1Result) {
      ATH_MSG_FATAL ( "Could not retrieve RoIBResult!" );
      return -1;
    }
    const std::vector<ROIB::CTPRoI> tav = l1Result->cTPResult().TAV();
    for (const ROIB::CTPRoI& word : tav)
      triggerWords.push_back(word.roIWord());
 
    for (const std::pair<const TString, unsigned int>& p : m_runData->triggerConfig()) {
      while (triggerWords.size() <= p.second/32) triggerWords.push_back(0);
      if (m_trigDec->isPassed(p.first.Data())) {
    triggerWords[p.second/32] |= (0x1 << (p.second % 32));
      //msg() << MSG::INFO << "Trigger line " << p.first.Data() << " passed" << endmsg;
      }
    }
    //msg() << MSG::INFO << "Trigger words : ";
    //for (unsigned int i = 0; i < triggerWords.size(); i++) msg() << MSG::INFO << triggerWords[i] << " ";
    //msg() << MSG::INFO << endmsg;
  }
  
  eventData = new EventData(event, 0, lumiBlock, bunchXing);
  if (m_runData->run() == 0) m_runData->setRun(run);
  eventData->setRunData(m_runData.get());
  eventData->setTriggerData(triggerWords);
  if (m_doRoIs) {
    //msg() << MSG::INFO << "Filling RoI list" << endmsg;
    for (const Trig::ChainGroup* group : m_triggerGroups) {
      std::vector<Trig::Feature<TrigRoiDescriptor> > roIs = group->features().get<TrigRoiDescriptor>();
      for (const Trig::Feature<TrigRoiDescriptor>& roI : roIs) {
    //msg() << MSG::INFO << "Found an roi for chain ";
        //for (unsigned int i = 0; i < group->getListOfTriggers().size(); i++) cout << group->getListOfTriggers()[i] << " ";
        //cout << "@ " << roI.cptr()->eta() << ", " << roI.cptr()->phi() << ", TE = " 
    //   << roI.te()->getId() << " " << Trig::getTEName(*roI.te()) << " with label " << roI.label() << endmsg;
    eventData->addRoI(roI.cptr()->eta(), roI.cptr()->phi(), group->getListOfTriggers()[0].c_str(), roI.label().c_str());
    //msg() << MSG::INFO << "nRoIs so far = " << eventData->nRoIs() << endmsg;
      }
    }
  }
  return m_samples->addEvent(eventData);
}