17 #include "GaudiKernel/MsgStream.h"
18 #include "GaudiKernel/SystemOfUnits.h"
26 float cot(
const float x) {
27 const float arg = M_PI_2 -
x;
29 if (std::abs(
arg - M_PI_2) <= FLT_EPSILON || std::abs(
arg + M_PI_2) <= FLT_EPSILON)
return 1.e8;
36 m_descriptor(descriptor) {
57 for (
int ci = 0; ci < cycles; ++ci) {
60 float theta = std::atan2(
x,
y);
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;
78 if (binmin < 0) binmin = 0;
81 std::cout <<
" filling " <<
x <<
" y " << binmin <<
" " << binmax <<
" w " <<
weight <<
" center " << bincenter <<
" range "
84 for (
int n = binmin;
n <= binmax; ++
n) {
87 if (
w < 0 && (
int)
val < -
w)
101 if (
hits.empty())
return;
105 for (
int ci = 0; ci < cycles; ++ci) {
110 int prevbinmin = 10000;
113 HitVec::const_iterator
it =
hits.begin();
114 HitVec::const_iterator it_end =
hits.end();
115 for (;
it != it_end; ++
it) {
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);
122 if (binmin >=
m_nbins || binmax < 0)
continue;
126 <<
" layer " << (*it)->layer <<
" binmin " << binmin <<
" max " << binmax;
127 if ((*it)->debugInfo()) {
131 <<
" slay " <<
db1->sublayer << std::endl;
133 std::cout << std::endl;
136 if (prevbinmax == -1) {
137 if (
m_debug) std::cout <<
" first range " << binmin <<
" " << binmax << std::endl;
140 prevlayer = (*it)->layer;
144 if (binmin < prevbinmin && prevlayer == (*it)->layer) {
145 std::cout <<
"Error hits are out of order: min " << binmin <<
" max " << binmax <<
" lay " << (*it)->layer << std::endl;
149 if (prevbinmax < binmin || prevlayer != (*it)->layer) {
151 std::cout <<
" filling range " << prevbinmin <<
" " << prevbinmax <<
" new min " << binmin <<
" " << binmax
152 <<
" weight " << (*it)->w << std::endl;
153 for (
int n = prevbinmin;
n <= prevbinmax; ++
n) {
155 int w = 1000 * (*it)->w;
156 if (subtract)
w *= -1;
157 if (
w < 0 && (
int)
val < -
w)
169 prevlayer = (*it)->layer;
174 std::cout <<
" updating range " << prevbinmin <<
" " << prevbinmax <<
" hit " << binmin <<
" " << binmax
175 <<
" new " << prevbinmin <<
" " << binmax << std::endl;
179 if (prevbinmax != -1) {
181 std::cout <<
" filling last range " << prevbinmin <<
" " << prevbinmax <<
" weight " <<
hits.back()->w << std::endl;
182 for (
int n = prevbinmin;
n <= prevbinmax; ++
n) {
184 int w = 1000 *
hits.back()->w;
185 if (subtract)
w *= -1;
186 if (
w < 0 && (
int)
val < -
w)
201 if (
hits.empty())
return;
203 std::vector<int> layerCounts(
m_nbins, 0);
204 int sign = subtract ? -1000 : 1000;
207 for (
int ci = 0; ci < cycles; ++ci) {
209 int prevlayer =
hits.front()->layer;
212 HitVec::const_iterator
it =
hits.begin();
213 HitVec::const_iterator it_end =
hits.end();
214 for (;
it != it_end; ++
it) {
216 if (prevlayer != (*it)->layer) {
218 if (subtract && -layerCounts[
i] >=
static_cast<int>(
m_histos[ci][
i]))
224 prevlayer = (*it)->layer;
228 std::pair<int, int> minMax =
range((*it)->x, (*it)->ymin, (*it)->ymax, ci);
229 int binmin = minMax.first;
230 int binmax = minMax.second;
233 if (binmin >=
m_nbins)
continue;
234 if (binmax < 0)
continue;
237 if (binmin < 0) binmin = 0;
242 std::cout <<
" cycle(theta layers) " << ci <<
" filling hit " << (*it)->x <<
" refpos "
244 << (*it)->layer <<
" weight " << (*it)->w <<
" binmin " << binmin <<
" max " << binmax;
247 if ((*it)->debugInfo()) {
251 <<
" bc " <<
db1->barcode << std::endl;
253 std::cout << std::endl;
257 for (; binmin <= binmax; ++binmin) layerCounts[binmin] =
weight;
261 if (subtract && -layerCounts[
i] >=
static_cast<int>(
m_histos[ci][
i]))
273 std::vector<TH1*>
hists;
280 for (
int ci = 0; ci < cycles; ++ci) {
291 const float preMaxCut =
selector.getMinCutValue();
296 maximum.refregion = DetRegIdx::DetectorRegionUnknown;
297 maximum.refchIndex = ChIdx::ChUnknown;
307 if (preMaxCut < 0)
return false;
312 int imaxval = preMaxCut * 1000;
315 for (
int ci = 0; ci < cycles; ++ci) {
317 1. - 0.01 * std::abs(ci - cycles / 2.0);
320 if (
val < imaxval)
continue;
329 if (posb == -1)
return false;
332 if (tmax <
selector.getCutValue(candidatePos))
return false;
349 unsigned int sidemax = 0.7 *
imax;
351 for (
int n = posb != 0 ? posb - 1 : posb;
n >= 0; --
n) {
371 HitVec::const_iterator
it =
hits.begin();
372 HitVec::const_iterator it_end =
hits.end();
373 for (;
it != it_end; ++
it) {
375 std::pair<int, int> minMax =
range((*it)->x, (*it)->ymin, (*it)->ymax,
maximum.bintheta);
377 std::cout <<
" associating hit: x " << (*it)->x <<
" ymin " << (*it)->ymin <<
" ymax " << (*it)->ymax <<
" layer "
378 << (*it)->layer <<
" range " << minMax.first <<
" " << minMax.second <<
" max range " <<
maximum.binposmin
379 <<
" " <<
maximum.binposmax <<
" prd " << (*it)->prd <<
" tgc " << (*it)->tgc;
380 if ((*it)->debugInfo()) std::cout <<
" trigC " << (*it)->debugInfo()->trigConfirm;
381 std::cout << std::endl;
383 if (minMax.first >
maximum.binposmax)
continue;
384 if (minMax.second <
maximum.binposmin)
continue;
387 if ((*it)->debugInfo() && (*it)->debugInfo()->trigConfirm > 0) ++
maximum.triggerConfirmed;
392 unsigned int max = 0;
397 for (
int ci = 0; ci < cycles; ++ci) {
398 int relbin = ci - (cycles - 1) / 2;
401 float theta = std::atan2(
x,
y);
410 if (binmin >=
m_nbins)
continue;
411 if (binmax < 0)
continue;
413 if (binmin < 0) binmin = 0;
415 for (
int n = binmin;
n < binmax; ++
n) {
424 return std::make_pair(0.001 *
max, -dthetaOffset);
428 unsigned int max = 0;
429 if (
x == 0)
return 0;
434 int binmin = bincenter - scanRange;
435 int binmax = bincenter + scanRange;
436 if (binmin >=
m_nbins)
return 0;
437 if (binmax < 0)
return 0;
440 if (binmin < 0) binmin = 0;
449 for (
int n = binmin;
n < binmax; ++
n) {
450 for (
int ci = cyclesmin; ci < cycles; ++ci) {
457 if (
range == 900) std::cout <<
" layerConfirmation " <<
max <<
" bin " <<
maxbin <<
" entry " << bincenter << std::endl;
462 unsigned int max = 0;
463 if (
x == 0)
return std::make_pair(0., 0.);
464 if (thetaBin >=
m_histos.size())
return std::make_pair(0., 0.);
468 int binmin = bincenter - scanRange;
469 int binmax = bincenter + scanRange;
470 if (binmin >=
m_nbins)
return std::make_pair(0., 0.);
471 if (binmax < 0)
return std::make_pair(0., 0.);
474 if (binmin < 0) binmin = 0;
476 for (
int n = binmin;
n < binmax; ++
n) {
482 std::cout <<
" layerConfirmation " <<
max <<
" bin " <<
maxbin <<
" entry " << bincenter <<
" val " <<
yval(
max) << std::endl;
491 throw std::runtime_error(
492 Form(
"File: %s, Line: %d\nMuonLayerHough::range() - dtheta is not positive (%.4f)", __FILE__, __LINE__, dtheta));
493 float dthetaOffset = 2 * dtheta * (bintheta - (cycles - 1) / 2.);
495 float theta1 = std::atan2(
x,
y1) - dthetaOffset;
496 float z01 =
dx * cot(theta1 + dtheta) +
y1;
497 float z11 =
dx * cot(theta1 - dtheta) +
y1;
501 float theta2 = std::atan2(
x,
y2) - dthetaOffset;
502 float z02 =
dx * cot(theta2 + dtheta) +
y2;
503 float z12 =
dx * cot(theta2 - dtheta) +
y2;
515 const int lower_bin = std::floor(flt_lower_bin);
516 const int upper_bin = std::floor(flt_upper_bin);
518 return std::make_pair(lower_bin, upper_bin);
523 double ref_z =
ref.getGlobalZ();
524 double ref_r =
ref.getGlobalR();
527 float theta_ref =
ref.getGlobalTheta();
528 if (!doparabolic || ref_z == 0 || theta_ref == 0) {
530 return ex_z - ex_r / ref_r * ref_z;
532 return ex_r - ex_z * ref_r / ref_z;
536 float extrapolated_diff = 9999;
537 const float tan_theta_ref =
std::tan(theta_ref);
538 const float invtan_theta_ref = 1. / tan_theta_ref;
539 float r_start =
static_cast<int>(
ref.hough->m_descriptor.chIndex) % 2 > 0
542 float z_start = 8500.;
543 float z_end = 12500.;
544 float r_SL = ref_r + (ex_z - ref_z) * tan_theta_ref;
545 float z_SL = ref_z + (ex_r - ref_r) * invtan_theta_ref;
547 if (!
ref.isEndcap()) {
548 float rgeo = ref_r * ref_r - r_start * r_start;
549 float rhoInv = (invtan_theta_ref * ref_r - ref_z) / rgeo;
551 expected = ref_z + (ex_r - ref_r) * invtan_theta_ref + (ex_r - ref_r) * (ex_r - ref_r) * rhoInv;
552 extrapolated_diff = ex_z - expected;
556 z_start = ref_z + (r_start - ref_r) * invtan_theta_ref + (r_start - ref_r) * (r_start - ref_r) * rhoInv;
557 float zgeo = ref_z * ref_z - z_start * z_start;
558 float rho = (tan_theta_ref * ref_z - ref_r) / zgeo;
559 expected = ref_r + (ex_z - ref_z) * tan_theta_ref + (ex_z - ref_z) * (ex_z - ref_z) *
rho;
560 extrapolated_diff = ex_r - expected;
564 if (tan_theta_ref < 0) {
569 if (std::abs(ref_z) < std::abs(
z_end)) {
572 if (std::abs(ex_z) > std::abs(z_start)) {
578 float r_end = ref_r + (
z_end - ref_z) * tan_theta_ref;
579 float zgeo = z_start * z_start -
z_end *
z_end;
580 float rhoInv = (r_end -
z_end * tan_theta_ref) / zgeo;
581 float tantheta = tan_theta_ref - 2 * (
z_end - z_start) * rhoInv;
582 expected = ex_z * tantheta + (ex_z - z_start) * (ex_z - z_start) * rhoInv;
585 extrapolated_diff = ex_r - expected;
588 extrapolated_diff = ex_z - expected;
591 return extrapolated_diff;