HGCalTriggerGeometryHexImp1.cc

#include "FWCore/ParameterSet/interface/FileInPath.h"

#include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryGenericMapping.h"
#include "DataFormats/ForwardDetId/interface/HGCalDetId.h"

#include <vector>
#include <iostream>
#include <fstream>

const int HGCalDetId::kHGCalCellMask;

class HGCalTriggerGeometryHexImp1 : public HGCalTriggerGeometryGenericMapping
{
    public:
        HGCalTriggerGeometryHexImp1(const edm::ParameterSet& conf);

        void initialize(const edm::ESHandle<CaloGeometry>& ) final;
        void initialize(const edm::ESHandle<HGCalGeometry>&,
                const edm::ESHandle<HGCalGeometry>&,
                const edm::ESHandle<HGCalGeometry>&) final;

    private:
        edm::FileInPath l1tCellsMapping_;
        edm::FileInPath l1tModulesMapping_;

        void fillMaps();
        void buildTriggerCellsAndModules();
};


/*****************************************************************/
HGCalTriggerGeometryHexImp1::HGCalTriggerGeometryHexImp1(const edm::ParameterSet& conf):
    HGCalTriggerGeometryGenericMapping(conf),
    l1tCellsMapping_(conf.getParameter<edm::FileInPath>("L1TCellsMapping")),
    l1tModulesMapping_(conf.getParameter<edm::FileInPath>("L1TModulesMapping"))
/*****************************************************************/
{
}

/*****************************************************************/
void HGCalTriggerGeometryHexImp1::initialize(const edm::ESHandle<CaloGeometry>& calo_geometry)
/*****************************************************************/
{
    edm::LogWarning("HGCalTriggerGeometry") << "WARNING: This HGCal trigger geometry is incomplete.\n"\
                                            << "WARNING: There is no neighbor information.\n";

    setCaloGeometry(calo_geometry);
    fillMaps();
    buildTriggerCellsAndModules();

}

/*****************************************************************/
 void HGCalTriggerGeometryHexImp1::initialize(const edm::ESHandle<HGCalGeometry>& hgc_ee_geometry,
        const edm::ESHandle<HGCalGeometry>& hgc_hsi_geometry,
        const edm::ESHandle<HGCalGeometry>& hgc_hsc_geometry
        )
/*****************************************************************/
{
    throw cms::Exception("BadGeometry")
        << "HGCalTriggerGeometryHexImp1 geometry cannot be initialized with the V9 HGCAL geometry";
}



/*****************************************************************/
void HGCalTriggerGeometryHexImp1::fillMaps()
/*****************************************************************/
{
    //
    // read module mapping file
    std::unordered_map<short, short> wafer_to_module_ee;
    std::unordered_map<short, short> wafer_to_module_fh;
    std::unordered_map<short, std::map<short,short>> module_to_wafers_ee;
    std::unordered_map<short, std::map<short,short>> module_to_wafers_fh;
    std::ifstream l1tModulesMappingStream(l1tModulesMapping_.fullPath());
    if(!l1tModulesMappingStream.is_open()) edm::LogError("HGCalTriggerGeometry") << "Cannot open L1TModulesMapping file\n";
    short subdet  = 0;
    short wafer   = 0;
    short module  = 0;
    for(; l1tModulesMappingStream>>subdet>>wafer>>module; )
    {
        if(subdet==3)
        {
            wafer_to_module_ee.emplace(wafer,module);
            auto itr_insert = module_to_wafers_ee.emplace(module, std::map<short,short>());
            itr_insert.first->second.emplace(wafer, 0);
        }
        else if(subdet==4)
        {
            wafer_to_module_fh.emplace(wafer,module);
            auto itr_insert = module_to_wafers_fh.emplace(module, std::map<short,short>());
            itr_insert.first->second.emplace(wafer, 0);
        }
        else edm::LogWarning("HGCalTriggerGeometry") << "Unsupported subdetector number ("<<subdet<<") in L1TModulesMapping file\n";
    }
    if(!l1tModulesMappingStream.eof()) edm::LogWarning("HGCalTriggerGeometry") << "Error reading L1TModulesMapping '"<<wafer<<" "<<module<<"' \n";
    l1tModulesMappingStream.close();
    // read trigger cell mapping file
    std::map<std::pair<short,short>, short> cells_to_trigger_cells;
    std::unordered_map<short, short> number_trigger_cells_in_wafers; // the map key is the wafer type
    std::ifstream l1tCellsMappingStream(l1tCellsMapping_.fullPath());
    if(!l1tCellsMappingStream.is_open()) edm::LogError("HGCalTriggerGeometry") << "Cannot open L1TCellsMapping file\n";
    short waferType   = 0;
    short cell        = 0;
    short triggerCell = 0;
    for(; l1tCellsMappingStream>>waferType>>cell>>triggerCell; )
    {
        cells_to_trigger_cells.emplace(std::make_pair((waferType?1:-1),cell), triggerCell);
        auto itr_insert = number_trigger_cells_in_wafers.emplace((waferType?1:-1), 0);
        if(triggerCell+1 > itr_insert.first->second) itr_insert.first->second = triggerCell+1;
    }
    if(!l1tCellsMappingStream.eof()) edm::LogWarning("HGCalTriggerGeometry") << "Error reading L1TCellsMapping'"<<waferType<<" "<<cell<<" "<<triggerCell<<"' \n";
    l1tCellsMappingStream.close();
    // For each wafer compute the trigger cell offset according to the wafer position inside
    // the module. The first wafer will have an offset equal to 0, the second an offset equal to the
    // number of trigger cells in the first wafer, etc.
    short offset = 0;
    for(auto& module_wafers : module_to_wafers_ee)
    {
        offset = 0;
        for(auto& wafer_offset : module_wafers.second)
        {
            wafer_offset.second = offset;
            int wafer_type = (eeTopology().dddConstants().waferTypeT(wafer_offset.first)==1?1:-1);
            offset += number_trigger_cells_in_wafers.at(wafer_type);
        }
    }
    for(auto& module_wafers : module_to_wafers_fh)
    {
        offset = 0;
        for(auto& wafer_offset : module_wafers.second)
        {
            wafer_offset.second = offset;
            int wafer_type = (fhTopology().dddConstants().waferTypeT(wafer_offset.first)==1?1:-1);
            offset += number_trigger_cells_in_wafers.at(wafer_type);
        }
    }
    //
    // Loop over HGC cells
    // EE
    for(const auto& id : eeGeometry()->getValidGeomDetIds())
    {
        if(id.rawId()==0) continue;
        HGCalDetId waferDetId(id); 
        short module      = wafer_to_module_ee[waferDetId.wafer()];
        short triggercell_offset = module_to_wafers_ee.at(module).at(waferDetId.wafer());
        int nCells = eeTopology().dddConstants().numberCellsHexagon(waferDetId.wafer());
        for(int c=0;c<nCells;c++)
        {
            short triggerCellId = cells_to_trigger_cells[std::make_pair(waferDetId.waferType(),c)];
            // Fill cell -> trigger cell mapping
            HGCalDetId cellDetId(ForwardSubdetector(waferDetId.subdetId()), waferDetId.zside(), waferDetId.layer(), waferDetId.waferType(), waferDetId.wafer(), c); 
            // The module id is used instead of the wafer id for the trigger cells
            // Since there are several wafers per module, an offset is applied on the HGCalDetId::cell field
            if(triggercell_offset+triggerCellId >= HGCalDetId::kHGCalCellMask) edm::LogError("HGCalTriggerGeometry") << "Trigger cell id requested with a cell field larger than available in HGCalDetId (" << triggercell_offset+triggerCellId << " >= " << HGCalDetId::kHGCalCellMask << ")\n";
            HGCalDetId triggerCellDetId(ForwardSubdetector(waferDetId.subdetId()), waferDetId.zside(), waferDetId.layer(), waferDetId.waferType(), module, triggercell_offset + triggerCellId);
            cells_to_trigger_cells_.emplace(cellDetId, triggerCellDetId);
            // Fill trigger cell -> module mapping
            HGCalDetId moduleDetId(ForwardSubdetector(waferDetId.subdetId()), waferDetId.zside(), waferDetId.layer(), waferDetId.waferType(), module, HGCalDetId::kHGCalCellMask);
            trigger_cells_to_modules_.emplace(triggerCellDetId, moduleDetId); 
        }
    }
    // FH
    for(const auto& id : fhGeometry()->getValidGeomDetIds())
    {
        if(id.rawId()==0) continue;
        HGCalDetId waferDetId(id); 
        short module      = wafer_to_module_fh[waferDetId.wafer()];
        short triggercell_offset = module_to_wafers_fh.at(module).at(waferDetId.wafer());
        int nCells = fhTopology().dddConstants().numberCellsHexagon(waferDetId.wafer());
        for(int c=0;c<nCells;c++)
        {
            short triggerCellId = cells_to_trigger_cells[std::make_pair(waferDetId.waferType(),c)];
            // Fill cell -> trigger cell mapping
            HGCalDetId cellDetId(ForwardSubdetector(waferDetId.subdetId()), waferDetId.zside(), waferDetId.layer(), waferDetId.waferType(), waferDetId.wafer(), c); 
            // The module id is used instead of the wafer id for the trigger cells
            // Since there are several wafers per module, an offset is applied on the HGCalDetId::cell field
            if(triggercell_offset+triggerCellId >= HGCalDetId::kHGCalCellMask) edm::LogError("HGCalTriggerGeometry") << "Trigger cell id requested with a cell field larger than available in HGCalDetId (" << triggercell_offset+triggerCellId << " >= " << HGCalDetId::kHGCalCellMask << ")\n";
            HGCalDetId triggerCellDetId(ForwardSubdetector(waferDetId.subdetId()), waferDetId.zside(), waferDetId.layer(), waferDetId.waferType(), module, triggercell_offset + triggerCellId);
            cells_to_trigger_cells_.emplace(cellDetId, triggerCellDetId);
            // Fill trigger cell -> module mapping
            HGCalDetId moduleDetId(ForwardSubdetector(waferDetId.subdetId()), waferDetId.zside(), waferDetId.layer(), waferDetId.waferType(), module, HGCalDetId::kHGCalCellMask);
            trigger_cells_to_modules_.emplace(triggerCellDetId, moduleDetId); 
        }
    }
}



/*****************************************************************/
void HGCalTriggerGeometryHexImp1::buildTriggerCellsAndModules()
/*****************************************************************/
{
    //
    // Build trigger cells and fill map
    typedef HGCalTriggerGeometry::TriggerCell::list_type list_cells;
    // make list of cells in trigger cells
    std::unordered_map<unsigned, list_cells> trigger_cells_to_cells;
    for(const auto& cell_triggerCell : cells_to_trigger_cells_)
    {
        unsigned cell        = cell_triggerCell.first;
        unsigned triggerCell = cell_triggerCell.second;
        auto itr_exist       = trigger_cells_to_cells.emplace(triggerCell, list_cells());
        itr_exist.first->second.emplace(cell);
        //trigger_cells_to_cells.at(triggerCell).emplace(cell);
    }
    for(const auto& triggerCell_cells : trigger_cells_to_cells)
    {
        unsigned triggerCellId = triggerCell_cells.first;
        list_cells cellIds = triggerCell_cells.second;
        // Position: barycenter of the trigger cell.
        Basic3DVector<float> triggerCellVector(0.,0.,0.);
        for(const auto& cell : cellIds)
        {
            HGCalDetId cellDetId(cell);
            triggerCellVector += (cellDetId.subdetId()==ForwardSubdetector::HGCEE ? eeGeometry()->getPosition(cellDetId) :  fhGeometry()->getPosition(cellDetId)).basicVector();
        }
        GlobalPoint triggerCellPoint( triggerCellVector/cellIds.size() );
        const auto& triggerCellToModuleItr = trigger_cells_to_modules_.find(triggerCellId);
        unsigned moduleId = (triggerCellToModuleItr!=trigger_cells_to_modules_.end() ? triggerCellToModuleItr->second : 0); // 0 if the trigger cell doesn't belong to a module
        //unsigned moduleId = trigger_cells_to_modules_.at(triggercellId);
        // FIXME: empty neighbours
        trigger_cells_.emplace(triggerCellId, std::make_unique<const HGCalTriggerGeometry::TriggerCell>(triggerCellId, moduleId, triggerCellPoint, list_cells(), cellIds));
    }
    //
    // Build modules and fill map
    typedef HGCalTriggerGeometry::Module::list_type list_triggerCells;
    typedef HGCalTriggerGeometry::Module::tc_map_type tc_map_to_cells;
    // make list of trigger cells in modules
    std::unordered_map<unsigned, list_triggerCells> modules_to_trigger_cells;
    for(const auto& triggerCell_module : trigger_cells_to_modules_)
    {
        unsigned triggerCell = triggerCell_module.first;
        unsigned module      = triggerCell_module.second;
        auto itr_exist       = modules_to_trigger_cells.emplace(module, list_triggerCells());
        itr_exist.first->second.emplace(triggerCell);
        //modules_to_trigger_cells.at(module).emplace(triggerCell);
    }
    for(const auto& module_triggerCell : modules_to_trigger_cells)
    {
        unsigned moduleId = module_triggerCell.first;
        list_triggerCells triggerCellIds = module_triggerCell.second;
        tc_map_to_cells cellsInTriggerCells;
        // Position: barycenter of the module, from trigger cell positions
        Basic3DVector<float> moduleVector(0.,0.,0.);
        for(const auto& triggerCell : triggerCellIds)
        {
            const auto& cells_in_tc = trigger_cells_to_cells.at(triggerCell);
            for( const unsigned cell : cells_in_tc ) 
            {
              cellsInTriggerCells.emplace(triggerCell,cell);
            }
            moduleVector += trigger_cells_.at(triggerCell)->position().basicVector();
        }
        GlobalPoint modulePoint( moduleVector/triggerCellIds.size() );
        // FIXME: empty neighbours
        modules_.emplace(moduleId, std::make_unique<const HGCalTriggerGeometry::Module>(moduleId, modulePoint, list_triggerCells(), triggerCellIds, cellsInTriggerCells));
    }
}


DEFINE_EDM_PLUGIN(HGCalTriggerGeometryFactory, 
        HGCalTriggerGeometryHexImp1,
        "HGCalTriggerGeometryHexImp1");