LArRodBlockPhysicsV5.cxx

Go to the documentation of this file.

//Dear emacs, this is -*- c++ -*-
 
/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
// Implementation of a LArRODBlockStructure class
// This version contains LArDigits in fixed gain.
// See .h file for more details.
 
#include "GaudiKernel/MsgStream.h"
#include "AthenaKernel/getMessageSvc.h"
#include "LArRawEvent/LArDigit.h"
#include "LArByteStream/LArRodBlockPhysicsV5.h"
#include "StoreGate/StoreGateSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
#include <iostream>
 
 
 
 
LArRodBlockPhysicsV5::LArRodBlockPhysicsV5(IMessageSvc* msgSvc)
  : LArRodBlockStructure(msgSvc, BlockType()),
    m_onlineHelper(nullptr)
{
  m_iHeadBlockSize=endtag/2; // The implicit cast rounds down to the right size 
  m_fixedGain=CaloGain::LARNGAIN;
  LArRodBlockPhysicsV5::resetPointers();
  m_EnergyThreshold1=100;
  m_EnergyThreshold2=0;
  m_requiredNSamples = 0;
  m_OffTimeCut=0;  //FIXME: Nowhere set to a sensible value ???
  m_numberHotCellOffTime=0;
  // retrieve onlineHelper
  SmartIF<StoreGateSvc> detStore{Gaudi::svcLocator()->service("DetectorStore")};
  if (!detStore) {
    m_logstr << MSG::ERROR << "Unable to locate DetectorStore" << endmsg;
    std::abort();
  }
  StatusCode sc = detStore->retrieve(m_onlineHelper, "LArOnlineID");
  if (sc.isFailure()) {
    m_logstr << MSG::ERROR << "Could not get LArOnlineID helper !" << endmsg;
    std::abort();
  }
}
 
 
void LArRodBlockPhysicsV5::resetPointers() 
{
  m_EnergyIndex=0;
  m_TimeQualityIndex=0;
  m_DigitsIndex=0;
  m_DigitsChannel=0;
  m_RawDataIndex=0;
 
  m_GainPointer=0;
  m_MaskTimeQualityPointer=0;
  m_MaskDigitsPointer=0;
  m_RaddPointer=0;
  m_EnergyPointer=0;
  m_SumPointer=0;
  m_TimeQualityPointer=0;
  m_DigitsPointer=0;
  m_RawDataPointer=0;
  m_numberHotCell=0;
}
 
bool LArRodBlockPhysicsV5::setPointers()
{
  if (m_FebBlockSize>m_iHeadBlockSize)
    {
      int off = -8;
      int ns   = getHeader16(NSamples) & 0xff;
      if (m_requiredNSamples > 0 && m_requiredNSamples != ns) return false;
      int radd = (ns+1)/2;
      int dim1 = getHeader16(ResultsDim1);
      int off1 = getHeader16(ResultsOff1);
      int off2 = getHeader16(ResultsOff2);
      int dim2 = getHeader16(ResultsDim2);
      int off3 = getHeader16(RawDataBlkOff);
      int dim3 = getHeader16(RawDataBlkDim);
 
      if (off1 && dim1+off1+off<m_FebBlockSize) {
    off1 += off;
    if (dim1>=8) 
      m_GainPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+off1);
    if (dim1>=12) 
      m_MaskTimeQualityPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+off1+8);
    if (dim1>=16) 
      m_MaskDigitsPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+off1+12);
    if (dim1>=16+radd) 
      m_RaddPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off1+16);
    if (dim1>=80+radd)
      m_EnergyPointer=reinterpret_cast<const uint16_t*> (m_FebBlock+off1+16+radd);
    if (dim1>=83+radd) 
      m_SumPointer=reinterpret_cast<const int32_t*>(m_FebBlock+off1+80+radd);
    if (dim1>83+radd) 
      m_TimeQualityPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off1+83+radd);
    off1 -= off;
      }
      if (off2 && dim2+off2+off<m_FebBlockSize) {
    m_DigitsPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off2+off);
      }
      if (off3 && dim3+off3+off<m_FebBlockSize) {
    m_RawDataPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off3+off);
      }
 
      // Check for offsets problems
      uint32_t problem = 0;
      int n1, n2;
      int n1_tmp, n2_tmp;
      int off1_tmp, dim1_tmp;
      int off2_tmp, dim2_tmp;
      int off3_tmp, dim3_tmp;
      if(off1==0) {
    n1 = n2 = 0;
    n1_tmp = n2_tmp =0;
    off1_tmp = dim1_tmp = 0; 
    off2_tmp = dim2_tmp = 0; 
    off3_tmp = dim3_tmp = 0; 
      }
      else { 
    m_RaddPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+26);
    m_MaskTimeQualityPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+18);
    m_MaskDigitsPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+22);
    n1 = getNbSweetCells1();
    n2 = getNbSweetCells2();
    n1_tmp = getNbSweetCells1FromMask();
    n2_tmp = getNbSweetCells2FromMask();
    off1_tmp = 10-off; 
    dim1_tmp = 83+(ns+1)/2+n1; 
    if ( m_requiredNSamples > 0 )
        dim1_tmp = 83 +(m_requiredNSamples+1)/2+n1;
    off2_tmp = off1_tmp+dim1_tmp;
    dim2_tmp = (n2*ns+1)/2;
    off3_tmp = off2_tmp+dim2_tmp;
    dim3_tmp = getNumberOfWords()-3-off3_tmp-off;
    if(dim2_tmp==0) off2_tmp = 0;
    if(dim3_tmp==0) off3_tmp = 0;
      }
 
      if(off1 != off1_tmp) problem=1;
      if(dim1 != dim1_tmp) problem=2;
      if(off2 != off2_tmp) problem=3;
      if(dim2 != dim2_tmp) problem=4;
      if(off3 != off3_tmp) problem=5;
      if(dim3 != dim3_tmp) problem=6;
      if(n1 != n1_tmp)     problem=7;
      if(n2 != n2_tmp)     problem=8;
      if (m_requiredNSamples > 0 && 
          getHeader32(NGains) != (uint32_t)0x10000 + m_requiredNSamples) problem=9;
      if(problem) { // Try to recompute offsets
        m_logstr << MSG::ERROR << "LArByteStreamProblem " << problem << endmsg;
        m_logstr << MSG::ERROR << "NSamples = " << std::dec << ns << endmsg;
        m_logstr << MSG::ERROR << "getHeader32(NGains) = " << std::hex << getHeader32(NGains) << endmsg;
    m_logstr << MSG::ERROR << "NWTot:   " << std::hex  << getNumberOfWords() << " n1=" << n1 << " (" << n1_tmp << ") n2=" << n2 << " (" << n2_tmp << ")" << endmsg;
    m_logstr << MSG::ERROR << "Found 1: " << off1 << " " << dim1 << endmsg;
    m_logstr << MSG::ERROR << "Found 2: " << off2 << " " << dim2 << endmsg;
    m_logstr << MSG::ERROR << "Found 3: " << off3 << " " << dim3 << std::dec << endmsg;
    
    if(n1==n1_tmp && n2==n2_tmp) {  // Check consistency of cells above threshold
      off1 = off1_tmp;
      dim1 = dim1_tmp;
      off2 = off2_tmp;
      dim2 = dim2_tmp;
      off3 = off3_tmp;
      dim3 = dim3_tmp;
      m_logstr << MSG::ERROR << "Recomputed 1: " << std::hex << off1 << " " << dim1 << endmsg;
      m_logstr << MSG::ERROR << "Recomputed 2: " << off2 << " " << dim2 << endmsg;
      m_logstr << MSG::ERROR << "Recomputed 3: " << off3 << " " << dim3 << std::dec << endmsg;
 
      if (off1 && dim1+off1+off<m_FebBlockSize) {
        off1 += off;
        if (dim1>=8) 
          m_GainPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+off1);
        if (dim1>=12) 
          m_MaskTimeQualityPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+off1+8);
        if (dim1>=16) 
          m_MaskDigitsPointer=reinterpret_cast<const uint32_t*>(m_FebBlock+off1+12);
        if (dim1>=16+radd) 
          m_RaddPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off1+16);
        if (dim1>=80+radd)
          m_EnergyPointer=reinterpret_cast<const uint16_t*> (m_FebBlock+off1+16+radd);
        if (dim1>=83+radd) 
          m_SumPointer=reinterpret_cast<const int32_t*>(m_FebBlock+off1+80+radd);
        if (dim1>83+radd) 
          m_TimeQualityPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off1+83+radd);
      }
      if (off2 && dim2+off2+off<m_FebBlockSize) {
        m_DigitsPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off2+off);
      }
      if (off3 && dim3+off3+off<m_FebBlockSize) {
        m_RawDataPointer=reinterpret_cast<const uint16_t*>(m_FebBlock+off3+off);
      }
    }
      }
 
      problem=0;
      // Recheck offsets
      if(off1< off2 && off1 + dim1 > off2) problem = 1;
      if(off1< off3 && off1 + dim1 > off3) problem = 2;
      if(off2< off1 && off2 + dim2 > off1) problem = 3;
      if(off2< off3 && off2 + dim2 > off3) problem = 4;
      if(off3< off1 && off3 + dim3 > off1) problem = 5;
      if(off3< off2 && off3 + dim3 > off2) problem = 6;
 
      if(problem) {
    resetPointers();
    m_logstr << MSG::ERROR << "LArByteStreamProblem " << problem << endmsg;
    m_logstr << MSG::ERROR << "Unrecoverable problem" << endmsg;
      }
   }
 
  return true;
}
 
int LArRodBlockPhysicsV5::getNextRawData(int& channelNumber, std::vector<short>& samples, uint32_t& gain)
{
#ifdef LARBSDBGOUTPUT
  //Debug output
  m_logstr << MSG::DEBUG << "Let s go in getNextRawData..." << endmsg;
  m_logstr << MSG::DEBUG << "GetNextRawData for FEB 0x" << MSG::hex << (uint32_t)getHeader32(FEBID) << MSG::dec << endmsg;
  m_logstr << MSG::DEBUG << "m_RawDataPointer=" << m_RawDataPointer << " m_RawDataIndex="<<  m_RawDataIndex
     << " m_channelsPerFEB=" << m_channelsPerFEB << endmsg;
#endif
 
  if (m_RawDataIndex>=m_channelsPerFEB) { //Already beyond maximal number of channels
#ifdef LARBSDBGOUTPUT
    m_logstr << MSG::DEBUG << "Maximum number of channels reached" << endmsg;
#endif
    return 0;
  }
  if (!m_RawDataPointer) { //Block does not exist
    // Try to get samples and gain from getNextDigits
    return getNextDigits(channelNumber,samples,gain);
  }
 
  // Get next channel
  unsigned rodChannelNumber=m_RawDataIndex;      // Index of Channel in ROD-Block
  channelNumber=((rodChannelNumber&0xe)<<2) + ((rodChannelNumber&0x1)<<6) + (rodChannelNumber>>4);    //channel number of the FEB
  uint32_t febgain;
  const unsigned int nsamples = getHeader16(NSamples) & 0xff;
  const unsigned int ngains   = getHeader16(NGains);
 
#ifdef LARBSDBGOUTPUT
  m_logstr << MSG::DEBUG << "This FEB has " << nsamples <<  " samples" << endmsg;
  m_logstr << MSG::DEBUG << "This FEB has " << ngains   <<  " gains" << endmsg;
#endif
 
  if(ngains==0 || nsamples==0) return 0;
  int s_size = nsamples+1;
  int offset = (10+nsamples)&0xfffc;
  int index;
  index = s_size*m_RawDataIndex + offset;
  uint16_t s[2];
  if((nsamples+1)&0x7) {
     s[0]    = m_RawDataPointer[index++]>>2;
     febgain = m_RawDataPointer[index++];
     samples.push_back(s[0]);
     for(unsigned int i=0;i<nsamples/2;i++) {
      s[1]    = m_RawDataPointer[index++]>>2;
      s[0]    = m_RawDataPointer[index++]>>2;
      samples.push_back(s[0]);
      samples.push_back(s[1]);
     }
   } // End of check for 5 samples
   else {
  if (!(m_RawDataIndex%2)) {
  s[0]    = m_RawDataPointer[index++]>>2;
  febgain = m_RawDataPointer[index++];
  samples.push_back(s[0]);
  for(unsigned int i=0;i<nsamples/2;i++) {
    s[1]    = m_RawDataPointer[index++]>>2;
    s[0]    = m_RawDataPointer[index++]>>2;
    samples.push_back(s[0]);
    samples.push_back(s[1]);
  }
  } else {
  for(unsigned int i=0;i<nsamples/2;i++) {
    s[1]    = m_RawDataPointer[index++]>>2;
    s[0]    = m_RawDataPointer[index++]>>2;
    samples.push_back(s[0]);
    samples.push_back(s[1]);
  }
  febgain = m_RawDataPointer[index++];
  s[0]    = m_RawDataPointer[index++]>>2;
  samples.push_back(s[0]);
  }
  } // End of >5 check
  gain=RawToOfflineGain(febgain);
 
#ifdef LARBSDBGOUTPUT
  m_logstr << MSG::DEBUG << " ===> ROD Channel = " << m_RawDataIndex << endmsg;
  m_logstr << MSG::DEBUG << " ===> FEB Channel = " << channelNumber << endmsg; 
  m_logstr << MSG::DEBUG << " ===> Gain        = " << gain << endmsg;
  for(int i=0;i<nsamples;i++)
    m_logstr << MSG::DEBUG << " ===> sample " << i << "    = " << samples[i] << endmsg;
  int n = m_RawDataIndex;
  int32_t e,t,q;
  uint32_t g;
  LArRodBlockPhysicsV5::getNextEnergy(n,e,t,q,g);
#endif
  ++m_RawDataIndex;
  unsigned  rearrangeFirstSample=0;
  if (m_rearrangeFirstSample)
    rearrangeFirstSample=m_rearrangeFirstSample; //Overwrite by jobOptions
  else
    rearrangeFirstSample=getFirstSampleIndex();
  if (rearrangeFirstSample && rearrangeFirstSample<samples.size()) //FIXME: Very ugly hack! See explanation in LArRodDecoder.h file
  {//Change e.g. 3 0 1 2 4 to 0 1 2 3 4
    short movedSample=samples[0];
    for (unsigned i=1;i<=rearrangeFirstSample;i++)
      samples[i-1]=samples[i];
    samples[rearrangeFirstSample]=movedSample;
  }
#ifdef LARBSDBGOUTPUT
  m_logstr << MSG::DEBUG << "GetNextRawData for FEB finished 0x" << MSG::hex << (uint32_t)getHeader32(FEBID) << MSG::dec << endmsg;
#endif
  return 1;
}
 
int LArRodBlockPhysicsV5::getNextDigits(int& channelNumber, std::vector<short>& samples, uint32_t& gain)
{
#ifdef LARBSDBGOUTPUT
  //Debug output
  m_logstr << MSG::DEBUG << "Let s go in getNextDigits..." << endmsg;
  m_logstr << MSG::DEBUG << "GetNextDigits for FEB 0x" << MSG::hex << (uint32_t)getHeader32(FEBID) << MSG::dec << endmsg;
  m_logstr << MSG::DEBUG << "m_DigitsPointer=" << m_DigitsPointer << " m_DigitsIndex="<<  m_DigitsIndex
     << " m_DigitsChannel="<<  m_DigitsChannel 
     << " m_channelsPerFEB=" << m_channelsPerFEB << endmsg;
#endif
 
  if (m_DigitsChannel>=m_channelsPerFEB) { //Already beyond maximal number of channels
#ifdef LARBSDBGOUTPUT
    m_logstr << MSG::DEBUG << "Maximum number of channels reached" << endmsg;
#endif
    return 0;
  }
  if (!m_DigitsPointer) { //Block does not exist
#ifdef LARBSDBGOUTPUT
    m_logstr << MSG::DEBUG << "No Digits Block in this FEB" << endmsg;
#endif
    return 0; 
  }
  if (!m_MaskDigitsPointer) { //Block does not exist
#ifdef LARBSDBGOUTPUT
    m_logstr << MSG::DEBUG << "No Mask Digits Block in this FEB" << endmsg;
#endif
    return 0; 
  }
 
  // Get Digits if the information is present according to summary block
  uint32_t hasDigits;
 
  hasDigits = (uint32_t) ((m_MaskDigitsPointer[m_DigitsChannel>>5] >> (m_DigitsChannel&0x1f)) &0x1);
  // Increment channel number until digits are found
  while(hasDigits==0) {
    m_DigitsChannel++;
    if (m_DigitsChannel>=m_channelsPerFEB) { //Already beyond maximal number of channels
#ifdef LARBSDBGOUTPUT
      m_logstr << MSG::DEBUG << "Maximum number of channels reached" << endmsg;
#endif
      return 0;
    }
    hasDigits = (uint32_t) ((m_MaskDigitsPointer[m_DigitsChannel>>5] >> (m_DigitsChannel&0x1f)) &0x1);
  }
 
  // Get next channel
  unsigned rodChannelNumber=m_DigitsChannel;      // Index of Channel in ROD-Block
  channelNumber=((rodChannelNumber&0xe)<<2) + ((rodChannelNumber&0x1)<<6) + (rodChannelNumber>>4);    //channel number of the FEB
  const unsigned int nsamples = getHeader16(NSamples) & 0xff;
 
  // gain in 2 bits of a 32 bits word 
  if(m_GainPointer) {
    gain = (uint32_t) ((m_GainPointer[m_DigitsChannel>>4] >> (m_DigitsChannel&0xf)*2) & 0x3);
    gain=RawToOfflineGain(gain);
  } else gain=0xffffffff;
 
#ifdef LARBSDBGOUTPUT
  m_logstr << MSG::DEBUG << "This FEB has " << nsamples <<  " samples" << endmsg;
#endif
 
  if(nsamples==0) return 0;
  int s_size = nsamples;
  int index;
  index = s_size*m_DigitsIndex;
  if(m_DigitsIndex&0x1) {
    samples.push_back(m_DigitsPointer[index-1]>>2);
    samples.push_back(m_DigitsPointer[index+2]>>2);
    samples.push_back(m_DigitsPointer[index+1]>>2);
    samples.push_back(m_DigitsPointer[index+4]>>2);
    samples.push_back(m_DigitsPointer[index+3]>>2);
    if(nsamples==7) {
      samples.push_back(m_DigitsPointer[index+6]>>2);
      samples.push_back(m_DigitsPointer[index+5]>>2);
    }
  } else {
    samples.push_back(m_DigitsPointer[index+1]>>2);
    samples.push_back(m_DigitsPointer[index+0]>>2);
    samples.push_back(m_DigitsPointer[index+3]>>2);
    samples.push_back(m_DigitsPointer[index+2]>>2);
    samples.push_back(m_DigitsPointer[index+5]>>2);
    if(nsamples==7) {
      samples.push_back(m_DigitsPointer[index+4]>>2);
      samples.push_back(m_DigitsPointer[index+7]>>2);
    }
  }
 
#ifdef LARBSDBGOUTPUT
  m_logstr << MSG::DEBUG << " ===> ROD Channel = " << m_DigitsChannel << endmsg;
  m_logstr << MSG::DEBUG << " ===> FEB Channel = " << channelNumber << endmsg; 
  m_logstr << MSG::DEBUG << " ===> Gain        = " << gain << endmsg;
  for(int i=0;i<nsamples;i++)
    m_logstr << MSG::DEBUG << " ===> sample " << i << "    = " << samples[i] << endmsg;
#endif
  m_DigitsIndex++;
  m_DigitsChannel++;
  unsigned  rearrangeFirstSample=0;
  if (m_rearrangeFirstSample)
    rearrangeFirstSample=m_rearrangeFirstSample; //Overwrite by jobOptions
  else
    rearrangeFirstSample=getFirstSampleIndex();
  if (rearrangeFirstSample && rearrangeFirstSample<samples.size()) //FIXME: Very ugly hack! See explanation in LArRodDecoder.h file
  {//Change e.g. 3 0 1 2 4 to 0 1 2 3 4
    short movedSample=samples[0];
    for (unsigned i=1;i<=rearrangeFirstSample;i++)
      samples[i-1]=samples[i];
    samples[rearrangeFirstSample]=movedSample;
  }
#ifdef LARBSDBGOUTPUT
  m_logstr << MSG::DEBUG << "GetNextDigits for FEB finished 0x" << MSG::hex << (uint32_t)getHeader32(FEBID) << MSG::dec << endmsg;
#endif
  return 1;
}
 
uint16_t LArRodBlockPhysicsV5::getNbSweetCells1()  const 
{ 
  if(!m_RaddPointer) return 0;
  return m_RaddPointer[1]>>8;
}
 
uint16_t LArRodBlockPhysicsV5::getNbSweetCells2()  const 
{ 
  if(!m_RaddPointer) return 0;
  return m_RaddPointer[1] & 0xff;
}
 
uint16_t LArRodBlockPhysicsV5::getNbSweetCells1FromMask()  const 
{ 
  if(!m_MaskTimeQualityPointer) return 0;
  int n=0;
  for(int i=0;i<4;i++)
    for(int j=0;j<32;j++)
      if((m_MaskTimeQualityPointer[i] >> j) &0x1) n++;
  return n;
}
 
uint16_t LArRodBlockPhysicsV5::getNbSweetCells2FromMask()  const 
{ 
  if(!m_MaskDigitsPointer) return 0;
  int n=0;
  for(int i=0;i<4;i++)
    for(int j=0;j<32;j++)
      if((m_MaskDigitsPointer[i] >> j) &0x1) n++;
  return n;
}
 
uint32_t LArRodBlockPhysicsV5::getNumberOfSamples() const
{ 
  return getHeader16(NSamples); 
}
 
uint32_t LArRodBlockPhysicsV5::getNumberOfGains() const
{ 
  return  getHeader16(NGains);
}
 
uint16_t LArRodBlockPhysicsV5::getResults1Size() const
{
  return getHeader16(ResultsDim1);
}
 
uint16_t LArRodBlockPhysicsV5::getResults2Size() const
{
  return getHeader16(ResultsDim2);
}
 
uint16_t LArRodBlockPhysicsV5::getRawDataSize() const
{
  return getHeader16(RawDataBlkDim);
}
 
uint32_t LArRodBlockPhysicsV5::getRadd(uint32_t adc, uint32_t sample)  const
{ 
  if(!m_RawDataPointer) {
    if(!m_RaddPointer) return 0;
    if(sample%2) sample+=2;
    return m_RaddPointer[sample];
  }
  int index;
  if(sample==0)         index=6;
  else if(sample & 0x1) index=7+sample-1;
  else                  index=7+sample+1;
  uint32_t x=m_RawDataPointer[index];
  if(adc>=8) return x>>8; 
  return x&0xff;
}
 
uint16_t LArRodBlockPhysicsV5::getCtrl1(uint32_t /*adc*/)  const
{ 
  if(!m_RawDataPointer) return 0;
  int index=5;
  uint16_t x=m_RawDataPointer[index];
  return x;
}
 
uint16_t LArRodBlockPhysicsV5::getCtrl2(uint32_t /*adc*/)  const
{ 
  if(!m_RawDataPointer) return 0;
  int index=4;
  uint16_t x=m_RawDataPointer[index];
  return x;
}
 
uint16_t LArRodBlockPhysicsV5::getCtrl3(uint32_t /*adc*/)  const
{ 
  if(!m_RawDataPointer) return 0;
  int index=7;
  uint16_t x=m_RawDataPointer[index];
  return x;
}
 
uint32_t LArRodBlockPhysicsV5::getStatus()  const
{ 
  if(getNumberOfWords()<EventStatus/2) return 0;
  uint32_t x=getHeader32(EventStatus);
  return x;
}
 
// start of encoding methods
void LArRodBlockPhysicsV5::initializeFragment(std::vector<uint32_t>& fragment ){
  m_pRODblock=&fragment; //remember pointer to fragment
  if (fragment.size()>m_iHeadBlockSize) {  //Got filled fragment
    unsigned int sizeRead=0;
    //Store existing data in the FEB-Map
    while (sizeRead<fragment.size()) {
      std::vector<uint32_t>::iterator FebIter;
      FebIter=fragment.begin()+sizeRead;     //Store pointer to current Feb-Header
      m_FebBlock=&(*FebIter); //Set m_FebBlock in order to use getHeader-functions.
      uint32_t currFEBid=getHeader32(FEBID); //Get this FEB-ID
      uint16_t currFebSize=getNumberOfWords(); //Size of this FEB-Block
      if (FebIter+currFebSize>fragment.end()) {
        fragment.clear(); //Clear existing vector
        return;
      }
      m_mFebBlocks[currFEBid].assign(FebIter,FebIter+currFebSize); //Copy data from ROD-fragment into FEB-Block
      sizeRead+=currFebSize+m_MiddleHeaderSize;  //6 is the middle header size
    } // end while
  }
  fragment.clear(); //Clear existing vector
}
 
//For writing: Initalizes a single FEB-Block
void LArRodBlockPhysicsV5::initializeFEB(const uint32_t id)
{
  m_vFragment=&(m_mFebBlocks[id]);
  if (m_vFragment->size()<m_iHeadBlockSize) //Got empty or spoiled fragment
  {
    m_vFragment->resize(m_iHeadBlockSize,0); //Initialize FEB-Header
    setHeader32(FEBID,id);                //Set Feb ID
    // At least 10 (head) + 16 (gain/sumblks) + 64 (energies)
    m_vFragment->reserve(90);
  }
 
  m_SumBlkBlockE1.resize(4);
  for(unsigned int i=0;i<4;i++) m_SumBlkBlockE1[i]=0x0;
  m_SumBlkBlockE2.resize(4);
  for(unsigned int i=0;i<4;i++) m_SumBlkBlockE2[i]=0x0;
  m_GainBlock.resize(8);
  for(unsigned int i=0;i<8;i++) m_GainBlock[i]=0x0;
  m_TimeQualityBlock.resize(6);
  for(unsigned int i=0;i<6;i++) m_TimeQualityBlock[i]=0x0;
  m_EnergyBlockEncode.resize(128);
  for(unsigned int i=0;i<128;i++) m_EnergyBlockEncode[i]=0x0;
  m_DigitsEncode.clear();
 
  resetPointers();
}
 
void LArRodBlockPhysicsV5::setNextEnergy(const int channel, const int32_t energy,
                                         const int32_t time, const int32_t quality, const uint32_t gain)
{
  int rcNb=FebToRodChannel(channel);
  //rcNb is supposed to be equal or bigger than m_EIndex.
  //In the latter case, we fill up the missing  channels with zero
  if (rcNb<m_EnergyIndex) {
    return;
  }
 
  //Fill up missing channels with zeros:
  while (m_EnergyIndex<rcNb)
    setNextEnergy((int16_t)0,(int16_t)32767,(int16_t)-32767,(uint32_t)0);
 
  // transform 32 bits data into 16 bits data
 
  uint16_t theenergy;
  uint32_t abse,EncodedE;
  int16_t thetime,thequality;
  int32_t sign;
 
  //Time is in 10 ps in ByteStream, hence the factor 10 to convert from ps
  thetime = (int16_t) time/10;
  thequality = (int16_t) quality;
 
  sign=(energy>=0?1:-1); // get sign of energy
  abse=(uint32_t)abs(energy);
 
  EncodedE=abse; // range 0
 
  if ((abse>8192)&&(abse<65536))
  {
    EncodedE=((abse>>3)|0x4000); // range 1 : drop last 3 bits and put range bits (bits 14 and 13 = 01)
  }
  else if ((abse>65535)&&(abse<524288))
  {
    EncodedE=((abse>>6)|0x8000); // range 2 : drop last 6 bits and put range bits (bits 14 and 13 = 10)
  }
  else if ((abse>524288))
  {
    EncodedE=((abse>>9)|0xc000); // range 3 : drop last 9 bits and put range bits (bits 14 and 13 = 11)
  }
 
  // treat sign now :
 
  if (sign<0) EncodedE |= 0x2000;
  theenergy = (uint16_t) EncodedE;
 
  // Add data...
 
  if (abse> m_EnergyThreshold1)
  {
    setNextEnergy(theenergy,thetime,thequality,gain);
  }
  else
  {
    setNextEnergy(theenergy,(int16_t)32767,(int16_t)-32767,gain);
  }
}
 
//Private function, expects channel number is rod-style ordering
void LArRodBlockPhysicsV5::setNextEnergy(const uint16_t energy,const int16_t time, const int16_t quality, const uint32_t gain)
{
  if (m_EnergyIndex>=m_channelsPerFEB)        //Use m_EIndex to count total number of channels
    {
      return;
    }
 
  // Energy
  int endianindex;
  if (m_EnergyIndex & 0x1) endianindex = m_EnergyIndex-1;
  else                    endianindex = m_EnergyIndex+1;
  m_EnergyBlockEncode[endianindex] = energy;
 
  // Find correct position
 
  // update summary block 
  // Gain is composed of two bits per cell
  uint16_t gain_idx=m_EnergyIndex>>4;
  uint16_t gain_bit=(m_EnergyIndex&0xf)*2;
  uint32_t gain1 = OfflineToRawGain(gain);
  m_GainBlock[gain_idx] |= (gain1 << gain_bit);
 
  // write Time and Chi2 for cells above HighEnergyCellCut threshold
 
  if (quality!=-32767) // Do write Time and Chi2 information
    {
      // count the number of hot cells
      m_numberHotCell++;
      // count the number of cells offtime
      if (abs(time)>m_OffTimeCut) m_numberHotCellOffTime++;
      uint16_t mask_idx=m_EnergyIndex>>5;
      uint16_t mask_bit=(m_EnergyIndex&0x1f);
      m_SumBlkBlockE1[mask_idx] |= (0x1 << mask_bit);
 
      m_TimeQualityBlock.push_back(*((uint16_t*)&time));
      m_TimeQualityBlock.push_back(*((uint16_t*)&quality));
    }
  m_EnergyIndex++; //Use m_EIndex to count the channels put in the Energy block
 
}
 
void LArRodBlockPhysicsV5::setRawData(const int chIdx, const std::vector<short>& samples , const uint32_t /* gain_not_used */ ){
 
  // First of all, set the bits
  int cchIdx = FebToRodChannel(chIdx);
  uint16_t mask_idx=cchIdx>>5;
  uint16_t mask_bit=(cchIdx&0x1f);
  m_SumBlkBlockE2[mask_idx] |= (0x1 << mask_bit);
  for(std::vector<short>::const_iterator i=samples.begin();i!=samples.end();++i){
    m_DigitsEncode.push_back((*i)<<2);
  }
}
 
void LArRodBlockPhysicsV5::finalizeFEB()
{
  //Complete non-complete Energy block
  while (m_EnergyIndex<m_channelsPerFEB)
    setNextEnergy((uint16_t)0,(int16_t)32767,(int32_t)-32767,(uint32_t)0);//E=0,t=32767,q=-32767,G=0
 
  uint16_t n;
  uint16_t nsamples=5;
  // checkSum value
  uint32_t sum=0;
 
  // Will Hardcode here for the moment FIXME. Minimal 1 sample
  setHeader16(NGains,1);
  setHeader16(NSamples,nsamples);
  // These will never be used form MC. Nice to put in here thought
  setHeader16(FEB_SN,0xfefe);
  setHeader16(FEB_SN_h,0xdede);
  setHeader16(InFPGAFormat,0x0);
  setHeader16(InFPGAFormat_h,0x2);
 
  // Gain block...
  n = m_GainBlock.size();
  //Check if Gain-Block exists and is not yet part of the fragment
  if (n)
  {
    for(unsigned int i=0;i<n;i++){
      m_vFragment->push_back(m_GainBlock[i]);
      sum+=m_GainBlock[i];
    }
  }
 
  // Cells above energy threshold E1
  n = m_SumBlkBlockE1.size();
  //Check if Summary Block exists and is not yet part of the fragment
  if (n)
  {
    for (unsigned i=0;i<n;i++){
      m_vFragment->push_back(m_SumBlkBlockE1[i]);
      sum+=m_SumBlkBlockE1[i];
    }
  }
 
  // Cells above energy threshold E2 (not included so far)
  n = m_SumBlkBlockE2.size();
  //Check if Summary Block exists and is not yet part of the fragment
  if (n)
  {
    for (unsigned i=0;i<n;i++){
      m_vFragment->push_back(m_SumBlkBlockE2[i]);
      sum+=m_SumBlkBlockE2[i];
    }
  }
 
  // fill info from counters
  // for moment just include 1 fake words (32 bits) to put radd
  uint32_t radd_nANC=0x0;
  // Second threshold missing (FIXME)
  radd_nANC = ((m_numberHotCell<<8))+(m_DigitsEncode.size()/nsamples);
  radd_nANC = (radd_nANC<<16);
  m_vFragment->push_back(radd_nANC);
  sum+=radd_nANC;
  // Need to include radd nsamples-1
  // No need to include in sum's for now
  for( int i=0; i < (nsamples-1)/2; i++)
    m_vFragment->push_back(0x0);
 
  // Energy block...
  n = 128 ; // Fixed size m_EnergyBlock.size();
  // Block also include time, whenever necessary
  int size_of_block=80+(nsamples+1)/2+(m_TimeQualityBlock.size())/2;
  //Check if Energy-Block exists and is not yet part of the fragment
  if (n)
  {
    setHeader16(ResultsOff1,18);
    setHeader16(ResultsDim1,size_of_block);
    for(unsigned int i=0;i<n/2;i++) {
      // WARNING which one should be >>16 2*i or 2*i+1? To be tested
      uint32_t Encode = m_EnergyBlockEncode[2*i]+(m_EnergyBlockEncode[2*i+1]<<16);
      m_vFragment->push_back(Encode);
      sum+=Encode;
    }
  }
 
  // Magic numbers (3 or 6) for Ex, Ey and Ez
  n = m_TimeQualityBlock.size();
  //Check if Time and Quality Block exists and is not yet part of the fragment
  if (n){
    unsigned int imax = n/2;
    for (unsigned int i=0;i<imax;i++) {
      const auto low = static_cast<std::uint32_t>(m_TimeQualityBlock[i * 2]);
      const auto high = static_cast<std::uint32_t>(m_TimeQualityBlock[i * 2 + 1]);
      const std::uint32_t to_push = low | (high << 16);
      m_vFragment->push_back(to_push);
      sum += to_push;
    }
  }
  // Now include digits
  n = m_DigitsEncode.size();
  if ( n ) {
    // First make sure it is not and odd number to store
    if ( m_DigitsEncode.size() & 0x1 ) m_DigitsEncode.push_back(0x0);
    unsigned int imax=m_DigitsEncode.size()/2;
    for (unsigned int i=0;i<imax;i++) {
      const auto low = static_cast<std::uint32_t>(m_DigitsEncode[i * 2 + 1]);
      const auto high = static_cast<std::uint32_t>(m_DigitsEncode[i * 2]);
      const std::uint32_t to_push = low | (high << 16);
      m_vFragment->push_back(to_push);
      sum += to_push;
    } 
    setHeader16(ResultsDim2,m_DigitsEncode.size()/2);
    setHeader16(ResultsOff2,18+size_of_block);
  } // End of check for format
 
  // Need to add header to check sum
  for(size_t ii=0;ii<endtag/2;ii++){
    sum+=((*m_vFragment)[ii]);
  }
  // Three final magic words
  m_vFragment->push_back(0x0); // For the moment
  m_vFragment->push_back(0x12345678); // For the moment
  m_vFragment->push_back(sum& 0x7fffffff);
   
  setHeader32(NWTot,m_vFragment->size());
}
 
 
void  LArRodBlockPhysicsV5::concatinateFEBs()
{
  FEBMAPTYPE::const_iterator feb_it_b=m_mFebBlocks.begin();
  FEBMAPTYPE::const_iterator feb_it_e=m_mFebBlocks.end();
  FEBMAPTYPE::const_iterator feb_it;
  for (feb_it=feb_it_b;feb_it!=feb_it_e;++feb_it) {
    if (feb_it!=feb_it_b) //Not first Feb
      m_pRODblock->resize( m_pRODblock->size()+m_MiddleHeaderSize);
 
    //Add feb data to rod data block
    m_pRODblock->insert (m_pRODblock->end(),
                         feb_it->second.begin(), feb_it->second.end());
  } //end for feb_it
 
  m_mFebBlocks.clear();
}
 
//Sort functions & ordering relation:
template<class RAWDATA>
bool LArRodBlockPhysicsV5::operator () 
  (const RAWDATA* ch1, const RAWDATA* ch2) const
{
  HWIdentifier id1 = ch1->channelID(); 
  HWIdentifier id2 = ch2->channelID(); 
  
  HWIdentifier febId1= m_onlineHelper->feb_Id(id1); 
  HWIdentifier febId2= m_onlineHelper->feb_Id(id2); 
 
  if(febId1 == febId2 ){
    int cId1 =  m_onlineHelper->channel(id1); 
    int cId2 =  m_onlineHelper->channel(id2);
    return FebToRodChannel(cId1) < FebToRodChannel(cId2);
  } 
 
  return febId1 < febId2 ; 
}