CaloWaveformDigiAlg.cxx

#include "CaloWaveformDigiAlg.h"

#include "Identifier/Identifier.h"

#include "FaserCaloSimEvent/CaloHitIdHelper.h"

#include <map>
#include <utility>

CaloWaveformDigiAlg::CaloWaveformDigiAlg(const std::string& name, 
					 ISvcLocator* pSvcLocator)
  : AthReentrantAlgorithm(name, pSvcLocator)
{ 

}

StatusCode 
CaloWaveformDigiAlg::initialize() {
  ATH_MSG_INFO(name() << "::initalize()");

  // Initalize tools
  ATH_CHECK( m_digiTool.retrieve() );
  ATH_CHECK( m_mappingTool.retrieve() );
  ATH_CHECK( m_digiCondTool.retrieve() );

  // Set key to read waveform from
  ATH_CHECK( m_caloHitContainerKey.initialize() );

  // Set key to write container
  ATH_CHECK( m_waveformContainerKey.initialize() );

  // Set up  helper
  ATH_CHECK(detStore()->retrieve(m_ecalID, "EcalID"));

  // Need to set this on first event
  m_kernel = 0;

  return StatusCode::SUCCESS;
}

StatusCode 
CaloWaveformDigiAlg::finalize() {
  ATH_MSG_INFO(name() << "::finalize()");

  if (m_kernel)
    delete m_kernel;

  return StatusCode::SUCCESS;
}

StatusCode 
CaloWaveformDigiAlg::execute(const EventContext& ctx) const {
  ATH_MSG_DEBUG("Executing");
  ATH_MSG_DEBUG("Run: " << ctx.eventID().run_number() << " Event: " << ctx.eventID().event_number());

  // Find the input HIT collection
  SG::ReadHandle<CaloHitCollection> caloHitHandle(m_caloHitContainerKey, ctx);

  ATH_CHECK( caloHitHandle.isValid() );
  ATH_MSG_DEBUG("Found ReadHandle for CaloHitCollection " << m_caloHitContainerKey);

  // Find the output waveform container
  SG::WriteHandle<RawWaveformContainer> waveformContainerHandle(m_waveformContainerKey, ctx);
  ATH_CHECK( waveformContainerHandle.record( std::make_unique<RawWaveformContainer>()) );

  ATH_MSG_DEBUG("WaveformsContainer '" << waveformContainerHandle.name() << "' initialized");

  if (caloHitHandle->size() == 0) {
    ATH_MSG_DEBUG("CaloHitCollection found with zero length!");
    return StatusCode::SUCCESS;
  }

  if (!m_kernel) {

    ATH_MSG_INFO(name() << ": initialize waveform digitization kernel");
    ATH_MSG_INFO(" Norm:  " << m_digiCondTool->cb_norm(ctx));
    ATH_MSG_INFO(" Mean:  " << m_digiCondTool->cb_mean(ctx));
    ATH_MSG_INFO(" Sigma: " << m_digiCondTool->cb_sigma(ctx));
    ATH_MSG_INFO(" Alpha: " << m_digiCondTool->cb_alpha(ctx));
    ATH_MSG_INFO(" N:     " << m_digiCondTool->cb_n(ctx));

    // Set up waveform shape
    m_kernel = new TF1("PDF", "[4] * ROOT::Math::crystalball_pdf(x, [0],[1],[2],[3])", 0, 1200);

    m_kernel->SetParameter(0, m_digiCondTool->cb_alpha(ctx));
    m_kernel->SetParameter(1, m_digiCondTool->cb_n(ctx));
    m_kernel->SetParameter(2, m_digiCondTool->cb_sigma(ctx));
    m_kernel->SetParameter(3, m_digiCondTool->cb_mean(ctx));
    m_kernel->SetParameter(4, m_digiCondTool->cb_norm(ctx));
    
    // Pre-evaluate time kernel for each bin
    m_timekernel = m_digiTool->evaluate_timekernel(m_kernel);

    // Also save the baseline parameters
    m_base_mean = m_digiCondTool->base_mean(ctx); 
    m_base_rms  = m_digiCondTool->base_rms(ctx); 
  }

  // Create structure to store pulse for each channel
  std::map<Identifier, std::vector<uint16_t>> waveforms = m_digiTool->create_waveform_map(m_ecalID);

  // Sum energy for each channel (i.e. identifier)
  std::map<unsigned int, float> esum;  
  for (const auto& hit : *caloHitHandle) { 
    esum[hit.identify()] += hit.energyLoss();
  }

  // Loop over time samples
  for (const auto& tk : m_timekernel) {
    std::map<unsigned int, float> counts;    

    // Convolve hit energy with evaluated kernel and sum for each hit id (i.e. channel)
    //for (const auto& hit : *caloHitHandle) { 
    //  counts[hit.identify()] += tk * hit.energyLoss();
    //}

    // Convolve summed energy with evaluated kernel for each hit id (i.e. channel)
    for (const auto& e : esum) {
      counts[e.first] = tk * e.second;
    }

    // Subtract count from basleine and add result to correct waveform vector
    for (const auto& c : counts) {

      unsigned int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
      int value = baseline - c.second;

      if (value < 0) {
	ATH_MSG_WARNING("Found pulse " << c.second << " larger than baseline " << c.first);
	value = 0; // Protect against scaling signal above baseline
      }

      // Convert hit id to Identifier and store 
      Identifier id = CaloHitIdHelper::GetHelper()->getIdentifier(c.first);
      waveforms[id].push_back(value);
    }
  }


  // This is a bit of a hack to make sure all waveforms have
  // at least baseline entries.  Should really add this to the 
  // logic above
  for (const auto& w : waveforms) {
    if (w.second.size() > 0) continue;

    // Waveform was empty, fill with baseline
    int channel = m_mappingTool->getChannelMapping(w.first);
    ATH_MSG_DEBUG("Writing baseline into empty waveform in channel "<< channel);
    int i = m_digiTool->nsamples();
    while(i--) {  // Use while to avoid unused variable warning with for
      int baseline = m_digiTool->generate_baseline(m_base_mean, m_base_rms);
      waveforms[w.first].push_back(baseline);
    }
  }

  //m_chrono->chronoStop("Digit");
  //m_chrono->chronoStart("Write");

  // Loop over wavefrom vectors to make and store waveform
  unsigned int nsamples = m_digiTool->nsamples();
  for (const auto& w : waveforms) {
    RawWaveform* wfm = new RawWaveform();
    wfm->setWaveform(0, w.second);
    wfm->setIdentifier(w.first);
    wfm->setChannel(m_mappingTool->getChannelMapping(w.first));
    wfm->setSamples(nsamples);
    waveformContainerHandle->push_back(wfm);
  }

  //m_chrono->chronoStop("Write");

  ATH_MSG_DEBUG("WaveformsHitContainer " << waveformContainerHandle.name() << "' filled with "<< waveformContainerHandle->size() <<" items");

  return StatusCode::SUCCESS;
}