SCTLorentzMonAlg.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
*/
 
#include "SCTLorentzMonAlg.h"
 
#include "AthenaMonitoringKernel/Monitored.h"
#include "InDetIdentifier/SCT_ID.h"
#include "InDetReadoutGeometry/SiDetectorElement.h"
#include "InDetPrepRawData/SiCluster.h"
#include "InDetRIO_OnTrack/SiClusterOnTrack.h"
#include "TrkTrackSummary/TrackSummary.h"
#include "TrkEventUtils/PRDtoTrackMap.h"
 
#include <cmath>
#include <memory>
#include <vector>
 
SCTLorentzMonAlg::SCTLorentzMonAlg(const std::string& name, ISvcLocator* pSvcLocator)
  :AthMonitorAlgorithm(name,pSvcLocator) {
}
 
StatusCode SCTLorentzMonAlg::initialize() {
  ATH_CHECK(detStore()->retrieve(m_pSCTHelper, "SCT_ID"));
  ATH_CHECK(m_trackSummaryTool.retrieve());
  ATH_CHECK(m_tracksName.initialize());
  ATH_CHECK(m_SCTDetEleCollKey.initialize());
  ATH_CHECK(m_assoTool.retrieve(DisableTool{!m_rejectSharedHits} ));
 
  ATH_CHECK(m_holeSearchTool.retrieve());
  ATH_MSG_INFO("Retrieved hole search tool " << m_holeSearchTool);
 
  if ( m_beamSpotKey.initialize().isFailure() ) {
    ATH_MSG_WARNING("Failed to retrieve beamspot service " << m_beamSpotKey << " - will use nominal beamspot at (0,0,0)");
    m_hasBeamCondSvc = false;
  } 
  else {
    m_hasBeamCondSvc = true;
    ATH_MSG_DEBUG("Retrieved service " << m_beamSpotKey);
  }
 
 
  return AthMonitorAlgorithm::initialize();
}
 
StatusCode SCTLorentzMonAlg::fillHistograms(const EventContext& ctx) const {
  static const int layer100[] = {
    2, 2, 3, 2, 2, 2, 0, 2, 3, 2, 0, 2, 3, 2, 3, 2, 0, 2, 3, 0, 2, 0, 2, 3, 2, 2, 2, 0, 0, 0, 0, 0, 0, 3, 0, 3, 2, 0, 2,
    2, 0, 3, 3, 3, 0, 2, 2, 2, 2, 2, 2, 2, 3, 2, 2, 3, 3, 2, 2, 2, 2, 2, 3, 3, 2, 3, 2, 2, 2, 3, 3, 3, 2, 2, 2, 2, 3, 3,
    2, 3, 2, 3, 3, 2, 3, 2, 2, 2, 2, 2, 2, 2
  };
  static const int phi100[] = {
    29, 29, 6, 13, 23, 13, 14, 29, 9, 29, 14, 29, 9, 29, 39, 32, 21, 32, 13, 22, 32, 22, 32, 13, 32, 32, 32, 20, 20, 20,
    20, 20, 20, 13, 21, 17, 33, 5, 33, 33, 31, 6, 19, 47, 21, 37, 37, 37, 37, 33, 37, 37, 24, 33, 33, 47, 19, 33, 33,
    37, 37, 37, 55, 9, 38, 24, 37, 38, 8, 9, 9, 26, 38, 38, 38, 38, 39, 39, 38, 11, 45, 54, 54, 24, 31, 14, 47, 45, 47,
    47, 47, 47
  };
  static const int eta100[] = {
    3, -4, -6, 2, 6, 3, -5, -1, 6, -2, -6, -5, 5, -3, 2, 6, -3, 5, 5, 3, 4, 2, 2, 2, -1, -3, -4, 1, -1, -2, -3, -4, 4,
    -1, -5, 6, 2, 4, 3, 1, 6, -2, 6, 3, -6, -1, 2, 1, 3, -5, 4, 5, -3, -4, -3, -5, -2, -1, -2, -3, -2, -4, -3, 2, 3, -6,
    -5, 4, 6, 1, -6, 1, 1, -5, -4, -3, -3, -5, -2, 1, 5, 5, 4, 4, 5, 4, -1, -5, 3, 4, 1, -5
  };
  constexpr unsigned int layer100_n{sizeof(layer100) / sizeof(*layer100)};
  constexpr unsigned int phi100_n{sizeof(phi100) / sizeof(*phi100)};
  constexpr unsigned int eta100_n{sizeof(eta100) / sizeof(*eta100)};
  constexpr bool theseArraysAreEqualInLength{(layer100_n == phi100_n) and (phi100_n == eta100_n)};
  static_assert(theseArraysAreEqualInLength, "Coordinate arrays for <100> wafers are not of equal length");
 
  ATH_MSG_DEBUG("enters fillHistograms");
 
  SG::ReadCondHandle<InDetDD::SiDetectorElementCollection> sctDetEle{m_SCTDetEleCollKey, ctx};
  const InDetDD::SiDetectorElementCollection* elements{sctDetEle.retrieve()};
  if (elements==nullptr) {
    ATH_MSG_WARNING(m_SCTDetEleCollKey.fullKey() << " could not be retrieved");
    return StatusCode::SUCCESS;
  }
  SG::ReadHandle<TrackCollection> tracks{m_tracksName, ctx};
  if (not tracks.isValid()) {
    ATH_MSG_WARNING(" TrackCollection not found: Exit SCTLorentzTool" << m_tracksName.key());
    return StatusCode::SUCCESS;
  }
 
  // Prepare AssociationTool
  std::unique_ptr<Trk::PRDtoTrackMap> prd_to_track_map;
  if (m_rejectSharedHits) {
     prd_to_track_map = m_assoTool->createPRDtoTrackMap();
     for (const Trk::Track* track : *tracks) {
        ATH_CHECK( m_assoTool->addPRDs(*prd_to_track_map,*track) );
     }
  }
 
  float beamSpotX = 0.;
  float beamSpotY = 0.;
  float beamSpotZ = 0.;
  float beamTiltX = 0.;
  float beamTiltY = 0.;
 
  if (m_hasBeamCondSvc) {
    auto beamSpotHandle = SG::ReadCondHandle(m_beamSpotKey, ctx);
    Amg::Vector3D bpos = beamSpotHandle->beamPos();
    beamSpotX = bpos.x();
    beamSpotY = bpos.y();
    beamSpotZ = bpos.z();
    beamTiltX = beamSpotHandle->beamTilt(0);
    beamTiltY = beamSpotHandle->beamTilt(1);
    ATH_MSG_DEBUG ("Beamspot: x0 = " << beamSpotX << ", y0 = " << beamSpotY << ", z0 = " << beamSpotZ << ", tiltX = " << beamTiltX << ", tiltY = " << beamTiltY);
  }
 
  for (const Trk::Track* track: *tracks) {
    if (track==nullptr) {
      ATH_MSG_ERROR("no pointer to track!!!");
      continue;
    }
 
    if (track->perigeeParameters() == nullptr) {
      continue;
    }
    // Get a track with holes
    std::unique_ptr<const Trk::Track> trackWithHoles(m_holeSearchTool->getTrackWithHoles(*track));
 
    if (not trackWithHoles) {
      ATH_MSG_WARNING("trackWithHoles pointer is invalid");
      continue;
    }
 
    ATH_MSG_VERBOSE("Found " << trackWithHoles->trackStateOnSurfaces()->size() << " states on track");
 
    // Get pointer to track state on surfaces
    const Trk::TrackStates* trackStates{trackWithHoles->trackStateOnSurfaces()};
    if (trackStates==nullptr) {
      ATH_MSG_WARNING("for current track, TrackStateOnSurfaces == Null, no data will be written for this track");
      continue;
    }
 
    const Trk::TrackSummary* summary{trackWithHoles->trackSummary()};
    std::unique_ptr<Trk::TrackSummary> mySummary;
    if (summary==nullptr) {
      mySummary = m_trackSummaryTool->summary(ctx, *trackWithHoles);
      summary = mySummary.get();
      if (summary==nullptr) {
        ATH_MSG_WARNING("Trk::TrackSummary is null and cannot be created by " << m_trackSummaryTool.name());
        continue;
      }
    }
 
    for (const Trk::TrackStateOnSurface* tsos: *trackStates) {
      int nStrip = 0;
      Identifier sct_id(0);
 
      if (tsos->type(Trk::TrackStateOnSurface::Measurement)) {
    // This is hit.
    const InDet::SiClusterOnTrack* clus{dynamic_cast<const InDet::SiClusterOnTrack*>(tsos->measurementOnTrack())};
        if (clus) { // Is it a SiCluster? If yes...
          // Reject shared hits if you want
          if (prd_to_track_map and prd_to_track_map->isShared(*(clus->prepRawData())) ) {
            continue;
          }
 
          const InDet::SiCluster* RawDataClus{dynamic_cast<const InDet::SiCluster*>(clus->prepRawData())};
          if (RawDataClus==nullptr) {
            continue; // Continue if dynamic_cast returns null
          }
          if (RawDataClus->detectorElement()->isSCT()) {
        sct_id = clus->identify();
        // find cluster size
        const std::vector<Identifier>& rdoList{RawDataClus->rdoList()};
        nStrip = static_cast<int>(rdoList.size());
      } else {
        continue;
      }
    } else {
      // This is not an Si cluster.
      continue;
    }
      } else if (tsos->type(Trk::TrackStateOnSurface::Hole)) {
    // This is a hole.
    const Trk::MeasurementBase* mesb{tsos->measurementOnTrack()};
    if (mesb and mesb->associatedSurface().associatedDetectorElement()) {
      sct_id = mesb->associatedSurface().associatedDetectorElement()->identify();
    } else if (tsos->trackParameters()) {
      sct_id = tsos->trackParameters()->associatedSurface().associatedDetectorElementIdentifier();
    } else {
      // sct_id wasn't correctly retrieved.
      continue;
    }
      } else {
    continue;
      }
 
      if (!m_pSCTHelper->is_sct(sct_id)) {
    continue;
      }
      
      const int bec{m_pSCTHelper->barrel_ec(sct_id)};
      const int layer{m_pSCTHelper->layer_disk(sct_id)};
      const int side{m_pSCTHelper->side(sct_id)};
      const int eta{m_pSCTHelper->eta_module(sct_id)};
      const int phi{m_pSCTHelper->phi_module(sct_id)};
 
      // Check if the silicon surface is 100.
      SiliconSurface surface{surface111};
      for (unsigned int i{0}; i < layer100_n; i++) {
    if ((layer100[i] == layer) and (eta100[i] == eta) and (phi100[i] == phi)) {
      surface = surface100;
      break;
    }
      }
 
      const Trk::TrackParameters* trkp{dynamic_cast<const Trk::TrackParameters*>(tsos->trackParameters())};
      if (trkp==nullptr) {
    ATH_MSG_WARNING(" Null pointer to MeasuredTrackParameters");
    continue;
      }
 
      const Trk::Perigee* perigee{track->perigeeParameters()};
      
      if (perigee) {
    // Get angle to wafer surface
    float phiToWafer{90.f};
    float thetaToWafer(90.f);
    float sinAlpha{0.f}; // for barrel, which is the only thing considered here
    float pTrack[3]; // 3 is for x, y, z.
    pTrack[0] = trkp->momentum().x();
    pTrack[1] = trkp->momentum().y();
    pTrack[2] = trkp->momentum().z();
    float etaTrack = trkp->eta();
 
    int iflag{findAnglesToWaferSurface(pTrack, sinAlpha, sct_id, elements, thetaToWafer, phiToWafer)};
    if (iflag < 0) {
      ATH_MSG_WARNING("Error in finding track angles to wafer surface");
      continue; // Let's think about this (later)... continue, break or return?
    }
 
        float beamX          = 0;
        float beamY          = 0;
        float d0bscorr       = -999;
 
        // correct the track parameters for the beamspot position
        beamX = beamSpotX + tan(beamTiltX) * (perigee->parameters()[Trk::z0]-beamSpotZ);
        beamY = beamSpotY + tan(beamTiltY) * (perigee->parameters()[Trk::z0]-beamSpotZ);
        d0bscorr = perigee->parameters()[Trk::d0] - ( -sin(perigee->parameters()[Trk::phi])*beamX + cos(perigee->parameters()[Trk::phi])*beamY );
 
 
    // The cuts for the physics runs follow ATL-COM-INDET-2021-011.
    if (trkp->momentum().perp() > 500. and
        summary->get(Trk::numberOfSCTHits) > 6 and
            (dataType() == AthMonitorAlgorithm::DataType_t::cosmics or
             (std::abs(d0bscorr) < 1. and trackWithHoles->perigeeParameters()->parameters()[Trk::qOverP] < 0.)
            )
        ) {
    // Fill profile
      std::string xVar{"phiToWafer"};
      std::string yVar{"nStrip"};
      if(bec == 0){ // Barrel
        xVar += "_" + std::to_string(layer);
        yVar += "_" + std::to_string(layer);
        if(surface == surface100){
          xVar += "_100";
          yVar += "_100";
        }
        if(surface == surface111){
          xVar += "_111";
          yVar += "_111";
        }
      } else { // Endcaps
        if (bec == -2) {
          xVar += "_ECC";
          yVar += "_ECC";
        } else {
          xVar += "_ECA";
          yVar += "_ECA";
        }
        xVar += std::to_string(layer);
        yVar += std::to_string(layer);
        if (eta == 0) {
          xVar += "_outer";
          yVar += "_outer";
        } else if (eta == 1) {
          xVar += "_middle";
          yVar += "_middle";
        } else {
          xVar += "_inner";
          yVar += "_inner";
        }
      } // Common for barrel and endcaps
      if(side == side0){
        xVar += "_0";
        yVar += "_0";
      }
      if(side == side1){
        xVar += "_1";
        yVar += "_1";
      }
 
      auto phiToWaferAcc{Monitored::Scalar<float>(xVar, phiToWafer)};
      auto nStripAcc{Monitored::Scalar<int>(yVar, nStrip)};
      auto isCentralAcc{Monitored::Scalar<bool>("isCentral", (std::abs(etaTrack) < 0.5))};
      fill("SCTLorentzMonitor", phiToWaferAcc, nStripAcc, isCentralAcc);
    }// end if passesCuts
      } // end if SCT..
    }// end of loop on TrackStatesonSurface (they can be SiClusters, TRTHits,..)
  } // end of loop on tracks
 
  return StatusCode::SUCCESS;
}
 
int
SCTLorentzMonAlg::findAnglesToWaferSurface(const float (&vec)[3], // 3 is for x, y, z.
                                           const float& sinAlpha, const Identifier& id,
                                           const InDetDD::SiDetectorElementCollection* elements,
                                           float& theta, float& phi) const {
  int iflag{-1};
 
  phi = 90.;
  theta = 90.;
 
  const Identifier waferId{m_pSCTHelper->wafer_id(id)};
  const IdentifierHash waferHash{m_pSCTHelper->wafer_hash(waferId)};
  const InDetDD::SiDetectorElement* element{elements->getDetectorElement(waferHash)};
  if (element==nullptr) {
    ATH_MSG_ERROR("findAnglesToWaferSurface:  failed to find detector element for id=" <<
                  m_pSCTHelper->show_to_string(id));
    return iflag;
  }
 
  float cosAlpha{std::sqrt(1.0f - sinAlpha * sinAlpha)};
  double phix{ cosAlpha * element->phiAxis().x() + sinAlpha * element->phiAxis().y()};
  double phiy{-sinAlpha * element->phiAxis().x() + cosAlpha * element->phiAxis().y()};
 
  double pNormal{vec[0] * element->normal().x()  + vec[1] * element->normal().y()  + vec[2] * element->normal().z()};
  double pEta{   vec[0] * element->etaAxis().x() + vec[1] * element->etaAxis().y() + vec[2] * element->etaAxis().z()};
  double pPhi{   vec[0] * phix                   + vec[1] * phiy                   + vec[2] * element->phiAxis().z()};
 
  if (pPhi < 0.) {
    phi = -90.;
  }
  if (pEta < 0.) {
    theta = -90.;
  }
  if (pNormal != 0.) {
    phi = std::atan(pPhi / pNormal) / CLHEP::deg;
    theta = std::atan(pEta / pNormal) / CLHEP::deg;
  }
  iflag = 1;
  return iflag;
}