LArRDOAnalysis.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
 
#include "LArRDOAnalysis.h"
#include "StoreGate/ReadHandle.h"
#include "TTree.h"
#include "TH1F.h"
#include <format>
 
StatusCode LArRDOAnalysis::initialize() {
  ATH_MSG_DEBUG( "Initializing LArRDOAnalysis" );
 
  // This will check that the properties were initialized
  // properly by job configuration.
  ATH_CHECK( m_inputRawChannelKey.initialize(SG::AllowEmpty) );
  ATH_CHECK( m_inputTTL1HADKey.initialize(SG::AllowEmpty) );
  ATH_CHECK( m_inputTTL1EMKey.initialize(SG::AllowEmpty) );
  ATH_CHECK( m_inputDigitKey.initialize(SG::AllowEmpty) );
 
  if(m_doNtuple){
    m_tree = new TTree(m_ntupleTreeName.value().c_str(), "LArRDOAna");
    std::string fullNtupleName =std::format("{}{}{}",m_ntupleFileName.value(), m_ntupleDirName.value(), 
                                                    m_ntupleTreeName.value());
    ATH_CHECK(histSvc()->regTree(fullNtupleName, m_tree));
    m_tree->Branch("larID", &m_larID);
    m_tree->Branch("energy", &m_energy);
    m_tree->Branch("time", &m_time);
    m_tree->Branch("qual", &m_qual);
    m_tree->Branch("prov", &m_prov);
    m_tree->Branch("gain", &m_gain);
    m_tree->Branch("hadOnID", &m_hadOnID);
    m_tree->Branch("hadOffID", &m_hadOffID);
    m_tree->Branch("hadSamples", &m_hadSamples);
    m_tree->Branch("emOnID", &m_emOnID);
    m_tree->Branch("emOffID", &m_emOffID);
    m_tree->Branch("emSamples", &m_emSamples);
    m_tree->Branch("digiID", &m_digiID);
    m_tree->Branch("digiGain", &m_digiGain);
    m_tree->Branch("digiSamples", &m_digiSamples);
  }
  
  m_h_larID = new TH1F("h_larID", "LAr ID", 100, 0, 5e18);
  m_h_larID->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_larID->GetName(), m_h_larID));
 
  m_h_energy = new TH1F("h_energy", "LAr energy", 100, -1e5, 5e5);
  m_h_energy->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_energy->GetName(), m_h_energy));
 
  m_h_time = new TH1F("h_time", "LAr time", 100, -1e7, 1e7);
  m_h_time->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_time->GetName(), m_h_time));
 
  m_h_qual = new TH1F("h_qual", "LAr quality", 100, 0, 70000);
  m_h_qual->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_qual->GetName(), m_h_qual));
 
  m_h_prov = new TH1F("h_prov", "LAr provenance", 100, 0, 9000);
  m_h_prov->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_prov->GetName(), m_h_prov));
 
  m_h_gain = new TH1F("h_gain", "LAr gain", 100, 0, 5);
  m_h_gain->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_gain->GetName(), m_h_gain));
 
  m_h_hadOnID = new TH1F("h_hadOnID", "Had LAr TTL1 online ID", 100, 0, 3e19);
  m_h_hadOnID->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hadOnID->GetName(), m_h_hadOnID));
 
  m_h_hadOffID = new TH1F("h_hadOffID", "Had LAr TTL1 offline ID", 100, 0, 3e19);
  m_h_hadOffID->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hadOffID->GetName(), m_h_hadOffID));
 
  m_h_hadSamples = new TH1F("h_hadSamples", "Had LAr TTL1 sample values", 100, -15000, 35000);
  m_h_hadSamples->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_hadSamples->GetName(), m_h_hadSamples));
 
  m_h_emOnID = new TH1F("h_emOnID", "EM LAr TTL1 online ID", 100, 0, 3e19);
  m_h_emOnID->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_emOnID->GetName(), m_h_emOnID));
 
  m_h_emOffID = new TH1F("h_emOffID", "EM LAr TTL1 offline ID", 100, 0, 3e19);
  m_h_emOffID->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_emOffID->GetName(), m_h_emOffID));
 
  m_h_emSamples = new TH1F("h_emSamples", "EM LAr TTL1 sample values", 100, -1e5, 3.5e5);
  m_h_emSamples->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_emSamples->GetName(), m_h_emSamples));
 
  m_h_digiID = new TH1F("h_digiID", "LAr digit ID", 100, 0, 5e18);
  m_h_digiID->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_digiID->GetName(), m_h_digiID));
 
  m_h_digiGain = new TH1F("h_digiGain", "LAr digit gain", 100, 0, 5);
  m_h_digiGain->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_digiGain->GetName(), m_h_digiGain));
 
  m_h_digiSamples = new TH1F("h_digiSamples", "LAr digit sample values", 100, 0, 5000);
  m_h_digiSamples->StatOverflows();
  ATH_CHECK(histSvc()->regHist(m_path + m_h_digiSamples->GetName(), m_h_digiSamples));
 
  return StatusCode::SUCCESS;
}
 
StatusCode LArRDOAnalysis::execute() {
  ATH_MSG_DEBUG( "In LArRDOAnalysis::execute()" );
  
  if(m_tree){
    m_larID->clear();
    m_energy->clear();
    m_time->clear();
    m_qual->clear();
    m_prov->clear();
    m_gain->clear();
    m_hadOnID->clear();
    m_hadOffID->clear();
    m_hadSamples->clear();
    m_emOnID->clear();
    m_emOffID->clear();
    m_emSamples->clear();
    m_digiID->clear();
    m_digiGain->clear();
    m_digiSamples->clear();
  }
 
  const EventContext& ctx{Gaudi::Hive::currentContext()};
  const LArRawChannelContainer* p_larRawCont{nullptr};
  ATH_CHECK(SG::get(p_larRawCont, m_inputRawChannelKey, ctx));
 
  const LArTTL1Container* p_larTTL1Cont_had{nullptr};
  ATH_CHECK(SG::get(p_larTTL1Cont_had, m_inputTTL1HADKey, ctx));
 
  const LArTTL1Container* p_larTTL1Cont_em{nullptr};
  ATH_CHECK(SG::get(p_larTTL1Cont_em, m_inputTTL1EMKey, ctx));
 
  const LArDigitContainer* p_larDigiCont{nullptr};
  ATH_CHECK(SG::get(p_larDigiCont, m_inputDigitKey, ctx));
  
 
  if (p_larRawCont) {
      // loop over LAr raw channels container
      LArRawChannelContainer::const_iterator lar_itr(p_larRawCont->begin());
      const LArRawChannelContainer::const_iterator lar_end(p_larRawCont->end());
      for ( ; lar_itr != lar_end; ++lar_itr ) {
          const HWIdentifier larID(lar_itr->identify());
          const int rawEnergy(lar_itr->energy());
          const int rawTime(lar_itr->time());
          const uint16_t rawQual(lar_itr->quality());
          const uint16_t rawProv(lar_itr->provenance());
          CaloGain::CaloGain larGain(lar_itr->gain());
 
          const unsigned long long larID_int = larID.get_compact();
          const int larGain_int = (int)larGain;
          if(m_tree){
            m_larID->push_back(larID_int);
            m_energy->push_back(rawEnergy);
            m_time->push_back(rawTime);
            m_qual->push_back(rawQual);
            m_prov->push_back(rawProv);
            m_gain->push_back(larGain_int);
          }
 
          m_h_larID->Fill(larID_int);
          m_h_energy->Fill(rawEnergy);
          m_h_time->Fill(rawTime);
          m_h_qual->Fill(rawQual);
          m_h_prov->Fill(rawProv);
          m_h_gain->Fill(larGain_int);
      }
  }
  if (p_larTTL1Cont_had) {
    LArTTL1Container::const_iterator ttl1Had_itr(p_larTTL1Cont_had->begin());
    const LArTTL1Container::const_iterator ttl1Had_end(p_larTTL1Cont_had->end());
    for ( ; ttl1Had_itr != ttl1Had_end; ++ttl1Had_itr ) {
          const HWIdentifier& hadOnID((*ttl1Had_itr)->ttOnlineID());
          const Identifier& hadOffID((*ttl1Had_itr)->ttOfflineID());
          const std::vector<float>& hadSamples = (*ttl1Had_itr)->samples();
 
          const unsigned long long hadOnID_int = hadOnID.get_compact();
          const unsigned long long hadOffID_int = hadOffID.get_compact();
          if(m_tree) m_hadOnID->push_back(hadOnID_int);
          if(m_tree) m_hadOffID->push_back(hadOffID_int);
          for (std::vector<float>::size_type i = 0; i != hadSamples.size(); ++i) {
            if(m_tree) m_hadSamples->push_back(hadSamples.at(i));
            m_h_hadSamples->Fill(hadSamples.at(i));
          }
 
          m_h_hadOnID->Fill(hadOnID_int);
          m_h_hadOffID->Fill(hadOffID_int);
    }
  }
 
  if (p_larTTL1Cont_em) {
        LArTTL1Container::const_iterator ttl1EM_itr(p_larTTL1Cont_em->begin());
        const LArTTL1Container::const_iterator ttl1EM_end(p_larTTL1Cont_em->end());
        for ( ; ttl1EM_itr != ttl1EM_end; ++ttl1EM_itr ) {
          const HWIdentifier& emOnID((*ttl1EM_itr)->ttOnlineID());
          const Identifier& emOffID((*ttl1EM_itr)->ttOfflineID());
          const std::vector<float>& emSamples((*ttl1EM_itr)->samples());
 
          const unsigned long long emOnID_int = emOnID.get_compact();
          const unsigned long long emOffID_int = emOffID.get_compact();
          if(m_tree) m_emOnID->push_back(emOnID_int);
          if(m_tree) m_emOffID->push_back(emOffID_int);
          for (std::vector<float>::size_type j = 0; j != emSamples.size(); ++j) {
            if(m_tree) m_emSamples->push_back(emSamples.at(j));
            m_h_emSamples->Fill(emSamples.at(j));
          }
 
          m_h_emOnID->Fill(emOnID_int);
          m_h_emOffID->Fill(emOffID_int);
        }
  }
 
 
  // LAr Digits
  if (p_larDigiCont) {
    LArDigitContainer::const_iterator digi_itr(p_larDigiCont->begin());
    const LArDigitContainer::const_iterator digi_end(p_larDigiCont->end());
    for ( ; digi_itr != digi_end; ++digi_itr ) {
      const HWIdentifier& digiID((*digi_itr)->hardwareID());
      CaloGain::CaloGain digiGain((*digi_itr)->gain());
      const std::vector<short>& digiSamples = (*digi_itr)->samples();
 
      const unsigned long long digiID_int = digiID.get_compact();
      const int digiGain_int = (int)digiGain;
      if(m_tree) m_digiID->push_back(digiID_int);
      if(m_tree) m_digiGain->push_back(digiGain_int);
      for (std::vector<short>::size_type k = 0; k != digiSamples.size(); ++k) {
        if(m_tree) m_digiSamples->push_back(digiSamples.at(k));
        m_h_digiSamples->Fill(digiSamples.at(k));
      }
 
      m_h_digiID->Fill(digiID_int);
      m_h_digiGain->Fill(digiGain_int);
   }
 }
  
 
 if(m_tree and m_doNtuple){ 
  m_tree->Fill();
 }
  
  return StatusCode::SUCCESS;
}