ATLAS Offline Software
FPGATrackSimLLPDoubletHoughTransformTool.cxx
Go to the documentation of this file.
1 // Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
2 
3 #include "TH2.h"
5 
6 // simple 2D vector with basic arithmetics
7 typedef double fp_t; // floating point type alias, can change it to float from double
8 
9 typedef std::array<fp_t, 2> pvec;
10 pvec operator+(const pvec& a, const pvec& b);
11 pvec operator-(const pvec& a, const pvec& b);
12 pvec operator*(const pvec& a, double scale);
13 
14 
15 pvec operator+(const pvec& a, const pvec& b) {
16  return { {a[0]+b[0], a[1]+b[1]} };
17 }
18 pvec operator-(const pvec& a, const pvec& b) {
19  return { {a[0]-b[0], a[1]-b[1]} };
20 
21 }
22 pvec operator*(const pvec& a, double scale) {
23  return { {a[0]*scale, a[1]*scale} };
24 }
25 
26 double length(const pvec& v) {
27  return std::hypot(v[0], v[1]);
28 }
29 
30 pvec rotate90( const pvec& v) {
31  return { -v[1], v[0]};
32 }
33 
34 double crossProduct( const pvec& a, const pvec& b ) {
35  return a[0]*b[1] - a[1]*b[0];
36 }
37 
38 
39 FPGATrackSimLLPDoubletHoughTransformTool::FPGATrackSimLLPDoubletHoughTransformTool(const std::string& algname, const std::string& name , const IInterface* ifc)
40  : base_class(algname, name, ifc)
41 {
42  declareInterface<IFPGATrackSimRoadFinderTool>(this);
43 }
44 
45 
47  if (m_imageSize_y %2 == 1) {
48  ATH_MSG_ERROR("Can not have odd number " << m_imageSize_y << " of bins in q/pT - will result in division by 0");
49  return StatusCode::FAILURE;
50  }
53 
54  return StatusCode::SUCCESS;
55  }
56 
58  ATH_MSG_INFO("Number of events procesed by the FPGATrackSimLLPDoubletHoughTransformTool " << m_eventsProcessed << " roads " << m_roadsGenerated);
59  return StatusCode::SUCCESS;
60 }
61 
62 StatusCode FPGATrackSimLLPDoubletHoughTransformTool::getRoads(const std::vector<std::shared_ptr<const FPGATrackSimHit>> & hits, std::vector<std::shared_ptr<const FPGATrackSimRoad>> & roads) {
63  roads.clear();
64  m_roads.clear();
66  Image image(m_imageSize_x, m_imageSize_y); // not quite optimal to allocate this memory in evey event, but following nominal HT
67  for (unsigned ihit1 = 0; ihit1 < hits.size(); ihit1++) {
68  std::shared_ptr<const FPGATrackSimHit> hit1 = hits[ihit1];
69  for (unsigned ihit2 = ihit1+1; ihit2 < hits.size(); ihit2++) {
70  std::shared_ptr<const FPGATrackSimHit> hit2 = hits[ihit2];
71  ATH_MSG_DEBUG("Hits pair R: " << hit1->getR() << " " << hit2->getR());
72  //TODO: replace with qualification by layer IDs
73  const double radiusDifference = hit2->getR() - hit1->getR();
74 
75  if ( not (m_acceptedDistanceBetweenLayersMin < radiusDifference && radiusDifference < m_acceptedDistanceBetweenLayersMax) )
76  continue;
77  if ( hit1->getLayer() == hit2->getLayer() )
78  continue;
79  ATH_CHECK(fillImage( hit1, hit2, image ));
80  }
81  }
82  if ( m_event < 20 )
84 
85  for (unsigned y = 0; y < m_imageSize_y; y++) {
86  for (unsigned x = 0; x < m_imageSize_x; x++) {
87  if (passThreshold(image, x, y) && isLocalMaxima( image, x, y ) ) {
88  ATH_MSG_DEBUG("Creating road with threshold " << image(x,y).first << " hits " << image(x,y).second.size()
89  << " at x: " << x << " d0 " << xtod0(x) << " y: "<< y << " q/pt " << ytoqoverpt(y));
91  m_roads.push_back(createRoad(image(x, y).second, x, y));
92  }
93  }
94  }
95  roads.reserve(m_roads.size());
96  for (FPGATrackSimRoad & r : m_roads) roads.emplace_back(std::make_shared<const FPGATrackSimRoad>(r));
97  m_event++;
98  return StatusCode::SUCCESS;
99 }
100 
101 
102 StatusCode FPGATrackSimLLPDoubletHoughTransformTool::fillImage(std::shared_ptr<const FPGATrackSimHit> &hit1, std::shared_ptr<const FPGATrackSimHit> &hit2, Image& image) const {
103  const pvec p1 {{hit1->getX(), hit1->getY()}};
104  const pvec p2 {{hit2->getX(), hit2->getY()}};
105  const pvec halfDiff = (p2 - p1)*0.5;
106  const fp_t halfLen = length(halfDiff);
107 
108  int xbefore = -1;
109 
110  for ( unsigned y = 0; y < m_imageSize_y; y++ ) {
111  const fp_t qoverpt = (y * m_step_y) + m_step_y*0.5 - m_qOverPt_range; // use middle f the bin
112  const fp_t radius = 1.0/(0.6*qoverpt); // magic transform constants
113  const fp_t scale = std::copysign( std::sqrt( std::pow(radius/halfLen, 2) - 1), radius );
114  const pvec rprime = rotate90(halfDiff) * scale;
115  const pvec center = p1 + halfDiff + rprime;
116  const fp_t d0 = (std::signbit(radius) ? -1.0 : 1.0)*(length(center) - abs(radius));
117 
118  // find the bin x position in the image
119  int x = (d0 + m_d0_range) / m_step_x;
120 
121  if ( 0 <= x && x < (int)m_imageSize_x) {
122  if (xbefore == -1) xbefore = x;
123  if ( m_continuous ) { // fill the bins along x starintg from the last one filled
124  const int xmin = (xbefore < x)? xbefore: x;
125  const int xmax = (xbefore < x)? x: xbefore;
126  for ( int xinterpolated = xmin; xinterpolated <= xmax; ++xinterpolated) {
127  image(xinterpolated, y).first++;
128  image(xinterpolated, y).second.insert( hit1 );
129  image(xinterpolated, y).second.insert( hit2 );
130  }
131  } else {
132  image(x, y).first++;
133  image(x, y).second.insert( hit1 );
134  image(x, y).second.insert( hit2 );
135  }
136  xbefore = x;
137  }
138  }
139  return StatusCode::SUCCESS;
140 }
141 
142 FPGATrackSimRoad FPGATrackSimLLPDoubletHoughTransformTool::createRoad(const std::unordered_set<std::shared_ptr<const FPGATrackSimHit>> & hits, unsigned /*x*/, unsigned /*y*/) const {
143  // Get the road hits
144  std::vector<std::shared_ptr<const FPGATrackSimHit>> road_hits;
145  layer_bitmask_t hitLayers = 0;
146  for (const auto & hit : hits)
147  {
148  road_hits.push_back(hit);
149  hitLayers |= 1 << hit->getLayer();
150  }
151 
152  auto sorted_hits = ::sortByLayer(road_hits);
153  sorted_hits.resize(8); // If no hits in last layer, return from sortByLayer will be too short
154 
156  r.setHitLayers(hitLayers);
157  r.setHits(std::move(sorted_hits));
158  return r;
159 }
160 
161 
163  const int count = image(x,y).first;
164  const float d0 = xtod0(x);
165  if ( std::abs(d0) < 50.0 && count >= m_threshold50 ) return true;
166  if ( std::abs(d0) >= 50.0 && count >= m_threshold ) return true;
167 
168  return false;
169 }
170 
172  const auto centerValue = image(x,y).first;
173  for ( int xaround = std::min(x-1, 0); xaround <= std::max((int)m_imageSize_x-1, x+1); xaround++ ) {
174  for ( int yaround = std::min(y-1, 0); yaround <= std::max((int)m_imageSize_y-1, y+1); yaround++ ) {
175  if ( image(xaround,yaround).first > centerValue )
176  return false;
177  }
178  }
179  return true;
180 }
181 
183 
184  TH2I h(("event"+std::to_string( m_event )).c_str(), "LLP Doublet Hough space;d_{0}[mm];q/pT [e/GeV]",
187 
188  for (unsigned x = 0; x < m_imageSize_x; x++)
189  for (unsigned y = 0; y < m_imageSize_y; y++)
190  h.SetBinContent(x+1, y+1, image(x, y).first); // +1 since root bins are 1-indexed
191 
192  h.SaveAs((name()+"_event_"+std::to_string( m_event )+".root").c_str(), "");
193 }
194 
beamspotman.r
def r
Definition: beamspotman.py:676
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
getMenu.algname
algname
Definition: getMenu.py:54
FPGATrackSimLLPDoubletHoughTransformTool::xtod0
float xtod0(int x) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:67
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
operator+
pvec operator+(const pvec &a, const pvec &b)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:15
FPGATrackSimLLPDoubletHoughTransformTool::passThreshold
bool passThreshold(Image const &image, int x, int y) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:162
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
FPGATrackSimLLPDoubletHoughTransformTool::m_eventsProcessed
int m_eventsProcessed
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:69
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
FPGATrackSimHit::getX
float getX() const
Definition: FPGATrackSimHit.h:137
crossProduct
double crossProduct(const pvec &a, const pvec &b)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:34
FPGATrackSimLLPDoubletHoughTransformTool::m_acceptedDistanceBetweenLayersMin
const double m_acceptedDistanceBetweenLayersMin
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:53
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
FPGATrackSimHit::getLayer
unsigned getLayer() const
Definition: FPGATrackSimHit.cxx:77
FPGATrackSimLLPDoubletHoughTransformTool::m_continuous
Gaudi::Property< bool > m_continuous
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:47
FPGATrackSimLLPDoubletHoughTransformTool::createRoad
FPGATrackSimRoad createRoad(const std::unordered_set< std::shared_ptr< const FPGATrackSimHit >> &hits, unsigned x, unsigned y) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:142
FPGATrackSimLLPDoubletHoughTransformTool::ytoqoverpt
float ytoqoverpt(int y) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:68
FPGATrackSimLLPDoubletHoughTransformTool::getRoads
virtual StatusCode getRoads(const std::vector< std::shared_ptr< const FPGATrackSimHit >> &hits, std::vector< std::shared_ptr< const FPGATrackSimRoad >> &roads) override
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:62
FPGATrackSimLLPDoubletHoughTransformTool::m_imageSize_y
Gaudi::Property< unsigned > m_imageSize_y
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:41
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
FPGATrackSimLLPDoubletHoughTransformTool::FPGATrackSimLLPDoubletHoughTransformTool
FPGATrackSimLLPDoubletHoughTransformTool(const std::string &, const std::string &, const IInterface *)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:39
x
#define x
XMLtoHeader.count
count
Definition: XMLtoHeader.py:85
FPGATrackSimLLPDoubletHoughTransformTool::m_imageSize_x
Gaudi::Property< unsigned > m_imageSize_x
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:40
fp_t
double fp_t
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:7
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
FPGATrackSimLLPDoubletHoughTransformTool::initialize
virtual StatusCode initialize() override
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:46
vector2D< std::pair< int, std::unordered_set< std::shared_ptr< const FPGATrackSimHit > > > >
FPGATrackSimLLPDoubletHoughTransformTool::m_threshold
Gaudi::Property< int > m_threshold
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:42
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
xmin
double xmin
Definition: listroot.cxx:60
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
extractSporadic.h
list h
Definition: extractSporadic.py:97
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
FPGATrackSimLLPDoubletHoughTransformTool::m_acceptedDistanceBetweenLayersMax
const double m_acceptedDistanceBetweenLayersMax
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:54
FPGATrackSimLLPDoubletHoughTransformTool::m_roadsGenerated
int m_roadsGenerated
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:70
TRT::Track::d0
@ d0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:62
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
FPGATrackSimHit::getY
float getY() const
Definition: FPGATrackSimHit.h:138
FPGATrackSimLLPDoubletHoughTransformTool::m_d0_range
Gaudi::Property< float > m_d0_range
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:45
operator*
pvec operator*(const pvec &a, double scale)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:22
FPGATrackSimLLPDoubletHoughTransformTool::drawImage
void drawImage(Image const &image) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:182
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
ActsTrk::to_string
std::string to_string(const DetectorType &type)
Definition: GeometryDefs.h:34
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
MyPlots.image
string image
Definition: MyPlots.py:43
FPGATrackSimLLPDoubletHoughTransformTool::m_event
unsigned m_event
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:64
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
FPGATrackSimLLPDoubletHoughTransformTool::m_threshold50
Gaudi::Property< int > m_threshold50
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:43
FPGATrackSimLLPDoubletHoughTransformTool::m_step_x
double m_step_x
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:56
python.PyAthena.v
v
Definition: PyAthena.py:154
FPGATrackSimLLPDoubletHoughTransformTool::finalize
virtual StatusCode finalize() override
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:57
FPGATrackSimLLPDoubletHoughTransformTool::m_qOverPt_range
Gaudi::Property< float > m_qOverPt_range
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:46
FPGATrackSimHit::getR
float getR() const
Definition: FPGATrackSimHit.h:140
a
TList * a
Definition: liststreamerinfos.cxx:10
y
#define y
h
layer_bitmask_t
uint32_t layer_bitmask_t
Definition: FPGATrackSimTypes.h:22
DeMoScan.first
bool first
Definition: DeMoScan.py:536
FPGATrackSimLLPDoubletHoughTransformTool::fillImage
StatusCode fillImage(std::shared_ptr< const FPGATrackSimHit > &hit1, std::shared_ptr< const FPGATrackSimHit > &hit2, Image &image) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:102
xmax
double xmax
Definition: listroot.cxx:61
FPGATrackSimLLPDoubletHoughTransformTool::m_roads
std::vector< FPGATrackSimRoad > m_roads
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:51
rotate90
pvec rotate90(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:30
operator-
pvec operator-(const pvec &a, const pvec &b)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:18
pvec
std::array< fp_t, 2 > pvec
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:9
FPGATrackSimLLPDoubletHoughTransformTool.h
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
FPGATrackSimLLPDoubletHoughTransformTool::isLocalMaxima
bool isLocalMaxima(Image const &image, int x, int y) const
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:171
length
double length(const pvec &v)
Definition: FPGATrackSimLLPDoubletHoughTransformTool.cxx:26
FPGATrackSimRoad
Definition: FPGATrackSimRoad.h:30
FPGATrackSimLLPDoubletHoughTransformTool::m_step_y
double m_step_y
Definition: FPGATrackSimLLPDoubletHoughTransformTool.h:57
sortByLayer
std::vector< std::vector< std::shared_ptr< const FPGATrackSimHit > > > sortByLayer(Container const &hits)
Definition: FPGATrackSimHit.h:270