diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/CMakeLists.txt b/Trigger/TrigTools/TrigInDetTrackFitter/CMakeLists.txt
index 15628d4f29ba6199818c3e3d8b468ce900cf0462..f8b430ae84ccf36f7434176a4ed0052575489f64 100644
--- a/Trigger/TrigTools/TrigInDetTrackFitter/CMakeLists.txt
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/CMakeLists.txt
@@ -1,4 +1,4 @@
-# Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
+# Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
# Declare the package name:
atlas_subdir( TrigInDetTrackFitter )
@@ -7,4 +7,4 @@ atlas_subdir( TrigInDetTrackFitter )
atlas_add_component( TrigInDetTrackFitter
src/*.cxx
src/components/*.cxx
- LINK_LIBRARIES AthenaBaseComps AtlasDetDescr GaudiKernel InDetIdentifier InDetPrepRawData InDetRIO_OnTrack MagFieldConditions MagFieldElements StoreGateLib TrigInDetEvent TrigInDetToolInterfacesLib TrkDistributedKalmanFilterLib TrkEventPrimitives TrkParameters TrkPrepRawData TrkRIO_OnTrack TrkSurfaces TrkToolInterfaces TrkTrack )
+ LINK_LIBRARIES AthenaBaseComps AtlasDetDescr GaudiKernel InDetIdentifier InDetPrepRawData InDetRIO_OnTrack MagFieldConditions MagFieldElements StoreGateLib InDetReadoutGeometry PixelReadoutGeometryLib SCT_ReadoutGeometry TrigInDetEvent TrigInDetToolInterfacesLib TrkDistributedKalmanFilterLib TrkEventPrimitives TrkParameters TrkPrepRawData TrkRIO_OnTrack TrkSurfaces TrkToolInterfaces TrkTrack )
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.cxx b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.cxx
index cd98a01cd074f714e21a9d6764fe79b6d537675b..bf826083b75aa81cac929eef28c42e7812c4b452 100644
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.cxx
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.cxx
@@ -10,6 +10,10 @@
#include "TrkParameters/TrackParameters.h"
#include "TrkPrepRawData/PrepRawData.h"
#include "TrkSpacePoint/SpacePoint.h"
+#include "PixelReadoutGeometry/PixelDetectorManager.h"
+#include "SCT_ReadoutGeometry/SCT_DetectorManager.h"
+#include "InDetIdentifier/PixelID.h"
+
#include "TrigInDetRoadPredictorTool.h"
@@ -19,6 +23,8 @@
#include "TrkSurfaces/Surface.h"
#include "InDetPrepRawData/PixelCluster.h"
+
+
TrigInDetRoadPredictorTool::TrigInDetRoadPredictorTool(const std::string& t,
const std::string& n,
const IInterface* p ): AthAlgTool(t,n,p)
@@ -28,15 +34,412 @@ TrigInDetRoadPredictorTool::TrigInDetRoadPredictorTool(const std::string& t,
StatusCode TrigInDetRoadPredictorTool::initialize() {
+ StatusCode sc = m_layerNumberTool.retrieve();
+ if(sc.isFailure()) {
+ ATH_MSG_ERROR("Could not retrieve "<<m_layerNumberTool);
+ return sc;
+ } else {
+ ATH_MSG_INFO("Detector layer structure has "<<m_layerNumberTool->maxNumberOfUniqueLayers()<<" unique layers");
+ }
+
+
ATH_CHECK( m_fieldCondObjInputKey.initialize());
+ ATH_CHECK( detStore()->retrieve(m_pixelManager, "ITkPixel") );
+ ATH_CHECK( detStore()->retrieve(m_stripManager, "ITkStrip") );
+ ATH_CHECK( detStore()->retrieve(m_pixelId, "PixelID") );
+ ATH_CHECK( detStore()->retrieve(m_stripId, "SCT_ID") );
+ buildDetectorDescription();
+
+
return StatusCode::SUCCESS;
}
+void TrigInDetRoadPredictorTool::addNewElement(unsigned int layerID, short phi_idx, short eta_idx, const InDetDD::SiDetectorElement* p) {
+
+ unsigned int hash = p->identifyHash();
+
+ DetectorElementDescription ded(hash);
+
+ //find corners in the global c.s.
+
+ float de_len = 0.5*p->design().length();
+ float de_wmax = 0.5*p->design().maxWidth();
+ float de_wmin = 0.5*p->design().minWidth();
+
+ float dPhi[4] = {de_wmax, de_wmin, -de_wmin, -de_wmax};//locX
+ float dEta[4] = {de_len, -de_len, -de_len, de_len}; //locY
+
+ const Amg::Vector3D& C = p->center();
+ const Amg::Vector3D& PhiAx = p->phiAxis();
+ const Amg::Vector3D& EtaAx = p->etaAxis();
+
+ for(int ic=0; ic<4; ic++) {
+ float x = C.x() + PhiAx.x()*dPhi[ic] + EtaAx.x()*dEta[ic];
+ float y = C.y() + PhiAx.y()*dPhi[ic] + EtaAx.y()*dEta[ic];
+ float z = C.z() + PhiAx.z()*dPhi[ic] + EtaAx.z()*dEta[ic];
+ ded.m_c[ic][0] = std::sqrt(x*x+y*y);//r
+ ded.m_c[ic][1] = z;
+ ded.m_c[ic][2] = std::atan2(y,x);//phi
+ }
+
+ auto& L = (*m_layerMap.find(layerID)).second;
+
+ short prim_idx = L.m_mappingType != 2 ? phi_idx : eta_idx;
+ short sec_idx = L.m_mappingType != 2 ? eta_idx : phi_idx;
+
+ if(L.m_mappingType == 0) {//barrel: primary index is phi, eta/z is secondary
+ ded.m_ref = C.z();
+ ded.m_index = sec_idx;
+ ded.m_minBound = ded.m_c[0][1];
+ ded.m_maxBound = ded.m_c[0][1];
+ for(int ic=1;ic<4;ic++) {
+ if(ded.m_minBound > ded.m_c[ic][1]) ded.m_minBound = ded.m_c[ic][1];
+ if(ded.m_maxBound < ded.m_c[ic][1]) ded.m_maxBound = ded.m_c[ic][1];
+ }
+ }
+
+ if(L.m_mappingType == 1) {//pixel endcap layers : primary index is phi, eta/R is secondary
+ ded.m_ref = C.perp();
+ ded.m_index = sec_idx;
+ ded.m_minBound = ded.m_c[0][0];
+ ded.m_maxBound = ded.m_c[0][0];
+ for(int ic=1;ic<4;ic++) {
+ if(ded.m_minBound > ded.m_c[ic][0]) ded.m_minBound = ded.m_c[ic][0];
+ if(ded.m_maxBound < ded.m_c[ic][0]) ded.m_maxBound = ded.m_c[ic][0];
+ }
+ }
+
+ if(L.m_mappingType == 2) {//strip endcaps: primary index is R, phi is secondary, layers are double
+ ded.m_ref = std::atan2(C.y(),C.x());
+ ded.m_index = sec_idx;
+ ded.m_minBound = ded.m_c[0][2];
+ ded.m_maxBound = ded.m_c[0][2];
+ for(int ic=1;ic<4;ic++) {
+ if(ded.m_minBound > ded.m_c[ic][2]) ded.m_minBound = ded.m_c[ic][2];
+ if(ded.m_maxBound < ded.m_c[ic][2]) ded.m_maxBound = ded.m_c[ic][2];
+ }
+ }
+
+ //put DetEl description into the corresponding collection
+
+ int detColIdx = 0;
+
+ if(L.m_nSubLayers == 2 && (hash % 2) != 0) {//double layer, odd detector elements
+ detColIdx = 1;
+ }
+
+ if(L.m_colls[ detColIdx].find(prim_idx) == L.m_colls[ detColIdx].end()) {
+
+ float primBounds[2] = {0,0};
+ if(L.m_mappingType == 2) {
+ primBounds[0] = p->rMin();
+ primBounds[1] = p->rMax();
+ }
+ else {
+ primBounds[0] = p->phiMin();
+ primBounds[1] = p->phiMax();
+ }
+
+ DetectorElementsCollection dc(prim_idx, primBounds[0], primBounds[1]);
+ L.m_colls[ detColIdx].insert(std::make_pair(prim_idx, dc));
+ }
+
+ (*L.m_colls[detColIdx].find(prim_idx)).second.m_vDE.push_back(ded);
+
+}
+
+
+void TrigInDetRoadPredictorTool::buildDetectorDescription() {
+
+ const std::vector<short>& vPixelL = *(m_layerNumberTool->pixelLayers());
+ const std::vector<TrigInDetSiLayer>& SiL = *(m_layerNumberTool->layerGeometry());
+
+ for(int hash = 0; hash<static_cast<int>(m_pixelId->wafer_hash_max()); hash++) {
+
+ Identifier offlineId = m_pixelId->wafer_id(hash);
+ short barrel_ec = m_pixelId->barrel_ec(offlineId);
+ if(std::abs(barrel_ec)>2) continue;//no DBM needed
+
+ unsigned int layerID = SiL.at(vPixelL.at(hash)).m_subdet;
+ unsigned int volumeID = layerID / 1000;
+
+ if(m_layerMap.find(layerID) == m_layerMap.end()) {
+
+ int nSubLayers = 1;
+ unsigned int mappingType = 1;
+
+ if(barrel_ec == 0) mappingType = 0;
+ else {
+ if(volumeID == 12 || volumeID == 14) mappingType = 2;
+ }
+ m_layerMap.try_emplace(layerID, layerID, nSubLayers, mappingType);
+ }
+
+ short phi_index = m_pixelId->phi_module(offlineId);
+ short eta_index = m_pixelId->eta_module(offlineId);
+
+ addNewElement(layerID, phi_index, eta_index, m_pixelManager->getDetectorElement(hash));
+ }
+
+ const std::vector<short>& vStripL = *(m_layerNumberTool->sctLayers());
+
+ for(int hash = 0; hash<static_cast<int>(m_stripId->wafer_hash_max()); hash++) {
+
+ Identifier offlineId = m_stripId->wafer_id(hash);
+ short barrel_ec = m_stripId->barrel_ec(offlineId);
+ if(std::abs(barrel_ec)>2) continue;
+
+ unsigned int layerID = SiL.at(vStripL.at(hash)).m_subdet;
+ unsigned int volumeID = layerID / 1000;
+
+ if(m_layerMap.find(layerID) == m_layerMap.end()) {
+ int nSubLayers = 2;//strip layers are double
+ unsigned int mappingType = 1;
+
+ if(barrel_ec == 0) mappingType = 0;
+ else {
+ if(volumeID == 12 || volumeID == 14) mappingType = 2;
+ }
+ m_layerMap.try_emplace(layerID, layerID, nSubLayers, mappingType);
+ }
+
+ short phi_index = m_stripId->phi_module(offlineId);
+ short eta_index = m_stripId->eta_module(offlineId);
+
+ addNewElement(layerID, phi_index, eta_index, m_stripManager->getDetectorElement(hash));
+ }
+
+ //building hit boxes
+
+ createHitBoxes();
+
+}
+
+void TrigInDetRoadPredictorTool::createHitBoxes() {
+
+ const float margin_r = 3.0;
+ const float margin_z = 3.0;
+
+ for(const auto& l : m_layerMap) {
+
+ unsigned int layerID = l.first;
+ unsigned int volumeID = layerID / 1000;
+
+ const auto& L = l.second;
+
+ float minR = 1e8;// vert0
+ float minZ = 1e8;// vert1
+ float maxR = -1e8;// vert2
+ float maxZ = -1e8;// vert3
+
+ float vert[4][2]; //z,r
+
+ memset(&vert[0][0],0,sizeof(vert));
+
+ for(int iSL = 0; iSL < L.m_nSubLayers;iSL++) {
+ for(const auto& deColl : L.m_colls[iSL]) {
+ for(const auto& de : deColl.second.m_vDE) {
+ for(int ic=0;ic<4;ic++) {
+ float r = de.m_c[ic][0];
+ float z = de.m_c[ic][1];
+ if(r < minR) {
+ minR = r;
+ vert[0][0] = z;
+ vert[0][1] = r;
+ }
+ if(z < minZ) {
+ minZ = z;
+ vert[1][0] = z;
+ vert[1][1] = r;
+ }
+ if(r > maxR) {
+ maxR = r;
+ vert[2][0] = z;
+ vert[2][1] = r;
+ }
+ if(z > maxZ) {
+ maxZ = z;
+ vert[3][0] = z;
+ vert[3][1] = r;
+ }
+ }
+ }
+ }
+ }
+ minR -= margin_r;
+ vert[0][1] -= margin_r;
+ minZ -= margin_z;
+ vert[1][0] -= margin_z;
+ maxR += margin_r;
+ vert[2][1] += margin_r;
+ maxZ += margin_z;
+ vert[3][0] += margin_z;
+
+ int new_layer_index = -1;
+ LayerBoundary lb;
+ new_layer_index = m_lBoundaries.size();
+ lb.m_index = new_layer_index;
+ lb.m_lay_id = layerID;
+ lb.m_nVertices = 5;
+
+ if(volumeID == 73 || volumeID == 75 || volumeID == 77) {//negative inclined
+ lb.m_z = {vert[3][0], vert[2][0], vert[1][0], vert[0][0], vert[3][0]};
+ lb.m_r = {vert[3][1], vert[2][1], vert[1][1], vert[0][1], vert[3][1]};
+ }
+ else if(volumeID == 93 || volumeID == 95 || volumeID == 97) {//positive inclined
+ lb.m_z = {vert[2][0], vert[1][0], vert[0][0], vert[3][0], vert[2][0]};
+ lb.m_r = {vert[2][1], vert[1][1], vert[0][1], vert[3][1], vert[2][1]};
+ }
+ else {//all other layers
+ lb.m_z = {maxZ, minZ, minZ, maxZ, maxZ};
+ lb.m_r = {maxR, maxR, minR, minR, maxR};
+ }
+
+ m_lBoundaries.push_back(lb);
+
+ bool volExists = false;
+
+ for(auto& v : m_vBoundaries) {
+ if(v.m_vol_id == (int)volumeID) {
+ volExists = true;//update corners
+ if(v.m_zr[0] > minZ) v.m_zr[0] = minZ;
+ if(v.m_zr[1] < maxZ) v.m_zr[1] = maxZ;
+ if(v.m_zr[2] > minR) v.m_zr[2] = minR;
+ if(v.m_zr[3] < maxR) v.m_zr[3] = maxR;
+ v.m_layers.push_back(new_layer_index);
+ break;
+ }
+ }
+ if(!volExists) {//add a new volume with 4 corners
+ VolumeBoundary vb;
+ vb.m_layers.push_back(new_layer_index);
+ vb.m_index = m_vBoundaries.size();
+ vb.m_vol_id = volumeID;
+ vb.m_zr[0] = minZ;
+ vb.m_zr[1] = maxZ;
+ vb.m_zr[2] = minR;
+ vb.m_zr[3] = maxR;
+ m_vBoundaries.push_back(vb);
+ }
+ }
+}
+
+void TrigInDetRoadPredictorTool::findDetectorElements(unsigned int layerID, const SearchInterval& searchArea,
+ std::vector<unsigned int>& vIDs, bool hasHit) const {
+
+ const float road_width_rz = hasHit ? m_min_rz_rw : m_max_rz_rw;
+ const float road_width_rphi = hasHit ? m_min_rphi_rw : m_max_rphi_rw;
+
+ float phi_test = searchArea.m_phi;
+ float phi_res = road_width_rphi/searchArea.m_r;
+
+ if(phi_res > m_max_phi_rw) phi_res = m_max_phi_rw;
+ if(phi_res < m_min_phi_rw) phi_res = m_min_phi_rw;
+
+ float phi_min = phi_test - phi_res;
+ float phi_max = phi_test + phi_res;
+
+ const auto& L = (*m_layerMap.find(layerID)).second;
+
+ if(L.m_mappingType != 2) { // primary index is Phi, secondary is Z or R
+
+ float rz_test_m = searchArea.getMinR();
+ float rz_test_p = searchArea.getMaxR();
+
+ if(L.m_mappingType == 0) {
+ rz_test_m = searchArea.getMinZ();
+ rz_test_p = searchArea.getMaxZ();
+ }
+
+ rz_test_m -= road_width_rz;
+ rz_test_p += road_width_rz;
+
+ for(int iSubL = 0; iSubL < L.m_nSubLayers;iSubL++) {
+
+ for(const auto& prim : L.m_colls[iSubL]) {
+
+ const auto& slice = prim.second;
+
+ float f1 = slice.m_minCoord;
+ float f2 = slice.m_maxCoord;
+
+ if(std::abs(f2 - f1) < M_PI) {
+ if(phi_max < f1) continue;
+ if(phi_min > f2) continue;
+ }
+ else {// +/- pi boundary
+ std::swap(f2, f1);
+ if(phi_test < 0 && phi_min > f1) continue;
+ if(phi_test > 0 && phi_max < f2) continue;
+ }
+
+ for(const auto& de : slice.m_vDE) {
+
+ float p1 = de.m_minBound;
+ float p2 = de.m_maxBound;
+
+ if(rz_test_p < p1) break;
+ if(rz_test_m > p2) continue;
+
+ if (! (rz_test_p < p1 || rz_test_m > p2)) vIDs.push_back(de.m_hash);
+ }
+ }
+ }
+ }
+ else { //Strip endcaps: primary R, secondary Phi
+
+ float r_test_m = searchArea.getMinR() - road_width_rz;
+ float r_test_p = searchArea.getMaxR() + road_width_rz;
+
+ for(int iSubL = 0; iSubL < L.m_nSubLayers;iSubL++) {
+
+ for(const auto& prim : L.m_colls[iSubL]) {
+
+ const auto& slice = prim.second;
+
+ float r1 = slice.m_minCoord;
+ float r2 = slice.m_maxCoord;
+
+ if(r_test_p < r1) break;
+ if(r_test_m > r2) continue;
+
+ if (! (r_test_p < r1 || r_test_m > r2)) {
+
+ for(const auto& de : slice.m_vDE) {
+
+ float f1 = de.m_minBound;
+ float f2 = de.m_maxBound;
+
+ if(f2 - f1 < M_PI) {
+ if(phi_max < f1) continue;
+ if(phi_min > f2) continue;
+ }
+ else {// +/- pi boundary
+ if(phi_test < 0 && phi_min > f1) continue;
+ if(phi_test > 0 && phi_max < f2) continue;
+ }
+ vIDs.push_back(de.m_hash);
+ }
+ }
+ }
+ }
+ }
+
+}
+
-int TrigInDetRoadPredictorTool::getRoad(const std::vector<const Trk::SpacePoint*>& seed, std::vector<const InDetDD::SiDetectorElement*>& road, const EventContext& ctx) const {
+
+int TrigInDetRoadPredictorTool::getRoad(const std::vector<const Trk::SpacePoint*>& seed,
+ std::vector<const InDetDD::SiDetectorElement*>& road,
+ const EventContext& ctx) const {
+
+ const float MAX_R = 1030.0;//detector envelope
+ const float MAX_Z = 3000.0;//detector envelope
+ const float maxCornerDist = 15.0;
+
//1. get magnetic field
MagField::AtlasFieldCache fieldCache;
@@ -51,17 +454,271 @@ int TrigInDetRoadPredictorTool::getRoad(const std::vector<const Trk::SpacePoint*
road.clear();
- unsigned int seedSize = seed.size();
+ unsigned int nSP = seed.size();
- if(seedSize < 3) return -2;
+ if(nSP < 3) return -2;
- //2. adding spacepoints' DEs
-
- for(unsigned int spIdx=0;spIdx<seedSize;spIdx++) {
- const Trk::PrepRawData* prd = seed.at(spIdx)->clusterList().first;
+ std::vector<unsigned int> seedHashes;
+
+ for(unsigned int spIdx=0;spIdx<nSP;spIdx++) {
+ const auto& sp = seed.at(spIdx);
+ const Trk::PrepRawData* prd = sp->clusterList().first;
const InDet::PixelCluster* pPixelHit = dynamic_cast<const InDet::PixelCluster*>(prd);
- road.push_back(pPixelHit->detectorElement());
+ unsigned int hash = pPixelHit->detectorElement()->identifyHash();
+ seedHashes.push_back(hash);
+ }
+
+ std::vector<std::array<float,2> > zr;
+ zr.resize(nSP+2);
+
+ for(unsigned int spIdx=0;spIdx<nSP;spIdx++) {
+ const auto& sp = seed.at(spIdx);
+ zr[spIdx+1][0] = sp->globalPosition().z();
+ zr[spIdx+1][1] = sp->globalPosition().perp();
+ }
+
+ //adding the first point at beamline
+
+ zr[0][0] = zr[1][0] - zr[1][1]*(zr[2][0]-zr[1][0])/(zr[2][1]-zr[1][1]);
+ zr[0][1] = 0;
+
+ //adding the last point at detector exit
+ float zlast = zr[nSP-1][0] + (MAX_R - zr[nSP-1][1])*(zr[nSP][0]-zr[nSP-1][0])/(zr[nSP][1]-zr[nSP-1][1]);
+ float rlast = MAX_R;
+
+ if (std::fabs(zlast) > MAX_Z) {
+ if(zlast > 0) zlast = MAX_Z;
+ else zlast = -MAX_Z;
+ rlast = zr[nSP-1][1] + (zlast - zr[nSP-1][0])*(zr[nSP][1]-zr[nSP-1][1])/(zr[nSP][0]-zr[nSP-1][0]);
+ }
+ zr[nSP+1][0] = zlast;
+ zr[nSP+1][1] = rlast;
+
+ std::map<unsigned int, SearchInterval> rzIntervals;
+
+ for(unsigned int k1 = 0;k1<zr.size()-1;k1++) {//loop over trajectory segments
+
+ unsigned int k2 = k1+1;
+
+ float z1 = zr[k1][0];
+ float r1 = zr[k1][1];
+ float z2 = zr[k2][0];
+ float r2 = zr[k2][1];
+ float dz21 = z2-z1;
+ float dr21 = r2-r1;
+ float L = std::sqrt(dz21*dz21 + dr21*dr21);
+ float invL = 1.0/L;
+ float sinF = dr21*invL;
+ float cosF = dz21*invL;
+
+ for(const auto& vb : m_vBoundaries) {
+
+ if (z1 < vb.m_zr[0] && z2 < vb.m_zr[0]) continue;
+ if (z1 > vb.m_zr[1] && z2 > vb.m_zr[1]) continue;
+ if (r1 < vb.m_zr[2] && r2 < vb.m_zr[2]) continue;
+ if (r1 > vb.m_zr[3] && r2 > vb.m_zr[3]) continue;
+
+ //corners:
+
+ float zc[4] = {vb.m_zr[0], vb.m_zr[1], vb.m_zr[1], vb.m_zr[0]};
+ float rc[4] = {vb.m_zr[2], vb.m_zr[2], vb.m_zr[3], vb.m_zr[3]};
+
+ int nUp(0), nDn(0);
+
+ float minDistances[4];
+
+ for (int ic=0;ic<4;ic++) {
+ minDistances[ic] = (rc[ic] - r1)*cosF - (zc[ic] - z1)*sinF;
+ }
+
+ float minH = std::abs(minDistances[0]);
+
+ for (int ic=0;ic<4;ic++) {
+ float h = minDistances[ic];
+ if (h <=0) nDn += 1;
+ else nUp += 1;
+ if(std::abs(h) < minH) minH = std::abs(h);
+ }
+
+ if (nUp == 4 || nDn == 4) {
+ if(minH > maxCornerDist) {//the closest corner is still too far
+ continue;
+ }
+ }
+
+ //search through layers inside the volume
+
+ for(auto lIdx : vb.m_layers) {
+ const auto& lb = m_lBoundaries.at(lIdx);
+ for(int s1=0;s1<lb.m_nVertices-1;s1++) {
+ int s2 = s1 + 1;
+
+ float h1 = (lb.m_r[s1] - r1)*cosF - (lb.m_z[s1] - z1)*sinF;
+ float h2 = (lb.m_r[s2] - r1)*cosF - (lb.m_z[s2] - z1)*sinF;
+ if (h1*h2 > 0) continue;
+ float l1 = (lb.m_z[s1] - z1)*cosF + (lb.m_r[s1] - r1)*sinF;
+ float l2 = (lb.m_z[s2] - z1)*cosF + (lb.m_r[s2] - r1)*sinF;
+ if (l1 < 0 && l2 < 0) continue;
+ if (l1 > L && l2 > L) continue;
+ float lx = (l1 + (l2-l1)*h1/(h1-h2))*invL;
+
+ if (lx < 0 || lx > 1) continue;
+ float zx = z2*lx + (1-lx)*z1;
+ float rx = r2*lx + (1-lx)*r1;
+
+ auto radItr = rzIntervals.find(lb.m_lay_id);
+
+ if(radItr == rzIntervals.end()) {
+ rzIntervals.insert(std::make_pair(lb.m_lay_id, SearchInterval(zx, rx)));
+ }
+ else {
+ (*radItr).second.addPoint(zx, rx);
+ }
+ }
+ }
+ }
}
- return 0;
+ //3. parabolic extrapolation in r-phi
+
+ unsigned int sp1Idx = 0;
+ unsigned int sp3Idx = nSP-1;
+ unsigned int sp2Idx = (sp3Idx+sp1Idx)/2;
+
+ //3a. Converting XY coords to UV coords
+
+ float uv_coords[2][2];
+
+ float dx = seed.at(sp3Idx)->globalPosition().x() - seed.at(sp2Idx)->globalPosition().x();
+ float dy = seed.at(sp3Idx)->globalPosition().y() - seed.at(sp2Idx)->globalPosition().y();
+
+ uv_coords[1][0] = -std::sqrt(dx*dx + dy*dy);
+ uv_coords[1][1] = 0.0;
+
+ float cos_theta = dx/(-uv_coords[1][0]);
+ float sin_theta = dy/(-uv_coords[1][0]);
+
+ float rot_matrix[2][2];
+ float inv_rot_matrix[2][2];
+
+ rot_matrix[0][0] = cos_theta;
+ rot_matrix[0][1] = sin_theta;
+ rot_matrix[1][0] = -sin_theta;
+ rot_matrix[1][1] = cos_theta;
+
+ inv_rot_matrix[0][0] = cos_theta;
+ inv_rot_matrix[0][1] = -sin_theta;
+ inv_rot_matrix[1][0] = sin_theta;
+ inv_rot_matrix[1][1] = cos_theta;
+
+ float sp3_coords[3];
+ sp3_coords[0] = seed.at(sp3Idx)->globalPosition().x();
+ sp3_coords[1] = seed.at(sp3Idx)->globalPosition().y();
+ sp3_coords[2] = seed.at(sp3Idx)->globalPosition().z();
+
+ //UV-coordinates of the XY origin (0,0)
+
+ float u_c = rot_matrix[0][0]*(0-sp3_coords[0]) + rot_matrix[0][1]*(0-sp3_coords[1]);
+ float v_c = rot_matrix[1][0]*(0-sp3_coords[0]) + rot_matrix[1][1]*(0-sp3_coords[1]);
+
+ int sign_up = u_c < 0 ? 1 : -1;
+
+ //transforming SP1
+
+ float dR1[2];
+
+ dR1[0] = seed.at(sp1Idx)->globalPosition().x() - sp3_coords[0];
+ dR1[1] = seed.at(sp1Idx)->globalPosition().y() - sp3_coords[1];
+
+ uv_coords[0][0] = rot_matrix[0][0]*dR1[0] + rot_matrix[0][1]*dR1[1];
+ uv_coords[0][1] = rot_matrix[1][0]*dR1[0] + rot_matrix[1][1]*dR1[1];
+
+ //3b. parameters of the parabola in the u-v c.s.
+
+ float a = uv_coords[0][1]/(uv_coords[0][0]*(uv_coords[0][0]-uv_coords[1][0]));
+ float b = -a*uv_coords[1][0]; // b = -a*x2
+
+ //3c. calculate impact point radii
+
+ std::vector<unsigned int> pixelHashIds;
+ std::vector<unsigned int> stripHashIds;
+
+ for(auto& ip : rzIntervals) {
+
+ float R = ip.second.getAverageRadius();
+
+ float R2 = R*R;
+
+ float dRv = R2-v_c*v_c;
+
+ if(dRv < 0) continue;//no intersection with the circle
+
+ float u_p = u_c + sign_up*std::sqrt(dRv);
+
+ float v_p = b*u_p + a*u_p*u_p;
+ float dv2 = (v_p-v_c)*(v_p-v_c);
+
+ if (R2-dv2<0) continue; // check for intersection between v=v_p and the circle
+
+ float u_star = u_c + sign_up*std::sqrt(R2-dv2);
+ float v_star = b*u_star + a*u_star*u_star;
+
+ float x_star = inv_rot_matrix[0][0]*u_star + inv_rot_matrix[0][1]*v_star + sp3_coords[0];
+ float y_star = inv_rot_matrix[1][0]*u_star + inv_rot_matrix[1][1]*v_star + sp3_coords[1];
+
+ ip.second.m_r = std::sqrt(x_star*x_star + y_star*y_star);
+ ip.second.m_phi = std::atan2(y_star, x_star);
+
+ if(ip.first > 15000) {//assuming Pixel-only seeds
+
+ bool hasHit = false;
+
+ for(unsigned int i=1;i<nSP-1;i++) {
+
+ float dpr = ip.second.m_r - zr[i][1];
+ float dpz = ip.second.m_z - zr[i][0];
+ float dist = std::sqrt(dpr*dpr + dpz*dpz);
+ if(dist < maxCornerDist) {
+ hasHit = true;
+ break;
+ }
+ }
+ findDetectorElements(ip.first, ip.second, pixelHashIds, hasHit);
+ }
+ else {
+ findDetectorElements(ip.first, ip.second, stripHashIds, false);
+ }
+ }
+
+ std::set<unsigned int> pixelHashSet(pixelHashIds.begin(), pixelHashIds.end());
+
+ for(auto id : seedHashes) {
+ if(pixelHashSet.find(id) == pixelHashSet.end()) pixelHashIds.push_back(id);
+ }
+
+ std::vector<std::pair<float, const InDetDD::SiDetectorElement*> > theRoad;
+
+ for(auto hash_id : pixelHashIds) {
+ const InDetDD::SiDetectorElement *p = m_pixelManager->getDetectorElement(hash_id);
+ if(p == nullptr) continue;
+ const Amg::Vector3D& C = p->center();
+ float dist = std::sqrt(C(0)*C(0) + C(1)*C(1) + C(2)*C(2));
+ theRoad.push_back(std::make_pair(dist,p));
+ }
+
+ for(auto hash_id : stripHashIds) {
+ const InDetDD::SiDetectorElement *p = m_stripManager->getDetectorElement(hash_id);
+ if(p == nullptr) continue;
+ const Amg::Vector3D& C = p->center();
+ float dist = std::sqrt(C(0)*C(0) + C(1)*C(1) + C(2)*C(2));
+ theRoad.push_back(std::make_pair(dist,p));
+ }
+
+ std::sort(theRoad.begin(), theRoad.end());
+
+ for(const auto & dp : theRoad) {
+ road.push_back(dp.second);
+ }
+
+ return (int)theRoad.size();
}
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.h b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.h
index 0796ba21f6c4db6b09e7b67f806c073c36d5d734..9ee21f1a9fce999eff8d5e57ba1849478c9b3904 100644
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.h
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetRoadPredictorTool.h
@@ -10,31 +10,147 @@
#include "GaudiKernel/ServiceHandle.h"
#include "TrigInDetToolInterfaces/ITrigInDetRoadPredictorTool.h"
+#include "TrigInDetToolInterfaces/ITrigL2LayerNumberTool.h"
// MagField cache
#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
#include "MagFieldElements/AtlasFieldCache.h"
+#include <fstream>
+#include <array>
+#include <map>
+
+namespace InDetDD {
+ class PixelDetectorManager;
+ class SCT_DetectorManager;
+}
+
namespace Trk {
class SpacePoint;
}
+class PixelID;
+class SCT_ID;
+
class TrigInDetRoadPredictorTool: public AthAlgTool, virtual public ITrigInDetRoadPredictorTool
{
public:
TrigInDetRoadPredictorTool( const std::string&, const std::string&, const IInterface* );
+
virtual StatusCode initialize() override;
- virtual int getRoad(const std::vector<const Trk::SpacePoint*>&, std::vector<const InDetDD::SiDetectorElement*>&, const EventContext&) const override;
+ virtual int getRoad(const std::vector<const Trk::SpacePoint*>&, std::vector<const InDetDD::SiDetectorElement*>&,
+ const EventContext&) const override;
private:
+
+ struct DetectorElementDescription {
+ DetectorElementDescription(unsigned int h) : m_hash(h) {};
+ unsigned int m_hash;
+ short m_index;
+ float m_ref;//z or Phi
+ float m_c[4][3];
+ float m_minBound, m_maxBound;
+ };
+
+ struct DetectorElementsCollection {
+ DetectorElementsCollection(short idx, float min, float max) : m_index(idx), m_minCoord(min), m_maxCoord(max) {};
+ short m_index;
+ float m_minCoord, m_maxCoord; //phi or R
+ std::vector<DetectorElementDescription> m_vDE;
+ };
+
+ struct LayerDescription {
+ LayerDescription(unsigned int id, int nSL, unsigned int mType) : m_id(id), m_nSubLayers(nSL), m_mappingType(mType) {};
+ unsigned int m_id;
+ int m_nSubLayers;
+ unsigned int m_mappingType;
+ std::map<short, DetectorElementsCollection> m_colls[2];
+ };
+
+
+ struct SearchInterval {
+ SearchInterval(float z, float r) : m_r(0.0), m_z(0.0), m_phi(0.0) {
+ m_vZx.push_back(z);
+ m_vRx.push_back(r);
+ }
+
+ void addPoint(float z, float r) {
+ m_vZx.push_back(z);
+ m_vRx.push_back(r);
+ }
+
+ float getAverageRadius() {
+ return std::accumulate(m_vRx.begin(), m_vRx.end(), 0)/m_vRx.size();
+ }
+
+ float getAverageZ() {
+ return std::accumulate(m_vZx.begin(), m_vZx.end(), 0)/m_vZx.size();
+ }
+
+ float getMinZ() const {
+ return *std::min_element( m_vZx.begin(), m_vZx.end());
+ }
+
+ float getMaxZ() const {
+ return *std::max_element( m_vZx.begin(), m_vZx.end());
+ }
+
+ float getMinR() const {
+ return *std::min_element( m_vRx.begin(), m_vRx.end());
+ }
+
+ float getMaxR() const {
+ return *std::max_element( m_vRx.begin(), m_vRx.end());
+ }
+
+ std::vector<float> m_vZx;
+ std::vector<float> m_vRx;
+ float m_r;
+ float m_z;
+ float m_phi;
+ };
- Gaudi::Property<int> m_nClustersMin {this, "nClustersMin", 7, "Minimum number of clusters on track"};
+ struct VolumeBoundary {
+ int m_index;
+ float m_zr[4];
+ int m_vol_id;
+ std::vector<int> m_layers;
+ };
+ struct LayerBoundary {
+ int m_index, m_lay_id, m_nVertices;
+ std::vector<float> m_z;
+ std::vector<float> m_r;
+ };
+
+ void buildDetectorDescription();
+ void createHitBoxes();
+
+ void addNewElement(unsigned int, short, short, const InDetDD::SiDetectorElement*);
+ void findDetectorElements(unsigned int, const SearchInterval&, std::vector<unsigned int>&, bool) const;
+
+ ToolHandle<ITrigL2LayerNumberTool> m_layerNumberTool {this, "LayerNumberTool", "TrigL2LayerNumberToolITk"};
+
+ Gaudi::Property<float> m_min_rz_rw {this, "MinRzRoadWidth", 3.0, "Minimum rz road width"};
+ Gaudi::Property<float> m_max_rz_rw {this, "MaxRzRoadWidth", 15.0, "Maximum rz road width"};
+ Gaudi::Property<float> m_min_rphi_rw {this, "MinRPhiRoadWidth", 3.0, "Minimum rphi road width"};
+ Gaudi::Property<float> m_max_rphi_rw {this, "MaxRPhiRoadWidth", 15.0, "Maximum rphi road width"};
+ Gaudi::Property<float> m_min_phi_rw {this, "MinPhiRoadWidth", 0.05, "Minimum phi road width"};
+ Gaudi::Property<float> m_max_phi_rw {this, "MaxPhiRoadWidth", 0.1, "Maximum phi road width"};
+
SG::ReadCondHandleKey<AtlasFieldCacheCondObj> m_fieldCondObjInputKey{
this, "AtlasFieldCacheCondObj", "fieldCondObj",
"Name of the Magnetic Field conditions object key"};
-
+
+ std::vector<VolumeBoundary> m_vBoundaries;
+ std::vector<LayerBoundary> m_lBoundaries;
+ std::map<unsigned int, LayerDescription> m_layerMap;
+
+ const InDetDD::PixelDetectorManager* m_pixelManager{nullptr};
+ const InDetDD::SCT_DetectorManager* m_stripManager{nullptr};
+ const PixelID* m_pixelId{nullptr};
+ const SCT_ID* m_stripId{nullptr};
};
#endif
diff --git a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFollowingTool.cxx b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFollowingTool.cxx
index baff86aa16cc2ccaeefde2ccd0ff4723f62677d2..4274ff26280fec60adff5a1b87847204fe599d6e 100644
--- a/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFollowingTool.cxx
+++ b/Trigger/TrigTools/TrigInDetTrackFitter/src/TrigInDetTrackFollowingTool.cxx
@@ -261,7 +261,7 @@ double TrigInDetTrackFollowingTool::processHit(const InDet::SCT_Cluster* pPRD, i
const Trk::PrepRawData* TrigInDetTrackFollowingTool::updateTrackState(const InDet::PixelCluster* pInputHit, const InDet::PixelClusterCollection* pColl, TrigFTF_ExtendedTrackState& ets) const {
const InDet::PixelCluster* bestHit = pInputHit;
-
+ if(pColl == nullptr && pInputHit == nullptr) return nullptr;
double resid[2];
double invcov[3];
@@ -317,7 +317,7 @@ const Trk::PrepRawData* TrigInDetTrackFollowingTool::updateTrackState(const InDe
const Trk::PrepRawData* TrigInDetTrackFollowingTool::updateTrackState(const InDet::SCT_Cluster* pInputHit, const InDet::SCT_ClusterCollection* pColl, int shape, TrigFTF_ExtendedTrackState& ets) const {
const InDet::SCT_Cluster* bestHit = pInputHit;
-
+ if(pColl == nullptr && pInputHit == nullptr) return nullptr;
double resid, invcov;
double H[2];//linearized observation matrix