Muon/L1MuonNSW/src/MMT_Diamond.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#include "MMT_Diamond.h"
#include "AthenaKernel/getMessageSvc.h"
#include "Acts/Utilities/AngleHelpers.hpp"
#include "Acts/Utilities/MathHelpers.hpp"
#include <algorithm>
#include <numbers>
 
namespace {
  constexpr int bc_wind = 4; // fixed time window (in bunch crossings) during which the algorithm collects ART hits
}
 
namespace L1Muon::L1MM {
 
MMT_Diamond::MMT_Diamond(const int diamXthreshold, const bool uv, const int diamUVthreshold, const int roadSize,
    const int olapEtaUp, const int olapEtaDown, const int olapStereoUp, const int olapStereoDown):
    AthMessaging(Athena::getMessageSvc(), "L1MM::MMT_Diamond"),
    m_uvflag{uv},
    m_roadSize{roadSize},
    m_roadSizeUpX{olapEtaUp},
    m_roadSizeDownX{olapEtaDown},
    m_roadSizeUpUV{olapStereoUp},
    m_roadSizeDownUV{olapStereoDown},
    m_xthr{diamXthreshold},
    m_uvthr{diamUVthreshold} {
}
 
void MMT_Diamond::createRoads(std::vector<MMT_Road>& roads, const bool isLarge, const bool isEta1) const {
  const char sec = (isLarge) ? 'L' : 'S';
  /*
   * This computation is done as follows: 1024 X roads, from 8192 strips in total (5120 for |eta|==1 and 3072 for |eta|==2)
   * MML: for i in [0,8] -> i*6 UV roads. Then: (1024-9)*6*9 UV roads
   * MMS: for i in [0,6] -> i*6 UV roads. Then: (1024-7)*6*7 UV roads
   */
  const unsigned int vecRoads = 35000;
  roads.reserve(vecRoads);
  const int nroad = kStripsPerLayer/m_roadSize;
  for (int i = 0; i < nroad; ++i) {
    const int div = kStripsEta1/m_roadSize;
    if((isEta1 && i>div) || (!isEta1 && i<=div)) continue;
    roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i);
 
    /*
     * The computation of "nuv" is:
     * B = (1./std::tan(1.5/180.*M_PI));
     * nuv = std::round( par->getlWidth() / (B * 0.4 * 2.)/this->getRoadSize() );
     * As getlWidth() has to deal with the full wedge, only two (fixed by construction) integer values are allowed, according to small (7) or large (9) wedge
     */
    if(m_uvflag) {
      const int nuv = (isLarge) ? 9 : 7;
      for (int uv = 1; uv <= nuv; uv++) {
        if (i-uv < 0) continue;
        roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i, i+uv, i-uv);
        roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i, i-uv, i+uv);
        roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i, i+uv-1, i-uv);
        roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i, i-uv, i+uv-1);
        roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i, i-uv+1, i+uv);
        roads.emplace_back(sec, m_roadSize, m_roadSizeUpX, m_roadSizeDownX, m_roadSizeUpUV, m_roadSizeDownUV, m_xthr, m_uvthr, i, i+uv, i-uv+1);
      }
    }
  }
}
 
void MMT_Diamond::findDiamonds(std::vector<MMT_Hit>& hits, std::vector<MMT_Road>& roads, std::vector<slope_t>& diamondSlopes,
                               const int sectorPhi, const double tanStereoAngle) const {
 
  // Comparison with lambda function (easier to implement)
  std::ranges::sort(hits, [](const MMT_Hit& h1, const MMT_Hit& h2){ return h1.getBC() < h2.getBC(); });
  const int bc_start = hits.front().getBC();
  const int bc_end = hits.front().getBC() + 16;
  ATH_MSG_DEBUG("Window Start: " << bc_start << " - Window End: " << bc_end);
 
  // Hits of the current bunch crossing (non-owning, they point into hits, which is not modified below)
  std::vector<const MMT_Hit*> hits_now;
  std::vector< std::pair<int, float> > vmm_same;
  std::vector< std::pair<int, int> > addc_same;
  std::vector<int> to_erase;
 
  // each road makes independent triggers, evaluated on each BC
  unsigned int ibc = 0;
  for (int bc = hits.front().getBC(); bc < bc_end; bc++) {
    // Cleaning stuff
    hits_now.clear();
 
    for (unsigned int j = ibc; j < hits.size(); j++) {
      if (hits[j].getBC() == bc) hits_now.push_back(&hits[j]);
      else if (hits[j].getBC() > bc) {
        ibc = j;
        break;
      }
    }
 
    // Simulate harware filters: ART hits + ADDC filter
    for (unsigned int ib = 0; ib < kPlanes; ib++) { //loop on plane from 0 to 7
      // VMM-ART-hit filter
      for (int j = 0; j < kVMMsPerLayer; j++) {
        vmm_same.clear();
        unsigned int k = 0;
        for (const auto &hit_pointer: hits_now) {
          if (static_cast<unsigned int>(hit_pointer->getPlane()) != ib) continue;
          if (hit_pointer->getVMM() == j){
            vmm_same.push_back( std::make_pair(k, hit_pointer->getTime()) );
          }
          k++;
        }
        if (vmm_same.size() > 1) {
          to_erase.clear();
          std::ranges::sort(vmm_same, {}, &std::pair<int, float>::second);
          for (auto pair: vmm_same) to_erase.push_back(pair.first);
          // reverse and erase
          std::sort(to_erase.rbegin(), to_erase.rend());
          to_erase.pop_back(); //removes the hit with the earliest time
          for (auto l : to_erase) {
            hits_now.erase(hits_now.begin() + l);
          }
        }
      }
      // ADDC-like filter
      for (int ia = 0; ia < kARTsPerLayer; ia++) { // From 0 to 3 (local index of the ART ASIC in the layer)
        addc_same.clear();
        for (unsigned int k = 0; k < hits_now.size(); k++) {
          if (static_cast<unsigned int>(hits_now[k]->getPlane()) != ib) continue;
          // the global strip index on the sector layer (getChannel returns the module's local strip index)
          if (hits_now[k]->getART() == ia) addc_same.emplace_back(k, hits_now[k]->getGlobalStrip());
        }
 
        if (addc_same.size() > kMaxHitsPerART) {
          // priority encode the hits by channel number; remember hits 8+
          to_erase.clear();
 
          std::ranges::sort(addc_same, {}, &std::pair<int, int>::second);
          for (unsigned int it = kMaxHitsPerART; it < addc_same.size(); it++) to_erase.push_back(addc_same[it].first);
 
          // reverse and erase
          std::sort(to_erase.rbegin(), to_erase.rend());
          for (auto l : to_erase) {
            hits_now.erase(hits_now.begin() + l);
          }
        }
      }
    } // loop on plane for VMM and ART ASIC filter
 
    for (auto &road : roads) {
 
      if (!road.getHitVector().empty()) road.incrementAge(bc_wind);
      if (!hits_now.empty()) road.addHits(hits_now);
 
      if (road.checkCoincidences(bc_wind) && bc >= (bc_start - 1)) {
 
        ATH_MSG_DEBUG("------------------------------------------------------------------");
        ATH_MSG_DEBUG("Coincidence FOUND @BC: " << bc);
        ATH_MSG_DEBUG("Road (x, u, v, count): (" << road.iRoadx() << ", " << road.iRoadu() << ", " << road.iRoadv() << ", " << road.countHits() << ")");
        ATH_MSG_DEBUG("------------------------------------------------------------------");
 
        std::vector<int> bcidVec;
        for (const auto &hit: road.getHitVector()) {
          bcidVec.push_back(hit.getBC());
        }
        std::ranges::sort(bcidVec);
 
        // evaluating mode of the BCID of the hits in the diamond
        // default setting in the firmware is the mode of the hits's bcid in the diamond
        int bcidVal=bcidVec[0], bcidCount=1, modeCount=1, bcidMode=bcidVec[0];
        for (unsigned int i=1; i<bcidVec.size(); i++){
          if (bcidVec[i] == bcidVal){
            bcidCount++;
          } else {
            bcidCount = 1;
            bcidVal = bcidVec[i];
          }
          if (bcidCount > modeCount) {
            modeCount = bcidCount;
            bcidMode = bcidVal;
          }
        }
 
        slope_t slope;
        slope.BC = bcidMode;
        slope.totalCount = road.countHits();
        slope.iRoad = road.iRoadx();
        slope.iRoadu = road.iRoadu();
        slope.iRoadv = road.iRoadv();
        slope.xCount = road.countXHits();
        slope.uCount = road.countUHits();
        slope.age = slope.BC - bc_start;
        slope.mxl = road.mxl();
        slope.my = road.avgSofXUV(LayerType::X); // defined as my in ATL-COM-UPGRADE-2015-033
        slope.uavg = road.avgSofXUV(LayerType::U);
        slope.vavg = road.avgSofXUV(LayerType::V);
        slope.mx = (slope.uavg-slope.vavg)/(2.*tanStereoAngle);
        const double theta = std::atan(std::hypot(slope.mx, slope.my));
        slope.theta = (slope.my > 0.) ? theta : std::numbers::pi - theta;
        slope.eta = Acts::AngleHelpers::etaFromTheta(slope.theta);
        slope.dtheta = (slope.mxl - slope.my)/(1. + slope.mxl*slope.my);
        slope.side = (slope.my > 0.) ? SideType::A : SideType::C;
        // Local phi relative to the sector centre, in (-pi/2, pi/2) (my is never 0)
        slope.phi = std::atan(slope.mx/slope.my);
        slope.phiShf = phiShift(sectorPhi, slope.phi, slope.side);
        slope.lowRes = road.evaluateLowRes();
 
        diamondSlopes.push_back(slope);
      }
    }
  }
}
 
double MMT_Diamond::phiShift(const int n, const double phi, const SideType side) const {
  const double Phi = -Acts::copySign(phi, side);  // local phi flipped on side A
  double shift = (n > 8) ? (16-n)*std::numbers::pi/8. : n*std::numbers::pi/8.;
  if (n < 8)       return (Phi + shift);
  else if (n == 8) return (Phi + ((Phi > 0.) ? -1. : 1.)*shift);
  else             return (Phi - shift);
}
 
}  // namespace L1Muon::L1MM