CaloTowerConstituentsMap.cc

Go to the documentation of this file.

#include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "Geometry/CaloTopology/interface/HcalTopology.h"
#include "Geometry/CaloTopology/interface/CaloTowerTopology.h"

#include <memory>

//#define EDM_ML_DEBUG

CaloTowerConstituentsMap::~CaloTowerConstituentsMap() {
  delete m_reverseItems.load();
  m_reverseItems = nullptr;
}
CaloTowerConstituentsMap::CaloTowerConstituentsMap(const HcalTopology * hcaltopo, const CaloTowerTopology * cttopo) :
  m_hcaltopo(hcaltopo),
  m_cttopo(cttopo),
  standardHB_(false),
  standardHE_(false),
  standardHF_(false),
  standardHO_(false),
  standardEB_(false),
  m_reverseItems(nullptr)
{
}

CaloTowerDetId CaloTowerConstituentsMap::towerOf(const DetId& id) const {
  CaloTowerDetId tid; // null to start with

  edm::SortedCollection<MapItem>::const_iterator i=m_items.find(id);
  if (i!=m_items.end()) tid=i->tower;

  //use hcaltopo when dealing with hcal detids
  if (tid.null()) {
    if (id.det()==DetId::Hcal) { 
      HcalDetId hid(id);
      if ( (hid.subdet()==HcalBarrel && standardHB_ )  ||
       (hid.subdet()==HcalEndcap && standardHE_ )  ||
       (hid.subdet()==HcalOuter  && standardHO_ )  ||
       (hid.subdet()==HcalForward && standardHF_) ) {
        tid = CaloTowerDetId(m_cttopo->convertHcaltoCT(hid.ietaAbs(),hid.subdet())*hid.zside(),hid.iphi());
      }      
    } else if (id.det()==DetId::Ecal) {
      EcalSubdetector esd=(EcalSubdetector)id.subdetId();
      if (esd==EcalBarrel && standardEB_) {
    EBDetId ebid(id);
    tid=CaloTowerDetId(ebid.tower_ieta(),ebid.tower_iphi());
      }
    }
  }

  return tid;
}

void CaloTowerConstituentsMap::assign(const DetId& cell, const CaloTowerDetId& tower) {
  if (m_items.find(cell)!=m_items.end()) {
    throw cms::Exception("CaloTowers") << "Cell with id " << std::hex << cell.rawId() << std::dec << " is already mapped to a CaloTower " << m_items.find(cell)->tower << std::endl;
  }
  m_items.emplace_back(MapItem(cell,tower));
}

void CaloTowerConstituentsMap::sort() {
  m_items.sort();
  
//  for (auto const & it : m_items)
//    std::cout << std::hex << it.cell.rawId() << " " << it.tower.rawId() << std::dec << std::endl;

}

std::vector<DetId> CaloTowerConstituentsMap::constituentsOf(const CaloTowerDetId& id) const {
#ifdef EDM_ML_DEBUG
  std::cout << "Get constituent of " << std::hex << id.rawId() << std::dec
        << " ID " << id << " ieta " << id.ieta() << std::endl;
#endif
  std::vector<DetId> items;
  if (id.ieta() == 0) return items;

  // build reverse map if needed
  if(!m_reverseItems.load(std::memory_order_acquire)) {
      std::unique_ptr<std::multimap<CaloTowerDetId,DetId>> ptr{new std::multimap<CaloTowerDetId,DetId>};
      for (auto m_item : m_items)
          ptr->insert(std::pair<CaloTowerDetId,DetId>(m_item.tower,m_item.cell));
      std::multimap<CaloTowerDetId,DetId>* expected = nullptr;
      if(m_reverseItems.compare_exchange_strong(expected, ptr.get(), std::memory_order_acq_rel)) {
          ptr.release();
      }
  }

  std::multimap<CaloTowerDetId,DetId>::const_iterator j;
  auto range=(*m_reverseItems.load(std::memory_order_acquire)).equal_range(id);
  for (j=range.first; j!=range.second; j++)
    items.emplace_back(j->second);

  // dealing with topo dependency...
  //use cttopo when dealing with calotower detids
  int nd, sd;
  int hcal_ieta = m_cttopo->convertCTtoHcal(id.ietaAbs());
  
  if (standardHB_) {
    if (id.ietaAbs()<=m_cttopo->lastHBRing()) {
      m_hcaltopo->depthBinInformation(HcalBarrel,hcal_ieta,id.iphi(),id.zside(),nd,sd);
      for (int i=0; i<nd; i++) {
    if (m_hcaltopo->getMergePositionFlag()) {
      HcalDetId hid = m_hcaltopo->mergedDepthDetId(HcalDetId(HcalBarrel,hcal_ieta*id.zside(),id.iphi(),i+sd));
      if (std::find(items.begin(),items.end(),hid) == items.end()) {
        items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
        std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
      }
    } else {
      HcalDetId hid(HcalBarrel,hcal_ieta*id.zside(),id.iphi(),i+sd);
      items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
      std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
    }
      }
    }
  }
  if (standardHO_) {
    if (id.ietaAbs()<=m_cttopo->lastHORing()) {
      m_hcaltopo->depthBinInformation(HcalOuter,hcal_ieta,id.iphi(),id.zside(),nd,sd);
      for (int i=0; i<nd; i++) {
    HcalDetId hid(HcalOuter,hcal_ieta*id.zside(),id.iphi(),i+sd);
        items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
    std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
      }
    }
  }
  if (standardHE_) {
    if (id.ietaAbs()>=m_cttopo->firstHERing() && id.ietaAbs()<=m_cttopo->lastHERing()) {
      m_hcaltopo->depthBinInformation(HcalEndcap,hcal_ieta,id.iphi(),id.zside(),nd,sd);
      for (int i=0; i<nd; i++) {
    if (m_hcaltopo->getMergePositionFlag()) {
      HcalDetId hid = m_hcaltopo->mergedDepthDetId(HcalDetId(HcalEndcap,hcal_ieta*id.zside(),id.iphi(),i+sd));
      if (std::find(items.begin(),items.end(),hid) == items.end()) {
        items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
        std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
      } 
    } else {
      HcalDetId hid(HcalEndcap,hcal_ieta*id.zside(),id.iphi(),i+sd);
      items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
      std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
    }
      }
    }
  }
  if (standardHF_) {
    if (id.ietaAbs()>=m_cttopo->firstHFRing() && id.ietaAbs()<=m_cttopo->lastHFRing()) { 
      m_hcaltopo->depthBinInformation(HcalForward,hcal_ieta,id.iphi(),id.zside(),nd,sd);
      for (int i=0; i<nd; i++) {
    HcalDetId hid(HcalForward,hcal_ieta*id.zside(),id.iphi(),i+sd);
        items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
    std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
      }
      // special handling for first HF tower
      if (id.ietaAbs() == m_cttopo->firstHFRing()) {
        int hcal_ieta2 = hcal_ieta-1;
        m_hcaltopo->depthBinInformation(HcalForward,hcal_ieta2,id.iphi(),id.zside(),nd,sd);
        for (int i=0; i<nd; i++) {
      HcalDetId hid(HcalForward,hcal_ieta2*id.zside(),id.iphi(),i+sd);
          items.emplace_back(hid);
#ifdef EDM_ML_DEBUG
      std::cout << id << " Depth " << i << ":" << i+sd << " " << hid <<"\n";
#endif
    }
      }
    }
  }
  if (standardEB_ && hcal_ieta<=EBDetId::MAX_IETA/5) {
    HcalDetId hid(HcalBarrel,hcal_ieta*id.zside(),id.iphi(),1); // for the limits
    int etaMin, etaMax;
    if (hid.zside() == -1) {
      etaMin = hid.crystal_ieta_high();
      etaMax = hid.crystal_ieta_low();
    } else {
      etaMin = hid.crystal_ieta_low();
      etaMax = hid.crystal_ieta_high();
    }
    for (int ie=etaMin; ie<=etaMax; ie++)
      for (int ip=hid.crystal_iphi_low(); ip<=hid.crystal_iphi_high(); ip++)
        items.emplace_back(EBDetId(ie,ip));
  }
  return items;
}

void CaloTowerConstituentsMap::useStandardHB(bool use) {
  standardHB_=use;
}
void CaloTowerConstituentsMap::useStandardHE(bool use) {
  standardHE_=use;
}
void CaloTowerConstituentsMap::useStandardHO(bool use) {
  standardHO_=use;
}
void CaloTowerConstituentsMap::useStandardHF(bool use) {
  standardHF_=use;
}
void CaloTowerConstituentsMap::useStandardEB(bool use) {
  standardEB_=use;
}