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()) {
129 std::cout <<
" sec " <<
db1->sector <<
" r " <<
db1->region <<
" type " <<
db1->type <<
" lay " <<
db1->layer
130 <<
" slay " <<
db1->sublayer << std::endl;
132 std::cout << std::endl;
135 if (prevbinmax == -1) {
136 if (
m_debug) std::cout <<
" first range " << binmin <<
" " << binmax << std::endl;
139 prevlayer = (*it)->layer;
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;
148 if (prevbinmax < binmin || prevlayer != (*it)->layer) {
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) {
154 int w = 1000 * (*it)->w;
155 if (subtract)
w *= -1;
156 if (
w < 0 && (
int)
val < -
w)
168 prevlayer = (*it)->layer;
173 std::cout <<
" updating range " << prevbinmin <<
" " << prevbinmax <<
" hit " << binmin <<
" " << binmax
174 <<
" new " << prevbinmin <<
" " << binmax << std::endl;
178 if (prevbinmax != -1) {
180 std::cout <<
" filling last range " << prevbinmin <<
" " << prevbinmax <<
" weight " <<
hits.back()->w << std::endl;
181 for (
int n = prevbinmin;
n <= prevbinmax; ++
n) {
183 int w = 1000 *
hits.back()->w;
184 if (subtract)
w *= -1;
185 if (
w < 0 && (
int)
val < -
w)
200 if (
hits.empty())
return;
202 std::vector<int> layerCounts(
m_nbins, 0);
203 int sign = subtract ? -1000 : 1000;
206 for (
int ci = 0; ci < cycles; ++ci) {
208 int prevlayer =
hits.front()->layer;
211 HitVec::const_iterator
it =
hits.begin();
212 HitVec::const_iterator it_end =
hits.end();
213 for (;
it != it_end; ++
it) {
215 if (prevlayer != (*it)->layer) {
217 if (subtract && -layerCounts[
i] >=
static_cast<int>(
m_histos[ci][
i]))
223 prevlayer = (*it)->layer;
227 std::pair<int, int> minMax =
range((*it)->x, (*it)->ymin, (*it)->ymax, ci);
228 int binmin = minMax.first;
229 int binmax = minMax.second;
232 if (binmin >=
m_nbins)
continue;
233 if (binmax < 0)
continue;
236 if (binmin < 0) binmin = 0;
241 std::cout <<
" cycle(theta layers) " << ci <<
" filling hit " << (*it)->x <<
" refpos "
243 << (*it)->layer <<
" weight " << (*it)->w <<
" binmin " << binmin <<
" max " << binmax;
246 if ((*it)->debugInfo()) {
248 std::cout <<
" sec " <<
db1->sector <<
" r " <<
db1->region <<
" type " <<
db1->type <<
" lay " <<
db1->layer
249 <<
" bc " <<
db1->barcode << std::endl;
251 std::cout << std::endl;
255 for (; binmin <= binmax; ++binmin) layerCounts[binmin] =
weight;
259 if (subtract && -layerCounts[
i] >=
static_cast<int>(
m_histos[ci][
i]))
271 std::vector<TH1*>
hists;
278 for (
int ci = 0; ci < cycles; ++ci) {
289 const float preMaxCut =
selector.getMinCutValue();
305 if (preMaxCut < 0)
return false;
310 int imaxval = preMaxCut * 1000;
313 for (
int ci = 0; ci < cycles; ++ci) {
315 1. - 0.01 * std::abs(ci - cycles / 2.0);
318 if (
val < imaxval)
continue;
327 if (posb == -1)
return false;
330 if (tmax <
selector.getCutValue(candidatePos))
return false;
347 unsigned int sidemax = 0.7 *
imax;
349 for (
int n = posb != 0 ? posb - 1 : posb;
n >= 0; --
n) {
369 HitVec::const_iterator
it =
hits.begin();
370 HitVec::const_iterator it_end =
hits.end();
371 for (;
it != it_end; ++
it) {
373 std::pair<int, int> minMax =
range((*it)->x, (*it)->ymin, (*it)->ymax,
maximum.bintheta);
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;
381 if (minMax.first >
maximum.binposmax)
continue;
382 if (minMax.second <
maximum.binposmin)
continue;
385 if ((*it)->debugInfo() && (*it)->debugInfo()->trigConfirm > 0) ++
maximum.triggerConfirmed;
390 unsigned int max = 0;
395 for (
int ci = 0; ci < cycles; ++ci) {
396 int relbin = ci - (cycles - 1) / 2;
399 float theta = std::atan2(
x,
y);
408 if (binmin >=
m_nbins)
continue;
409 if (binmax < 0)
continue;
411 if (binmin < 0) binmin = 0;
413 for (
int n = binmin;
n < binmax; ++
n) {
422 return std::make_pair(0.001 *
max, -dthetaOffset);
426 unsigned int max = 0;
427 if (
x == 0)
return 0;
432 int binmin = bincenter - scanRange;
433 int binmax = bincenter + scanRange;
434 if (binmin >=
m_nbins)
return 0;
435 if (binmax < 0)
return 0;
438 if (binmin < 0) binmin = 0;
447 for (
int n = binmin;
n < binmax; ++
n) {
448 for (
int ci = cyclesmin; ci < cycles; ++ci) {
455 if (
range == 900) std::cout <<
" layerConfirmation " <<
max <<
" bin " <<
maxbin <<
" entry " << bincenter << std::endl;
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.);
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.);
472 if (binmin < 0) binmin = 0;
474 for (
int n = binmin;
n < binmax; ++
n) {
480 std::cout <<
" layerConfirmation " <<
max <<
" bin " <<
maxbin <<
" entry " << bincenter <<
" val " <<
yval(
max) << std::endl;
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.);
493 float theta1 = std::atan2(
x,
y1) - dthetaOffset;
494 float z01 =
dx * cot(theta1 + dtheta) +
y1;
495 float z11 =
dx * cot(theta1 - dtheta) +
y1;
499 float theta2 = std::atan2(
x,
y2) - dthetaOffset;
500 float z02 =
dx * cot(theta2 + dtheta) +
y2;
501 float z12 =
dx * cot(theta2 - dtheta) +
y2;
513 const int lower_bin = std::floor(flt_lower_bin);
514 const int upper_bin = std::floor(flt_upper_bin);
516 return std::make_pair(lower_bin, upper_bin);
521 double ref_z =
ref.getGlobalZ();
522 double ref_r =
ref.getGlobalR();
525 float theta_ref =
ref.getGlobalTheta();
526 if (!doparabolic || ref_z == 0 || theta_ref == 0) {
528 return ex_z - ex_r / ref_r * ref_z;
530 return ex_r - ex_z * ref_r / ref_z;
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
540 float z_start = 8500.;
541 float z_end = 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;
545 if (!
ref.isEndcap()) {
546 float rgeo = ref_r * ref_r - r_start * r_start;
547 float rhoInv = (invtan_theta_ref * ref_r - ref_z) / rgeo;
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;
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;
562 if (tan_theta_ref < 0) {
567 if (std::abs(ref_z) < std::abs(
z_end)) {
570 if (std::abs(ex_z) > std::abs(z_start)) {
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;
583 extrapolated_diff = ex_r - expected;
586 extrapolated_diff = ex_z - expected;
589 return extrapolated_diff;