jFexEmulatedTowers.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
//***************************************************************************
//                           jFexEmulatedTowers  -  description
//                              -------------------
//      This reentrant algorithm is meant build jFEX Towers from LAr and Tile
//      Information about SCellContainer objetcs are in:
//          - https://gitlab.cern.ch/atlas/athena/-/blob/22.0/Calorimeter/CaloEvent/CaloEvent/CaloCell.h
//      
//     begin                : 01 11 2022
//     email                : sergi.rodriguez@cern.ch
//***************************************************************************/
 
 
#include "jFexEmulatedTowers.h"
#include "L1CaloFEXSim/jFEXCompression.h"
#include "StoreGate/WriteDecorHandle.h"
#include "StoreGate/ReadDecorHandle.h"
 
#include <iostream>
#include <fstream>
#include <sstream>
#include <algorithm>
#include <string>
 
namespace LVL1 {
 
jFexEmulatedTowers::jFexEmulatedTowers(const std::string& name, ISvcLocator* svc) : AthReentrantAlgorithm(name, svc){}
 
StatusCode jFexEmulatedTowers::initialize() {
    
    ATH_MSG_INFO( "Initializing L1CaloFEXAlgos/jFexEmulatedTowers algorithm with name: "<< name());
    ATH_MSG_INFO( "Writting into SG key: "<< m_jTowersWriteKey);
    ATH_MSG_INFO( "SCell masking: "<< m_apply_masking);
    ATH_MSG_INFO( "Thinnig towers: "<< m_doThinning);
    
    ATH_CHECK( m_SCellKey.initialize() );
    ATH_CHECK( m_triggerTowerKey.initialize() );
    ATH_CHECK( m_jTowersWriteKey.initialize() );    
 
    //Reading from CVMFS Fiber mapping
    ATH_CHECK(ReadFibersfromFile(m_FiberMapping));    
    
    //Reading from CVMFS Trigger Tower and their corresponding SCell ID
    ATH_CHECK(ReadSCfromFile(m_jFEX2Scellmapping));
    ATH_CHECK(ReadTilefromFile(m_jFEX2Tilemapping));
 
    // offline energy decoration setup
    // SCell to Cell mapping
    if (!m_CaloCellKey.empty()) {
      ATH_CHECK( m_scellIdTool.retrieve() );
      ATH_CHECK( detStore()->retrieve (m_caloCellIdHelper, "CaloCell_ID") );
    }
 
    // Calo cell container
    ATH_CHECK( m_CaloCellKey.initialize(SG::AllowEmpty) );
 
    // read Tile Tower decoration
    ATH_CHECK( m_readTileOfflineDecorKey.initialize(SG::AllowEmpty) );
 
    // Readout path cell sum Et decoration
    ATH_CHECK( m_caloCellSumETdecorKey.initialize(SG::AllowEmpty) );
 
    // tower energy decorations
    ATH_CHECK( m_jtowerEtMeVdecorKey.initialize(SG::AllowEmpty) );
    ATH_CHECK( m_jtowerEtTimingMeVdecorKey.initialize(SG::AllowEmpty) );
 
    // Check if ApplyTimingCut and ApplyTimingCutAll flags are both enabled
    if (m_applyTimingCut && m_applyTimingCutAll) {
        ATH_MSG_FATAL(
            "ApplyTimingCut and ApplyTimingCutAll flags cannot both be enabled");
        return StatusCode::FAILURE;
    }
 
 
    return StatusCode::SUCCESS;
}
 
StatusCode jFexEmulatedTowers::execute(const EventContext& ctx) const {
 
    //Reading the Scell container
    SG::ReadHandle<CaloCellContainer> ScellContainer(m_SCellKey, ctx);
    if(!ScellContainer.isValid()) {
        ATH_MSG_ERROR("Could not retrieve collection " << ScellContainer.key() );
        return StatusCode::FAILURE;
    }
    
    //Reading the TriggerTower container
    SG::ReadHandle<xAOD::TriggerTowerContainer> triggerTowerContainer(m_triggerTowerKey, ctx);
    if(!triggerTowerContainer.isValid()) {
        ATH_MSG_ERROR("Could not retrieve collection " << triggerTowerContainer.key() );
        return StatusCode::FAILURE;
    }
 
    // Reading the offline Calo cell container
    const CaloCellContainer* caloCells = nullptr;
    if (!m_CaloCellKey.empty()) {
      SG::ReadHandle<CaloCellContainer> caloCellContainer(m_CaloCellKey, ctx);
      if(!caloCellContainer.isValid()) {
        ATH_MSG_ERROR("Could not retrieve collection " << caloCellContainer.key() );
        return StatusCode::FAILURE;
      }
      caloCells = caloCellContainer.cptr();
    }
 
    //WriteHandle for jFEX EDMs
    //---jTower EDM
    SG::WriteHandle<xAOD::jFexTowerContainer> jTowersContainer(m_jTowersWriteKey, ctx);
    ATH_CHECK(jTowersContainer.record(std::make_unique<xAOD::jFexTowerContainer>(), std::make_unique<xAOD::jFexTowerAuxContainer>()));
    ATH_MSG_DEBUG("Recorded jFexEmulatedTower container with key " << jTowersContainer.key());
 
    if(ScellContainer->empty() || triggerTowerContainer->empty() ){
        ATH_MSG_WARNING("Cannot fill jTowers here, at least one container is empty. ScellContainer.size="<<ScellContainer->size() << " or triggerTowerContainer.size=" << triggerTowerContainer->size() );
        return StatusCode::SUCCESS;
    }
 
    std::optional<SG::ReadDecorHandle<xAOD::TriggerTowerContainer, std::vector<float>>> cellEtByLayer_handle;
    if (!m_CaloCellKey.empty()) {
      cellEtByLayer_handle.emplace(m_readTileOfflineDecorKey, ctx);
    }
    std::optional<SG::WriteDecorHandle<xAOD::jFexTowerContainer, float>> jTowerCaloCellSumEt;
    if (!m_caloCellSumETdecorKey.empty()) {
      jTowerCaloCellSumEt.emplace(m_caloCellSumETdecorKey, ctx);
    }
    std::optional<SG::WriteDecorHandle<xAOD::jFexTowerContainer, int>> jTowerEtMeV;
    if (!m_jtowerEtMeVdecorKey.empty()) {
      jTowerEtMeV.emplace(m_jtowerEtMeVdecorKey , ctx);
    }
    std::optional<SG::WriteDecorHandle<xAOD::jFexTowerContainer, int>> jTowerEtTimingMeV;
    if (!m_jtowerEtTimingMeVdecorKey.empty()) {
      jTowerEtTimingMeV.emplace(m_jtowerEtTimingMeVdecorKey , ctx);
    }
 
    // building Scell ID pointers
    std::unordered_map< uint64_t, const CaloCell*> map_ScellID2ptr;
    map_ScellID2ptr.reserve(ScellContainer->size());
    
    for(const CaloCell* scell : *ScellContainer){
        const uint64_t ID = scell->ID().get_compact();
        map_ScellID2ptr[ID] = scell;
    }
    
    // building Tile ID pointers
    std::unordered_map< uint32_t, const xAOD::TriggerTower*> map_TileID2ptr;
    
    for(const xAOD::TriggerTower* tower : *triggerTowerContainer){
        
        // keeping just tile information
        if(std::abs(tower->eta())>1.5 || tower->sampling()!=1) continue;
        
        map_TileID2ptr[tower->coolId()]=tower;
    }
    
    
    for( const auto& [key, element] : m_Firm2Tower_map){
        
        unsigned int jfex    = (key >> 16) & 0xf;
        unsigned int fpga    = (key >> 12) & 0xf;
        unsigned int channel = (key >> 4 ) & 0xff;
        unsigned int tower   = (key >> 0 ) & 0xf; 
               
        const auto [f_IDSimulation, eta, phi, f_source, f_iEta, f_iPhi] = element;
        
        //elements that need to be casted
        unsigned int IDSimulation = static_cast<int>(f_IDSimulation);
        unsigned int source       = static_cast<int>(f_source);
        unsigned int iEta         = static_cast<int>(f_iEta);
        unsigned int iPhi         = static_cast<int>(f_iPhi);
        
        uint16_t Total_Et_encoded = 0;
        uint16_t Total_Et_Timing_encoded = 0;
        char jTower_sat = 0;
 
    float caloCellSumET(0); // Readout Et for Tower
    int Total_Et_decoded(0);
    int Total_Et_Timing_decoded(0);
        if( source != 1 ){
            
            const std::unordered_map< uint32_t, std::vector<uint64_t> > * ptr_TTower2Cells;
            
            //HAD layer for HEC, FCAL2 and FCAL3
            if(source == 3 or source > 4){
                ptr_TTower2Cells = &m_map_TTower2SCellsHAD;
            } 
            else{
                ptr_TTower2Cells = &m_map_TTower2SCellsEM;
            }  
            
            //check that the jFEX Tower ID exists in the map
            auto it_TTower2SCells = (*ptr_TTower2Cells).find(IDSimulation);
            if(it_TTower2SCells == (*ptr_TTower2Cells).end()) {
                ATH_MSG_ERROR("jFEX ID: "<<IDSimulation<< " not found on map m_map_TTower2SCellsEM/HAD");
                return StatusCode::FAILURE;
            }
        
            int Total_Et = 0;
            int Total_Timing_Et = 0;
            float Total_Et_float = 0;
            bool invalid = m_apply_masking&&m_isDATA; // the isDATA is because there is no concept of invalid supercell in MC (the provenance bit is actually used for BCID in MC), so can never have an invalid jTower. 
            bool masked = m_apply_masking;
            for (auto const& SCellID : it_TTower2SCells->second ) {
                //check that the SCell Identifier exists in the map
                auto it_ScellID2ptr = map_ScellID2ptr.find(SCellID);
                if(it_ScellID2ptr == map_ScellID2ptr.end()) {
                    if(m_isDATA) ATH_MSG_DEBUG("Scell ID: 0x"<<std::hex<< (SCellID >> 32) <<std::dec<< " not found in the CaloCell Container, skipping");
                    // this is equivalent to treat the scell as a masked input, since masking takes precedence over invalidity
                    continue;
                }
 
                const CaloCell* myCell = it_ScellID2ptr->second;
                int val =  std::round(myCell->energy()/(12.5*std::cosh(myCell->eta()))); // 12.5 is b.c. energy is in units of 12.5MeV per count
                bool isMasked = m_apply_masking ? (myCell->provenance()&0x80) : false;
                bool isInvalid = (m_apply_masking&&m_isDATA) ? (myCell->provenance()&0x40) : false;
                bool isSaturated = (m_isDATA) ? myCell->quality() : false; // saturation algorithm not implemented in MC yet
                bool passTiming = ( myCell->provenance() & 0x200 );
 
                invalid &= isInvalid;
                masked &= isMasked;
        if (!isMasked) {
          jTower_sat |= isSaturated;
        }
 
                if( isMasked ) {
                    //if masked then Et = 0
                    val = 0;
                    //countMasked++;
                } else if( isInvalid ) {
                    val = 0;
                }
                // Apply timing cut to MC if the respective flags are set
                if (!m_isDATA && m_applyTimingCut && !passTiming && source != 3) {
                    // Applying timing cut everywhere but HEC if m_applyTimingCut set
                    val = 0;
                } else if (!m_isDATA && m_applyTimingCutAll && !passTiming) {
                    // Applying timing cut on all SCells if m_applyTimingCutAll set
                    val = 0;
                }
                
                Total_Et += val;
                if(val!=0) Total_Et_float += myCell->et();
 
        // tau timing cut
        Total_Timing_Et += (passTiming) ? val : 0;
 
        // Readout path energy sum
        if ( !m_CaloCellKey.empty() ) {
          const std::vector<Identifier> caloCellIds = m_scellIdTool->superCellToOfflineID(myCell->ID());
          float caloCellEt(0.);  // offline Et for the SCell
          //use findCell function of CaloCellContainer, which takes an identifier hash
          for(auto& caloCellId : caloCellIds) {
            const CaloCell* caloCell = caloCells->findCell(m_caloCellIdHelper->calo_cell_hash(caloCellId));
            if(!caloCell) { ATH_MSG_WARNING("Could not find cell"); continue; }
            caloCellEt += caloCell->e()*caloCell->sinTh();
          }
          caloCellSumET += caloCellEt;
        }
                
            } // SCell loop
 
            // now must convert Total_Et int value into fex value: multi-level encoding
            if(masked) {
                Total_Et_encoded = 0; // no data
        Total_Et_Timing_encoded = 0;
            } else if(invalid) {
                Total_Et_encoded = 4095; // invalid
                Total_Et_Timing_encoded = 4095;
            } else {
                Total_Et_encoded = jFEXCompression::Compress( Total_Et*12.5);
                Total_Et_Timing_encoded = jFEXCompression::Compress( Total_Timing_Et*12.5);
            }
 
            // leaving this commented while outstanding questions above about treatment of supercells
            // using floating point for MC until determine correct procedure for MC values re invalid/masking
            if(!m_isDATA) Total_Et_encoded = jFEXCompression::Compress( Total_Et_float, masked );
 
        Total_Et_decoded = jFEXCompression::Expand( Total_Et_encoded );
        Total_Et_Timing_decoded = jFEXCompression::Expand( Total_Et_Timing_encoded );
        }
        else{
            
            //check that the jFEX Tower ID exists in the map
            auto it_TTower2Tile = m_map_TTower2Tile.find(IDSimulation);
            if(it_TTower2Tile == m_map_TTower2Tile.end()) {
                ATH_MSG_ERROR("ID: "<<IDSimulation<< " not found on map m_map_TTower2Tile");
                return StatusCode::FAILURE;
            }
            
            uint32_t TileID = std::get<0>( it_TTower2Tile->second );
            
            //check that the Tile Identifier exists in the map
            auto it_TileID2ptr = map_TileID2ptr.find(TileID);
            if(it_TileID2ptr == map_TileID2ptr.end()) {
                if(m_isDATA) ATH_MSG_WARNING("Tile cool ID: "<<TileID<< " not found in the xAOD::TriggerTower, skipping");
                continue;
            }
 
        Total_Et_encoded = (it_TileID2ptr->second)->cpET();
 
        Total_Et_decoded = Total_Et_encoded * 500;
        Total_Et_Timing_decoded = Total_Et_encoded * 500;
 
        // get the offline reconstructed energy from xAODTriggerTower decoration
        if ( !m_CaloCellKey.empty() ) {
          if ( cellEtByLayer_handle->isAvailable() ) {
        const std::vector<float>& cellEtByLayer = (*cellEtByLayer_handle)(*it_TileID2ptr->second);
        // Stored in GeV, see TrigT1CaloCalibTools/src/L1CaloCells2TriggerTowers.cxx
        caloCellSumET = 1000.0 * std::accumulate(cellEtByLayer.begin(), cellEtByLayer.end(), 0.0);
          }
        }
        } // end of Tile
        std::vector<uint16_t> vtower_ET;
        vtower_ET.push_back(Total_Et_encoded);
        
        std::vector<char> vtower_SAT;
        
        //Needs to be updated with Saturation flag from LAr CaloCell container, not ready yet!
        vtower_SAT.push_back(jTower_sat);  
        
        jTowersContainer->push_back( std::make_unique<xAOD::jFexTower>() );
        jTowersContainer->back()->initialize(eta, phi, iEta, iPhi, IDSimulation, source, vtower_ET, jfex, fpga, channel, tower, vtower_SAT ); 
        
        if( m_doThinning && !jTowersContainer->back()->isCore() ){
            jTowersContainer->pop_back();
        } else {
      // decorate with offline energy
      if (!m_caloCellSumETdecorKey.empty()) {
        (*jTowerCaloCellSumEt) (*jTowersContainer->back()) = caloCellSumET;
      }
      // decorate tower energies in MeV without/with timing cut
      if (!m_jtowerEtMeVdecorKey.empty()) {
        (*jTowerEtMeV) (*jTowersContainer->back()) = Total_Et_decoded;
      }
      if (!m_jtowerEtTimingMeVdecorKey.empty()) {
        (*jTowerEtTimingMeV) (*jTowersContainer->back()) = Total_Et_Timing_decoded;
      }
    }
    } // firmware tower map
    
    // Return gracefully
    return StatusCode::SUCCESS;
}
 
 
StatusCode jFexEmulatedTowers::ReadFibersfromFile(const std::string & fileName){
 
    //opening file with ifstream
    std::ifstream file(fileName);
    
    if ( !file.is_open() ){
        ATH_MSG_ERROR("Could not open file:" << fileName);
        return StatusCode::FAILURE;
    }
    
    std::string line;
    //loading the mapping information
    while ( std::getline (file, line) ) {
 
        //removing the header of the file (it is just information!)
        if(line[0] == '#') continue;
        
        //Splitting line in different substrings
        std::stringstream oneLine(line);
        
        //reading elements
        std::vector<float> elements;
    elements.reserve(10);
        std::string element;
        while(std::getline(oneLine, element, ' '))
        {
            elements.push_back(std::stof(element));
        }
        
        // It should have 10 elements
        // ordered as:  jfex fpga channel towerNr source globalEtaIndex globalPhiIndex IDSimulation eta phi
        if(elements.size() != 10){
            ATH_MSG_ERROR("Unexpected number of elemennts (10 expected) in file: "<< fileName);
            return StatusCode::FAILURE;
        }
        // building array of  <IDSimulation, eta, phi, source, iEta, iPhi>
        std::array<float,6> aux_arr{ {elements.at(7),elements.at(8),elements.at(9),elements.at(4),elements.at(5),elements.at(6)} };
        
        //filling the map with the hash given by mapIndex()
        m_Firm2Tower_map[ mapIndex(elements.at(0),elements.at(1),elements.at(2),elements.at(3)) ] = aux_arr;
        
    }
    file.close();
 
    return StatusCode::SUCCESS;
}
 
 
constexpr unsigned int jFexEmulatedTowers::mapIndex(unsigned int jfex, unsigned int fpga, unsigned int channel, unsigned int tower) {
  // values from hardware: jfex=[0,5] 4 bits, fpga=[0,3] 4 bits, channel=[0,59] 8 bits, tower=[0,15] 4 bits
  return (jfex << 16) | (fpga << 12) | (channel << 4) | tower;
}
 
StatusCode  jFexEmulatedTowers::ReadSCfromFile(const std::string& fileName){
 
    //opening file with ifstream
    std::ifstream file(fileName);
    
    if ( !file.is_open() ){
        ATH_MSG_ERROR("Could not open file:" << fileName);
        return StatusCode::FAILURE;
    }
    
    std::string line;
    //loading the mapping information into an unordered_map <Fex Tower ID, vector of SCell IDs>
    while ( std::getline (file, line) ) {
        //removing the header of the file (it is just information!)
        if(line[0] == '#') continue;
 
    std::vector<uint64_t> SCellvectorEM;
    std::vector<uint64_t> SCellvectorHAD;
 
        //Splitting line in different substrings
        std::stringstream oneSCellID(line);
        
        //reading elements
        std::string substr = "";
        int TTID = 0;
        int elem = 0;
        
        while(std::getline(oneSCellID, substr, ' '))
        {
            ++elem;
            if(elem == 1){
               TTID =  std::stoi(substr);
            }
            else{
                //Check if it looks like a SCell Identifier
                if(isBadSCellID(substr)){
                    return StatusCode::FAILURE;
                }
                
                // converts hex number to unsigned long long int
                uint64_t scid_uint64 = std::strtoull(substr.c_str(), nullptr, 0);
                
                //empty slots are filled with 0xffffffffffffffff
                if(scid_uint64 == 0xffffffffffffffff) continue;
                
                //from element from 2 to 13 are EM SCells, element 14 is a HAD SCell
                if(elem<14) SCellvectorEM.push_back(scid_uint64);
                else        SCellvectorHAD.push_back(scid_uint64);          
            }
        }        
        
        m_map_TTower2SCellsEM[TTID] = std::move(SCellvectorEM);
        m_map_TTower2SCellsHAD[TTID] = std::move(SCellvectorHAD);
        
    }
    file.close();
 
    return StatusCode::SUCCESS;
}
 
bool jFexEmulatedTowers::isBadSCellID(const std::string& ID) const{
    
    // does it start with "0x"?, if so then is a GOOD SCell ID!
    if (ID.find("0x") == std::string::npos) {
        ATH_MSG_ERROR("Invalid SuperCell ID " << ID << ". Expecting hexadecimal number on the mapping file");
        return true;
    }
    return false;
}
 
 
 
 
StatusCode  jFexEmulatedTowers::ReadTilefromFile(const std::string& fileName){
    
    //opening file with ifstream
    std::ifstream file(fileName);
    
    if ( !file.is_open() ){
        ATH_MSG_ERROR("Could not open file:" << fileName);
        return StatusCode::FAILURE;
    }
    
    std::string line;
    //loading the mapping information into an unordered_map <Fex Tower ID, vector of SCell IDs>
    while ( std::getline (file, line) ) {
 
        //removing the header of the file (it is just information!)
        if(line[0] == '#') continue;
        
        //Splitting line in different substrings
        std::stringstream oneLine(line);
        
        std::vector<std::string> elements;
        std::string element = "";
        
        while(std::getline(oneLine, element, ' ')){
            elements.push_back(element);
        }
        
        if(elements.size() != 4){
            ATH_MSG_ERROR("Invalid number of element in " << line << ". Expecting 4 elements {jFexID, TileID, eta, phi}");
            return StatusCode::FAILURE;
        }
        
        uint32_t jFexID = std::stoi( elements.at(0) );
        uint32_t TileID = std::stoi( elements.at(1) );
        float eta       = std::stof( elements.at(2) );
        float phi       = std::stof( elements.at(3) );
        
        m_map_TTower2Tile[jFexID] = {TileID,eta,phi};
        
    }
    file.close();
    
    return StatusCode::SUCCESS;
}
 
 
 
 
 
}