ATLAS Offline Software
MuonLayerHough.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
7 #include <TH1.h>
8 #include <TString.h> // for Form
9 #include <memory.h>
10 
11 #include <bitset>
12 #include <cmath>
13 #include <cstdlib>
14 #include <iostream>
15 
17 #include "GaudiKernel/MsgStream.h"
18 #include "GaudiKernel/SystemOfUnits.h"
19 namespace {
20  // Avoid floating point exceptions arising from cases of x = 0 or x = PI
21  // by extending the inverse tan function towards a large number
22  // Avoid FPEs occuring in clang 10 , e.g.,
23  // FPEAuditor 2 1 WARNING FPE OVERFLOW in [Execute] of [MuGirlStauAlg] on event 257948896 0 0
24  // by large number
26  float cot(const float x) {
27  const float arg = M_PI_2 - x;
28  // Tan becomes singular at -pi/2 and pi/2
29  if (std::abs(arg - M_PI_2) <= FLT_EPSILON || std::abs(arg + M_PI_2) <= FLT_EPSILON) return 1.e8;
30  return std::tan(arg);
31  }
32 } // namespace
33 namespace MuonHough {
34 
36  m_descriptor(descriptor) {
37  // calculate the number of bins
38  m_nbins = (m_descriptor.yMaxRange - m_descriptor.yMinRange) / m_descriptor.yBinSize; // the same for all the cycles
40  m_invbinsize = 1. / m_binsize; // binsize in the floating plane
41 
42  // setup the histograms; determines the number of theta slices
43  m_histos.reserve(m_descriptor.nthetaSamples); // it contains all the theta layers
44  for (unsigned int i = 0; i < m_descriptor.nthetaSamples; ++i) m_histos.emplace_back(new unsigned int[m_nbins]);
45  reset();
46  }
47 
49  for (unsigned int i = 0; i < m_histos.size(); ++i) memset(m_histos[i].get(), 0, sizeof(unsigned int) * m_nbins);
50  max = 0;
51  maxhist = -1;
52  maxbin = -1;
53  }
54 
55  void MuonLayerHough::fill(float x, float y, float weight) {
56  int cycles = m_histos.size();
57  for (int ci = 0; ci < cycles; ++ci) {
58  float dtheta = m_descriptor.thetaStep;
59  float dthetaOffset = 2 * m_descriptor.thetaStep * (ci - (cycles - 1) / 2.);
60  float theta = std::atan2(x, y);
61  float zref = (m_descriptor.referencePosition - x) * cot(theta - dthetaOffset) + y;
62  float z0 = (m_descriptor.referencePosition - x) * cot(theta - dthetaOffset + dtheta) + y;
63  float z1 = (m_descriptor.referencePosition - x) * cot(theta - dthetaOffset - dtheta) + y;
64 
65  float zmin = z0 < z1 ? z0 : z1;
66  float zmax = z0 < z1 ? z1 : z0;
67  int bincenter = (zref - m_descriptor.yMinRange) * m_invbinsize;
68 
69  int binmin = (zmin - m_descriptor.yMinRange) * m_invbinsize;
70  int binmax = (zmax - m_descriptor.yMinRange) * m_invbinsize;
71  if (binmin - bincenter < -3) binmin = bincenter - 3;
72  if (binmin == bincenter) binmin -= 1;
73  if (binmax == bincenter) binmax += 1;
74  if (binmin >= m_nbins) continue;
75  if (binmax - bincenter > 3) binmax = bincenter + 3;
76  if (binmax < 0) continue;
77 
78  if (binmin < 0) binmin = 0;
79  if (binmax >= m_nbins) binmax = m_nbins - 1;
80  if (m_debug)
81  std::cout << " filling " << x << " y " << binmin << " " << binmax << " w " << weight << " center " << bincenter << " range "
83  << std::endl;
84  for (int n = binmin; n <= binmax; ++n) {
85  unsigned int& val = m_histos[ci][n];
86  int w = 1000 * weight;
87  if (w < 0 && (int)val < -w)
88  val = 0;
89  else
90  val += weight;
91  if (val > max) {
92  max = val;
93  maxhist = ci;
94  maxbin = n;
95  }
96  }
97  }
98  }
99 
100  void MuonLayerHough::fillLayer(const HitVec& hits, bool subtract) {
101  if (hits.empty()) return;
102 
103  // outer loop over cycles
104  int cycles = m_histos.size();
105  for (int ci = 0; ci < cycles; ++ci) {
106  // float dtheta = m_descriptor.thetaStep;
107  // float dthetaOffset = 2*m_descriptor.thetaStep*(ci-(cycles-1)/2.);
108 
109  int prevlayer = -1;
110  int prevbinmin = 10000;
111  int prevbinmax = -1;
112  // inner loop over hits
113  HitVec::const_iterator it = hits.begin();
114  HitVec::const_iterator it_end = hits.end();
115  for (; it != it_end; ++it) {
116  float x = (*it)->x;
117  float y1 = (*it)->ymin;
118  float y2 = (*it)->ymax;
119  std::pair<int, int> minMax = range((*it)->x, (*it)->ymin, (*it)->ymax, ci);
120  int binmin = std::max(minMax.first,0);
121  int binmax = std::min(minMax.second, m_nbins - 1);
122  if (binmin >= m_nbins || binmax < 0) continue;
123 
124  if (m_debug) {
125  std::cout << " filling hit " << x << " refpos " << m_descriptor.referencePosition << " ymin " << y1 << " ymax " << y2
126  << " layer " << (*it)->layer << " binmin " << binmin << " max " << binmax;
127  if ((*it)->debugInfo()) {
128  const HitDebugInfo* db1 = (*it)->debugInfo();
129  std::cout << " sec " << db1->sector << " r " << db1->region << " type " << db1->type << " lay " << db1->layer
130  << " slay " << db1->sublayer << std::endl;
131  } else
132  std::cout << std::endl;
133  }
134  // first hit within range
135  if (prevbinmax == -1) {
136  if (m_debug) std::cout << " first range " << binmin << " " << binmax << std::endl;
137  prevbinmin = binmin;
138  prevbinmax = binmax;
139  prevlayer = (*it)->layer;
140  continue;
141  }
142 
143  if (binmin < prevbinmin && prevlayer == (*it)->layer) {
144  std::cout << "Error hits are out of order: min " << binmin << " max " << binmax << " lay " << (*it)->layer << std::endl;
145  }
146  // if the max value of the previous hit is smaller than the current minvalue fill the histogram of the previous hit
147  // do the same when reached last hit
148  if (prevbinmax < binmin || prevlayer != (*it)->layer) {
149  if (m_debug)
150  std::cout << " filling range " << prevbinmin << " " << prevbinmax << " new min " << binmin << " " << binmax
151  << " weight " << (*it)->w << std::endl;
152  for (int n = prevbinmin; n <= prevbinmax; ++n) {
153  unsigned int& val = m_histos[ci][n];
154  int w = 1000 * (*it)->w;
155  if (subtract) w *= -1;
156  if (w < 0 && (int)val < -w)
157  val = 0;
158  else
159  val += w;
160  if (val > max) {
161  max = val;
162  maxhist = ci;
163  maxbin = n;
164  }
165  }
166  prevbinmin = binmin;
167  prevbinmax = binmax;
168  prevlayer = (*it)->layer;
169 
170  } else {
171  // update the maximum value of the window
172  if (m_debug)
173  std::cout << " updating range " << prevbinmin << " " << prevbinmax << " hit " << binmin << " " << binmax
174  << " new " << prevbinmin << " " << binmax << std::endl;
175  prevbinmax = binmax;
176  }
177  }
178  if (prevbinmax != -1) {
179  if (m_debug)
180  std::cout << " filling last range " << prevbinmin << " " << prevbinmax << " weight " << hits.back()->w << std::endl;
181  for (int n = prevbinmin; n <= prevbinmax; ++n) {
182  unsigned int& val = m_histos[ci][n];
183  int w = 1000 * hits.back()->w;
184  if (subtract) w *= -1;
185  if (w < 0 && (int)val < -w)
186  val = 0;
187  else
188  val += w;
189  if (val > max) {
190  max = val;
191  maxhist = ci;
192  maxbin = n;
193  }
194  }
195  }
196  }
197  }
198 
199  void MuonLayerHough::fillLayer2(const HitVec& hits, bool subtract) {
200  if (hits.empty()) return;
201 
202  std::vector<int> layerCounts(m_nbins, 0);
203  int sign = subtract ? -1000 : 1000;
204  // outer loop over cycles
205  int cycles = m_histos.size();
206  for (int ci = 0; ci < cycles; ++ci) {
207  // keep track of the previous layer
208  int prevlayer = hits.front()->layer;
209 
210  // inner loop over hits
211  HitVec::const_iterator it = hits.begin();
212  HitVec::const_iterator it_end = hits.end();
213  for (; it != it_end; ++it) {
214  // if we get to the next layer process the current one and fill the Hough space
215  if (prevlayer != (*it)->layer) {
216  for (int i = 0; i < m_nbins; ++i) {
217  if (subtract && -layerCounts[i] >= static_cast<int>(m_histos[ci][i]))
218  m_histos[ci][i] = 0;
219  else
220  m_histos[ci][i] += layerCounts[i];
221  layerCounts[i] = 0; // reset bin
222  }
223  prevlayer = (*it)->layer;
224  }
225 
226  // get bin range
227  std::pair<int, int> minMax = range((*it)->x, (*it)->ymin, (*it)->ymax, ci);
228  int binmin = minMax.first;
229  int binmax = minMax.second;
230 
231  // check wether we are within the Hough space
232  if (binmin >= m_nbins) continue;
233  if (binmax < 0) continue;
234 
235  // adjust boundaries if needed
236  if (binmin < 0) binmin = 0;
237  if (binmax >= m_nbins) binmax = m_nbins - 1;
238 
239  // output hit for debug purposes
240  if (m_debug) {
241  std::cout << " cycle(theta layers) " << ci << " filling hit " << (*it)->x << " refpos "
242  << m_descriptor.referencePosition << " ymin " << (*it)->ymin << " ymax " << (*it)->ymax << " layer "
243  << (*it)->layer << " weight " << (*it)->w << " binmin " << binmin << " max " << binmax;
244  std::cout << " zmin " << binmin * m_binsize + m_descriptor.yMinRange << " zmax "
245  << binmax * m_binsize + m_descriptor.yMinRange;
246  if ((*it)->debugInfo()) {
247  const HitDebugInfo* db1 = (*it)->debugInfo();
248  std::cout << " sec " << db1->sector << " r " << db1->region << " type " << db1->type << " lay " << db1->layer
249  << " bc " << db1->barcode << std::endl;
250  } else
251  std::cout << std::endl;
252  }
253  int weight = sign * (*it)->w; // set the hit weight
254  // set bits to true
255  for (; binmin <= binmax; ++binmin) layerCounts[binmin] = weight;
256  } // end of loopin gover hits
257  // if the last set of hits was not filled, fill them now; just for the last layer
258  for (int i = 0; i < m_nbins; ++i) {
259  if (subtract && -layerCounts[i] >= static_cast<int>(m_histos[ci][i]))
260  m_histos[ci][i] = 0;
261  else
262  m_histos[ci][i] += layerCounts[i];
263  // if( m_debug && layerCounts[i] != 0 ) std::cout << " filling layer " << prevlayer << " bin " << i << " val " <<
264  // layerCounts[i] << " tot " << m_histos[ci][i] << std::endl;
265  layerCounts[i] = 0; // reset bin
266  }
267  }
268  }
269 
270  std::vector<TH1*> MuonLayerHough::rootHistos(const std::string& prefix, const float* rmi, const float* rma) const {
271  std::vector<TH1*> hists;
272  hists.reserve(m_histos.size());
273 
274  float rmin = rmi ? *rmi : m_descriptor.yMinRange;
275  float rmax = rma ? *rma : m_descriptor.yMaxRange;
276 
277  int cycles = m_histos.size();
278  for (int ci = 0; ci < cycles; ++ci) {
279  TString hname = prefix + "_hist";
280  hname += ci;
281  TH1F* h = new TH1F(hname, hname, m_nbins, rmin, rmax); // register the new histograms here
282  for (int n = 0; n < m_nbins; ++n) h->SetBinContent(n + 1, m_histos[ci][n] * 0.001);
283  hists.push_back(h);
284  }
285  return hists;
286  }
287 
289  const float preMaxCut = selector.getMinCutValue(); // in this case it is likely to be 1.9
290  maximum.max = 0;
291  maximum.pos = 0;
292  maximum.theta = 0;
293  maximum.refpos = 0;
294  maximum.refregion = -1;
295  maximum.refchIndex = -1;
296  maximum.binpos = -1;
297  maximum.binposmin = -1;
298  maximum.binposmax = -1;
299  maximum.binsize = -1;
300  maximum.bintheta = -1;
301  maximum.triggerConfirmed = 0;
302  maximum.hits.clear();
303  maximum.hough = this;
304 
305  if (preMaxCut < 0) return false; // this is where there is no cut
306 
307  float tmax = 0;
308  int posb = -1;
309  int thetab = -1;
310  int imaxval = preMaxCut * 1000;
311  const int cycles = m_histos.size();
312  // loop over histograms and find maximum
313  for (int ci = 0; ci < cycles; ++ci) {
314  const float scale =
315  1. - 0.01 * std::abs(ci - cycles / 2.0); // small deweighting of non pointing bins; weighting more on the central part
316  for (int n = 0; n < m_nbins; ++n) {
317  const int val = m_histos[ci][n];
318  if (val < imaxval) continue;
319 
320  if (scale * val > tmax) {
321  tmax = scale * val;
322  posb = n;
323  thetab = ci;
324  }
325  }
326  }
327  if (posb == -1) return false; // didn't find a max
328 
329  const float candidatePos = m_descriptor.yMinRange + m_binsize * posb;
330  if (tmax < selector.getCutValue(candidatePos)) return false;
331 
332  maximum.max = tmax * 0.001;
333  maximum.pos = candidatePos;
334  maximum.theta =
335  -m_descriptor.thetaStep * (thetab - (m_histos.size() - 1) * 0.5) + std::atan2(m_descriptor.referencePosition, maximum.pos);
337  maximum.refregion = m_descriptor.region;
338  maximum.refchIndex = m_descriptor.chIndex;
339  maximum.binpos = posb;
340  maximum.binposmin = posb;
341  maximum.binposmax = posb;
342  maximum.binsize = m_binsize;
343  maximum.bintheta = thetab;
344 
345  // determin width of maximum; now to be 70% of the original height
346  unsigned int imax = m_histos[thetab][posb];
347  unsigned int sidemax = 0.7 * imax;
348  // loop down, catch case the maximum is in the first bin
349  for (int n = posb != 0 ? posb - 1 : posb; n >= 0; --n) {
350  if (m_histos[thetab][n] > sidemax) {
351  maximum.binposmin = n;
352  } else {
353  break;
354  }
355  }
356  for (int n = posb + 1; n < m_nbins; ++n) {
357  if (m_histos[thetab][n] > sidemax) {
358  maximum.binposmax = n;
359  } else {
360  break;
361  }
362  }
363  return true;
364  }
365 
367  if (maximum.bintheta == -1 || maximum.binposmax == -1 || maximum.binposmin == -1) return;
368  // loop over hits and find those that are compatible with the maximum
369  HitVec::const_iterator it = hits.begin();
370  HitVec::const_iterator it_end = hits.end();
371  for (; it != it_end; ++it) {
372  // calculate the bins associated with the hit and check whether any of they are part of the maximum
373  std::pair<int, int> minMax = range((*it)->x, (*it)->ymin, (*it)->ymax, maximum.bintheta);
374  if (m_debug) {
375  std::cout << " associating hit: x " << (*it)->x << " ymin " << (*it)->ymin << " ymax " << (*it)->ymax << " layer "
376  << (*it)->layer << " range " << minMax.first << " " << minMax.second << " max range " << maximum.binposmin
377  << " " << maximum.binposmax << " prd " << (*it)->prd << " tgc " << (*it)->tgc;
378  if ((*it)->debugInfo()) std::cout << " trigC " << (*it)->debugInfo()->trigConfirm;
379  std::cout << std::endl;
380  }
381  if (minMax.first > maximum.binposmax) continue; // minimum bin large than the maximum, drop
382  if (minMax.second < maximum.binposmin) continue; // maximum bin smaller than the minimum, drop
383  // keep everything else
384  maximum.hits.push_back(*it);
385  if ((*it)->debugInfo() && (*it)->debugInfo()->trigConfirm > 0) ++maximum.triggerConfirmed;
386  }
387  }
388 
389  std::pair<float, float> MuonLayerHough::maximum(float x, float y, int& posbin, int& thetabin) const {
390  unsigned int max = 0;
391  thetabin = -1;
392  posbin = -1;
393 
394  int cycles = m_histos.size();
395  for (int ci = 0; ci < cycles; ++ci) {
396  int relbin = ci - (cycles - 1) / 2;
397  float dtheta = m_descriptor.thetaStep;
398  float dthetaOffset = 2 * m_descriptor.thetaStep * (relbin); // if bintheta = cycles, this is the same as dtheta * (cycles + 1 )
399  float theta = std::atan2(x, y);
400  float z0 = (m_descriptor.referencePosition - x) * cot(theta - dthetaOffset + dtheta) +
401  y; // move the angle by a step, recalculate the new y value
402  float z1 = (m_descriptor.referencePosition - x) * cot(theta - dthetaOffset - dtheta) + y;
403 
404  float zmin = z0 < z1 ? z0 : z1;
405  float zmax = z0 < z1 ? z1 : z0;
406  int binmin = (zmin - m_descriptor.yMinRange) / m_binsize - 1;
407  int binmax = (zmax - m_descriptor.yMinRange) / m_binsize + 1;
408  if (binmin >= m_nbins) continue;
409  if (binmax < 0) continue;
410 
411  if (binmin < 0) binmin = 0;
412  if (binmax >= m_nbins) binmax = m_nbins - 1;
413  for (int n = binmin; n < binmax; ++n) {
414  if (max < m_histos[ci][n]) {
415  max = m_histos[ci][n];
416  posbin = n;
417  thetabin = ci;
418  }
419  }
420  }
421  float dthetaOffset = 2 * m_descriptor.thetaStep * (thetabin - (cycles - 1) / 2.);
422  return std::make_pair(0.001 * max, -dthetaOffset);
423  }
424 
425  float MuonLayerHough::layerConfirmation(float x, float y, float range) const {
426  unsigned int max = 0;
427  if (x == 0) return 0;
428 
429  float yloc = m_descriptor.referencePosition * y / x;
430  int bincenter = (yloc - m_descriptor.yMinRange) / m_binsize;
431  int scanRange = range / m_binsize;
432  int binmin = bincenter - scanRange;
433  int binmax = bincenter + scanRange;
434  if (binmin >= m_nbins) return 0;
435  if (binmax < 0) return 0;
436 
437  int maxbin = -1;
438  if (binmin < 0) binmin = 0;
439  if (binmax >= m_nbins) binmax = m_nbins - 1;
440  int cyclesmin = 0;
441  int cycles = m_histos.size();
442  if (cycles > 3) {
443  int c0 = cycles / 2;
444  cyclesmin = c0 - 1;
445  cycles = c0 + 1;
446  }
447  for (int n = binmin; n < binmax; ++n) {
448  for (int ci = cyclesmin; ci < cycles; ++ci) {
449  if (max < m_histos[ci][n]) {
450  max = m_histos[ci][n];
451  maxbin = n;
452  }
453  }
454  }
455  if (range == 900) std::cout << " layerConfirmation " << max << " bin " << maxbin << " entry " << bincenter << std::endl;
456  return 0.001 * max;
457  }
458 
459  std::pair<float, float> MuonLayerHough::layerConfirmation(unsigned int thetaBin, float x, float y, float range) const {
460  unsigned int max = 0;
461  if (x == 0) return std::make_pair(0., 0.);
462  if (thetaBin >= m_histos.size()) return std::make_pair(0., 0.);
463  float yloc = m_descriptor.referencePosition * y / x;
464  int bincenter = (yloc - m_descriptor.yMinRange) / m_binsize;
465  int scanRange = range / m_binsize;
466  int binmin = bincenter - scanRange;
467  int binmax = bincenter + scanRange;
468  if (binmin >= m_nbins) return std::make_pair(0., 0.);
469  if (binmax < 0) return std::make_pair(0., 0.);
470 
471  int maxbin = -1;
472  if (binmin < 0) binmin = 0;
473  if (binmax >= m_nbins) binmax = m_nbins - 1;
474  for (int n = binmin; n < binmax; ++n) {
475  if (max < m_histos[thetaBin][n]) {
476  max = m_histos[thetaBin][n];
477  maxbin = n;
478  }
479  }
480  std::cout << " layerConfirmation " << max << " bin " << maxbin << " entry " << bincenter << " val " << yval(max) << std::endl;
481  return std::make_pair(0.001 * max, yval(maxbin));
482  }
483 
484  std::pair<int, int> MuonLayerHough::range(const float x, const float y1, const float y2, const int bintheta) const {
485  int cycles = m_histos.size();
487  float dtheta = m_descriptor.thetaStep;
488  if (dtheta <= 0)
489  throw std::runtime_error(
490  Form("File: %s, Line: %d\nMuonLayerHough::range() - dtheta is not positive (%.4f)", __FILE__, __LINE__, dtheta));
491  float dthetaOffset = 2 * dtheta * (bintheta - (cycles - 1) / 2.);
492 
493  float theta1 = std::atan2(x, y1) - dthetaOffset;
494  float z01 = dx * cot(theta1 + dtheta) + y1;
495  float z11 = dx * cot(theta1 - dtheta) + y1;
496  float zmin1 = std::min(z01, z11);
497  float zmax1 = std::max(z01, z11);
498 
499  float theta2 = std::atan2(x, y2) - dthetaOffset;
500  float z02 = dx * cot(theta2 + dtheta) + y2;
501  float z12 = dx * cot(theta2 - dtheta) + y2;
502  float zmin2 = std::min(z02, z12);
503  float zmax2 = std::max(z02, z12);
504 
505  const float zmin = std::min(zmin1, zmin2);
506  const float zmax = std::max(zmax1, zmax2);
507 
510  constexpr float cavern_size = 100.*Gaudi::Units::meter;
511  const float flt_lower_bin = std::max(-cavern_size, (zmin - m_descriptor.yMinRange) * m_invbinsize);
512  const float flt_upper_bin = std::min(cavern_size, (zmax - m_descriptor.yMinRange) * m_invbinsize);
513  const int lower_bin = std::floor(flt_lower_bin);
514  const int upper_bin = std::floor(flt_upper_bin);
515 
516  return std::make_pair(lower_bin, upper_bin); // convert the output to bins
517  }
518 
519  float extrapolate(const MuonLayerHough::Maximum& ref, const MuonLayerHough::Maximum& ex, bool doparabolic) {
520  // z is always the precision plane. r is the reference plane, simple and robust
521  double ref_z = ref.getGlobalZ();
522  double ref_r = ref.getGlobalR();
523  double ex_z = ex.getGlobalZ();
524  double ex_r = ex.getGlobalR();
525  float theta_ref = ref.getGlobalTheta();
526  if (!doparabolic || ref_z == 0 || theta_ref == 0) { // do linear extrapolation
527  if (!ex.isEndcap()) { // if extrapolate to barell
528  return ex_z - ex_r / ref_r * ref_z;
529  } else { // if extrapolate to endcap
530  return ex_r - ex_z * ref_r / ref_z;
531  }
532  } else { // do parabolic
533  float expected = 0;
534  float extrapolated_diff = 9999;
535  const float tan_theta_ref = std::tan(theta_ref);
536  const float invtan_theta_ref = 1. / tan_theta_ref;
537  float r_start = ref.hough->m_descriptor.chIndex % 2 > 0
538  ? 4900.
539  : 5200.; // start of barrel B field; values could be further optimized; 5500.:6500.
540  float z_start = 8500.; // start of endcap B field; used to be 6500; should start at 8500
541  float z_end = 12500.; // end of endcap B field; used to be 12000; should end at 12500
542  float r_SL = ref_r + (ex_z - ref_z) * tan_theta_ref;
543  float z_SL = ref_z + (ex_r - ref_r) * invtan_theta_ref;
544  // start extrapolation
545  if (!ref.isEndcap()) { // this is starting with barrel chamber; BEE is 6
546  float rgeo = ref_r * ref_r - r_start * r_start;
547  float rhoInv = (invtan_theta_ref * ref_r - ref_z) / rgeo;
548  if (!ex.isEndcap()) { // this ends with barrel chamber; BEE is 6
549  expected = ref_z + (ex_r - ref_r) * invtan_theta_ref + (ex_r - ref_r) * (ex_r - ref_r) * rhoInv;
550  extrapolated_diff = ex_z - expected;
551  } // end with barrel to barrel extrapolation
552  else { // this is a barrel to endcap extrapolation, mostly in the transition region
553  // recalculate z start
554  z_start = ref_z + (r_start - ref_r) * invtan_theta_ref + (r_start - ref_r) * (r_start - ref_r) * rhoInv;
555  float zgeo = ref_z * ref_z - z_start * z_start;
556  float rho = (tan_theta_ref * ref_z - ref_r) / zgeo;
557  expected = ref_r + (ex_z - ref_z) * tan_theta_ref + (ex_z - ref_z) * (ex_z - ref_z) * rho;
558  extrapolated_diff = ex_r - expected;
559  }
560  } else { // this starts with endcap chamber;
561  // swap the starting position if on the other side
562  if (tan_theta_ref < 0) {
563  z_start = -z_start;
564  z_end = -z_end;
565  }
566  if (ex.isEndcap()) { // extrapolate to endcap
567  if (std::abs(ref_z) < std::abs(z_end)) { // extrapolate from EI or EE, have to use linear
568  expected = r_SL;
569  } else { // from EM or EO or EE
570  if (std::abs(ex_z) > std::abs(z_start)) { // extrapolate to EM or EO
571  // extrapolate to outer layer, just using theta of the middle measurement; only works if the theta measurement
572  // is correct can extrapolate with either the outside theta or middle theta; outside theta is better; farther
573  // away from the B field
574  expected = r_SL;
575  } else { // to EI
576  float r_end = ref_r + (z_end - ref_z) * tan_theta_ref;
577  float zgeo = z_start * z_start - z_end * z_end;
578  float rhoInv = (r_end - z_end * tan_theta_ref) / zgeo;
579  float tantheta = tan_theta_ref - 2 * (z_end - z_start) * rhoInv;
580  expected = ex_z * tantheta + (ex_z - z_start) * (ex_z - z_start) * rhoInv;
581  }
582  }
583  extrapolated_diff = ex_r - expected;
584  } else { // exrapolate to barrel; again verly likely to be in transition region, complicated B field; just use linear
585  expected = z_SL;
586  extrapolated_diff = ex_z - expected;
587  }
588  }
589  return extrapolated_diff;
590  } // end of parabolic extrapolation
591  }
592 } // namespace MuonHough
TRTCalib_Extractor.hits
hits
Definition: TRTCalib_Extractor.py:35
max
#define max(a, b)
Definition: cfImp.cxx:41
MuonHough::RegionDescriptor::region
Muon::MuonStationIndex::DetectorRegionIndex region
Definition: MuonLayerHough.h:41
MuonHough::HitVec
std::vector< std::shared_ptr< MuonHough::Hit > > HitVec
Definition: MuonLayerHough.h:21
MuonHough::MuonLayerHough::yval
float yval(int posBin) const
access to y coordinate of a given bin
Definition: MuonLayerHough.h:196
getMessageSvc.h
singleton-like access to IMessageSvc via open function and helper
MuonHough::MuonLayerHough::m_binsize
float m_binsize
Definition: MuonLayerHough.h:168
dqt_zlumi_pandas.hname
string hname
Definition: dqt_zlumi_pandas.py:279
PixelAthClusterMonAlgCfg.zmin
zmin
Definition: PixelAthClusterMonAlgCfg.py:169
MuonHough::MuonLayerHough::layerConfirmation
float layerConfirmation(const Hit &hit, float range=1000.) const
calculate the highest value of the hough transform within the specified range for the given hit
Definition: MuonLayerHough.h:198
MuonHough::MuonLayerHough::max
unsigned int max
Definition: MuonLayerHough.h:174
MuonHough::MuonLayerHough::m_descriptor
RegionDescriptor m_descriptor
Definition: MuonLayerHough.h:179
MuonHough::MuonLayerHough::m_nbins
int m_nbins
inverse binsize
Definition: MuonLayerHough.h:172
MuonHough::HitDebugInfo
struct containing additional debug information on the hits that is not needed for the actual alg but ...
Definition: MuonSpectrometer/MuonReconstruction/MuonRecUtils/MuonLayerHough/MuonLayerHough/Hit.h:26
MuonHough::MuonLayerHough::m_debug
bool m_debug
Definition: MuonLayerHough.h:177
MuonHough::RegionDescriptor::yMinRange
float yMinRange
Definition: MuonLayerHough.h:44
skel.it
it
Definition: skel.GENtoEVGEN.py:396
MuonHough::RegionDescriptor::thetaStep
float thetaStep
Definition: MuonLayerHough.h:47
MuonHough::RegionDescriptor::nthetaSamples
unsigned int nthetaSamples
Definition: MuonLayerHough.h:48
MuonHough::MuonLayerHoughSelector
Definition: MuonLayerHoughSelector.h:12
MuonHough::RegionDescriptor::yMaxRange
float yMaxRange
Definition: MuonLayerHough.h:45
MuonHough::RegionDescriptor::referencePosition
float referencePosition
Definition: MuonLayerHough.h:43
yodamerge_tmp.scale
scale
Definition: yodamerge_tmp.py:138
MuonHough::RegionDescriptor
struct containing all information to build a Hough transform for a given chamber index
Definition: MuonLayerHough.h:25
MuonHough::MuonLayerHough::Maximum::getGlobalR
float getGlobalR() const
Definition: MuonLayerHough.h:85
x
#define x
MuonHough::MuonLayerHough::maxhist
int maxhist
Definition: MuonLayerHough.h:175
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:189
MuonLayerHough.h
bsCompare.db1
int db1
Definition: bsCompare.py:40
hotSpotInTAG.c0
c0
Definition: hotSpotInTAG.py:192
MuonHough::MuonLayerHough::Maximum
struct representing the maximum in the hough space
Definition: MuonLayerHough.h:56
python.SystemOfUnits.meter
int meter
Definition: SystemOfUnits.py:61
MuonHough::MuonLayerHough::fill
void fill(const Hit &hit)
fill the hough space with a single position
Definition: MuonLayerHough.h:200
MuonHough::MuonLayerHough::range
std::pair< int, int > range(const float x, const float y1, const float y2, const int bintheta) const
calculates the first and last bin the hit should be filled in for a given theta bin
Definition: MuonLayerHough.cxx:484
lumiFormat.i
int i
Definition: lumiFormat.py:85
MuonHough::MuonLayerHough::maximum
std::pair< float, float > maximum(float x, float y, int &posbin, int &thetabin) const
returns a pair with the position and angle corresponing to the input x,y values
Definition: MuonLayerHough.cxx:389
beamspotman.n
n
Definition: beamspotman.py:731
MuonHough::MuonLayerHough::findMaximum
bool findMaximum(Maximum &maximum, const MuonLayerHoughSelector &selector) const
find the highest maximum that is above maxval
Definition: MuonLayerHough.cxx:288
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
MuonHough
Definition: MuonLayerHoughTool.h:41
MuonHough::RegionDescriptor::chIndex
Muon::MuonStationIndex::ChIndex chIndex
Definition: MuonLayerHough.h:42
plotBeamSpotVxVal.range
range
Definition: plotBeamSpotVxVal.py:195
checkCorrelInHIST.prefix
dictionary prefix
Definition: checkCorrelInHIST.py:391
PixelAthClusterMonAlgCfg.zmax
zmax
Definition: PixelAthClusterMonAlgCfg.py:169
sign
int sign(int a)
Definition: TRT_StrawNeighbourSvc.h:107
MuonHough::z_end
constexpr double z_end
z value whereafter no magnetic field / curvature
Definition: MuonHoughMathUtils.h:28
MakeTH3DFromTH2Ds.hists
hists
Definition: MakeTH3DFromTH2Ds.py:72
MuonHough::MuonLayerHough::Maximum::getGlobalZ
float getGlobalZ() const
Definition: MuonLayerHough.h:89
MuonHough::MuonLayerHough::associateHitsToMaximum
void associateHitsToMaximum(Maximum &maximum, const HitVec &hits) const
associates the list of input hits to the provided maximum
Definition: MuonLayerHough.cxx:366
MuonHough::MuonLayerHough::m_invbinsize
float m_invbinsize
binsize
Definition: MuonLayerHough.h:169
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
TRT::Track::z0
@ z0
Definition: InnerDetector/InDetCalibEvent/TRT_CalibData/TRT_CalibData/TrackInfo.h:63
MuonHough::MuonLayerHough::MuonLayerHough
MuonLayerHough(const RegionDescriptor &descriptor)
constructor
Definition: MuonLayerHough.cxx:35
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
min
#define min(a, b)
Definition: cfImp.cxx:40
create_dcsc_inputs_sqlite.arg
list arg
Definition: create_dcsc_inputs_sqlite.py:48
MuonHough::MuonLayerHough::maxbin
int maxbin
Definition: MuonLayerHough.h:176
MuonHough::RegionDescriptor::yBinSize
float yBinSize
Definition: MuonLayerHough.h:46
MuonHough::MuonLayerHough::reset
void reset()
reset the transform
Definition: MuonLayerHough.cxx:48
MuonHough::MuonLayerHough::fillLayer
void fillLayer(const HitVec &hits, bool substract=false)
fill the hough space with a vector of hits using the layer mode
Definition: MuonLayerHough.cxx:100
python.selector.AtlRunQuerySelectorLhcOlc.selector
selector
Definition: AtlRunQuerySelectorLhcOlc.py:611
y
#define y
h
ref
const boost::regex ref(r_ef)
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
MuonHough::MuonLayerHough::m_histos
std::vector< std::unique_ptr< unsigned int[]> > m_histos
Definition: MuonLayerHough.h:178
MuonHough::MuonLayerHough::rootHistos
std::vector< TH1 * > rootHistos(const std::string &prefix, const float *rmin=0, const float *rmax=0) const
returns a vector with all the histograms of the hough as TH1*
Definition: MuonLayerHough.cxx:270
makeTRTBarrelCans.dx
tuple dx
Definition: makeTRTBarrelCans.py:20
MuonHough::MuonLayerHough::fillLayer2
void fillLayer2(const HitVec &hits, bool subtract=false)
Definition: MuonLayerHough.cxx:199
get
T * get(TKey *tobj)
get a TObject* from a TKey* (why can't a TObject be a TKey?)
Definition: hcg.cxx:127
MuonHough::extrapolate
float extrapolate(const MuonLayerHough::Maximum &ref, const MuonLayerHough::Maximum &ex, bool doparabolic=false)
Definition: MuonLayerHough.cxx:519
python.TrigEgammaMonitorHelper.TH1F
def TH1F(name, title, nxbins, bins_par2, bins_par3=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:24
python.IoTestsLib.w
def w
Definition: IoTestsLib.py:200
MuonHough::MuonLayerHough::Maximum::isEndcap
bool isEndcap() const
Definition: MuonLayerHough.h:77
fitman.rho
rho
Definition: fitman.py:532