41 return "no such plane";
44 CscPlane findPlane(
int station,
bool measphi) {
50 }
else if (station == 2) {
63 std::vector<double> pfac;
65 if (qrat < 0.095)
return 1;
66 if (qrat > 1.01)
return 1;
68 pfac.push_back(-5.483);
69 pfac.push_back(130.6);
70 pfac.push_back(-1296.);
71 pfac.push_back(5994.);
72 pfac.push_back(-10580.);
73 }
else if (qrat < 0.30) {
74 pfac.push_back(-0.9225);
75 pfac.push_back(5.569);
76 pfac.push_back(2.908);
77 pfac.push_back(-66.19);
78 pfac.push_back(108.5);
80 pfac.push_back(-0.4246);
81 pfac.push_back(2.619);
82 pfac.push_back(-3.642);
83 pfac.push_back(2.837);
84 pfac.push_back(-0.8916);
87 if (qrat < 0.105)
return 1;
88 if (qrat > 1.01)
return 1;
90 pfac.push_back(-2.823);
91 pfac.push_back(42.71);
92 pfac.push_back(-279.1);
93 pfac.push_back(879.2);
94 pfac.push_back(-1062.);
96 pfac.push_back(-0.5409);
97 pfac.push_back(3.110);
98 pfac.push_back(-4.630);
99 pfac.push_back(3.800);
100 pfac.push_back(-1.241);
108 for (
unsigned int ipow = 0; ipow < pfac.size(); ++ipow) {
109 dcordqrat += ipow * pfac[ipow] * term;
110 if (ipow) term *= qrat;
111 cor += pfac[ipow] * term;
128 int qrat_interpolation(
double qrmin,
const std::vector<double>& corvals,
double qrat,
double& cor,
double& dcordqrat) {
129 int nbin = corvals.size();
132 if (qrat < qrmin) qrat = qrmin;
134 double dqrat = 1.0 / nbin;
135 int bin =
int(qrat / dqrat);
138 double x1 = corvals[
bin];
139 if (
x1 == 0.0)
return 2;
140 double x0 =
bin == 0 ? 0.0 : corvals[
bin - 1];
141 double x2 =
bin == nbin - 1 ? 1.0 : corvals[
bin + 1];
144 double qrat0 = dqrat *
bin;
145 double qrat1 = qrat0 + dqrat;
148 if (qrmin > qrat0 && qrmin < qrat1) { qrat0 = qrmin; }
154 if (
bin == nbin - 1) {
156 double b = (
x1 - x0) / dqrat;
157 cor =
a +
b * (qrat - qrat0);
160 double d0 = qrat1 - qrat0;
166 cor =
a +
b * (qrat - qrat1) +
c * (qrat - qrat1) * (qrat - qrat1);
167 dcordqrat =
b + 2.0 *
c * (qrat - qrat1);
169 if (cor < 0.0) cor = 0.0;
194 int qrat_atanh(
const double a,
const double b,
double c,
const double x0,
double qrat,
double& cor,
double& dcordqrat) {
195 if (qrat <= 0)
return 1;
197 double qrmin =
a +
b * tanh(
c * (-1 - x0));
200 if (qrat < qrmin) qrat = qrmin;
201 double z = (qrat -
a) /
b;
202 cor = atanh(
z) /
c + x0;
203 dcordqrat = 1.0 / (1.0 -
z *
z) /
b /
c;
212 declareInterface<ICscClusterFitter>(
this);
297 return StatusCode::SUCCESS;
305 bool dofixed =
false;
306 bool docharge =
false;
315 if (dofixed || docharge) {
316 dnames.emplace_back(
"scor1");
317 dnames.emplace_back(
"scor2");
318 dnames.emplace_back(
"scor");
321 dnames.emplace_back(
"dscor1");
322 dnames.emplace_back(
"dscor2");
323 dnames.emplace_back(
"dscor");
324 dnames.emplace_back(
"scordiff");
325 dnames.emplace_back(
"dscordiff");
341 unsigned int nstrip = sfits.size();
345 Muon::CscTimeStatus tstatus = (sfits[0].charge > sfits[1].charge) ? sfits[0].timeStatus : sfits[1].timeStatus;
347 }
else if (nstrip == 1) {
355 for (
unsigned int istrip = 0; istrip < nstrip; ++istrip) {
356 if (sfits[istrip].
strip ==
nullptr) {
370 if (MuonDetMgr ==
nullptr) {
371 ATH_MSG_ERROR(
"Null pointer to the MuonDetectorManager conditions object");
376 bool measphi =
m_idHelperSvc->cscIdHelper().CscIdHelper::measuresPhi(idStrip0);
379 unsigned int strip0 =
m_idHelperSvc->cscIdHelper().strip(idStrip0) - 1;
380 int station =
m_idHelperSvc->cscIdHelper().stationName(idStrip0) - 49;
382 CscPlane plane = findPlane(station, measphi);
384 ATH_MSG_WARNING(
"Invalid CSC plane: station=" << station <<
"; measphi=" << measphi);
391 for (
unsigned int istrip = 0; istrip < nstrip; ++istrip) {
392 Identifier id = sfits[istrip].strip->identify();
394 <<
m_idHelperSvc->cscIdHelper().strip(
id) <<
" " << sfits[istrip].charge);
398 unsigned int istrip_peak = 0;
400 for (
unsigned int istrip = 1; istrip < nstrip - 1; ++istrip) {
404 float qthis = sfit.
charge;
405 float qlast = sfits[istrip - 1].charge;
406 float qnext = sfits[istrip + 1].charge;
408 bool ispeak = qthis > qlast && qthis > qnext;
413 if (qthis == qnext) { ispeak = (qthis > qlast) && (istrip + 2 == nstrip || sfits[istrip + 2].
charge < qthis); }
419 if (qthis == qlast) {
420 ispeak = qthis > qnext;
430 istrip_peak = istrip;
433 ATH_MSG_VERBOSE(
" Peak is at " << istrip_peak <<
"th strip in cluster");
436 if (strip0 <= 0 || strip0 + nstrip > maxstrip) {
452 bool is_even = 2 * (nstrip / 2) == nstrip;
453 bool atcenter = istrip_peak == nstrip / 2 || (is_even && istrip_peak + 1 == nstrip / 2);
480 double savg = istrip_peak;
490 double q1 = sfits[istrip_peak - 1].charge;
491 double q0 = sfits[istrip_peak].charge;
492 double q2 = sfits[istrip_peak + 1].charge;
493 double qrat1 = q1 / q0;
494 double qrat2 = q2 / q0;
495 double dq1 = sfits[istrip_peak - 1].dcharge;
496 double dq0 = sfits[istrip_peak].dcharge;
497 double dq2 = sfits[istrip_peak + 1].dcharge;
501 double dscordqrat1 = 0.;
503 double dscordqrat2 = 0.;
506 if (posopt ==
"POLYNOMIAL") {
509 }
else if (posopt ==
"TABLE") {
511 const std::vector<double>* pcor =
nullptr;
529 }
else if (posopt ==
"ATANH") {
554 if (stat1 || stat2) {
560 ATH_MSG_VERBOSE(
" QRAT derivatives: " << dscordqrat1 <<
" " << dscordqrat2);
565 dscordqrat1 = -dscordqrat1;
570 if (erropt ==
"FIXED") {
571 double w1 = 1 / (dscordqrat1 * dscordqrat1);
572 double w2 = 1 / (dscordqrat2 * dscordqrat2);
573 scor = (w1 * scor1 + w2 * scor2) / (w1 + w2);
574 double scor_diff = std::abs(scor2 - scor1);
578 dmap[
"scor1"] = scor1;
579 dmap[
"scor2"] = scor2;
588 }
else if (erropt ==
"CHARGE") {
590 double x1 = dscordqrat1 * qrat1;
591 double x2 = dscordqrat2 * qrat2;
592 double dqq0 = dq0 / q0;
593 double dqq1 = dq1 / q1;
594 double dqq2 = dq2 / q2;
595 double rnum = (-
x1 *
x1 * dqq1 * dqq1 +
x2 *
x2 * dqq2 * dqq2 + (
x2 *
x2 -
x1 *
x1) * dqq0 * dqq0);
596 double rden = (
x1 *
x1 * dqq1 * dqq1 +
x2 *
x2 * dqq2 * dqq2 + (
x2 -
x1) * (
x2 -
x1) * dqq0 * dqq0);
597 double rfac = rnum / rden;
599 double dscor_diff = sqrt(rden);
600 double scor_diff = scor2 - scor1;
601 double scor_sig = std::abs(scor_diff) / dscor_diff;
603 double w1 = 0.5 * (1.0 +
rfac);
604 double w2 = 0.5 * (1.0 -
rfac);
605 scor = w1 * scor1 + w2 * scor2;
607 double ddscor1 = w1 *
x1 * dqq1;
608 double ddscor2 = w2 *
x2 * dqq2;
609 double ddscor0 = (w1 *
x1 + w2 *
x2) * dqq0;
610 double dscor = sqrt(ddscor1 * ddscor1 + ddscor2 * ddscor2 + ddscor0 * ddscor0);
611 dpos = pitch * dscor;
615 double dscor1 = std::abs(
x1) * sqrt(dqq1 * dqq1 + dqq0 * dqq0);
616 double dscor2 = std::abs(
x2) * sqrt(dqq2 * dqq2 + dqq0 * dqq0);
617 dmap[
"scor1"] = scor1;
618 dmap[
"dscor1"] = dscor1;
619 dmap[
"scor2"] = scor2;
620 dmap[
"dscor2"] = dscor2;
621 dmap[
"scordiff"] = scor_diff;
622 dmap[
"dscordiff"] = dscor_diff;
624 dmap[
"dscor"] = dscor;
628 ATH_MSG_VERBOSE(
" dq1, dq0, dq2: " << dq1 <<
" " << dq0 <<
" " << dq2);
650 savg += scor + strip0;
651 ATH_MSG_VERBOSE(
" QRAT corr " << splane(plane) <<
" nstrip=" << nstrip <<
" savg=" << savg <<
" qrat1=" << qrat1
652 <<
" qrat2=" << qrat2 <<
" scor1=" << scor1 <<
" scor2=" << scor2);
662 res.strip = istrip_peak;
663 res.position = pitch * (savg + 0.5 - 0.5 * maxstrip);
670 res.dposition = sqrt(dpos * dpos + dpostht * dpostht);
673 res.lstrip = nstrip - 1;
674 res.time = sfits[istrip_peak].time;
675 res.time_beforeT0Corr = sfits[istrip_peak].time_beforeT0Corr;
676 res.time_beforeBPCorr = sfits[istrip_peak].time_beforeBPCorr;
677 res.timeStatus = sfits[istrip_peak].timeStatus;
688 res.charge_beforeBPCorr =
689 sfits[istrip_peak].charge_beforeBPCorr + sfits[istrip_peak - 1].charge_beforeBPCorr + sfits[istrip_peak + 1].charge_beforeBPCorr;
693 ATH_MSG_VERBOSE(
" Position :: pos=" <<
res.position <<
" dpos:dtht=" << dpos <<
":" << dpostht <<
" ==>" <<
res.dposition
694 <<
" at tanth=" << tantheta);
713 int station =
m_idHelperSvc->cscIdHelper().stationName(idStrip0) - 49;
726 double newError = sqrt(dpos * dpos - old_dpostht * old_dpostht + new_dpostht * new_dpostht);
728 ATH_MSG_VERBOSE(
" Position :: pos=" <<
pos <<
" dpos:newdpos=" << dpos <<
" : " << newError <<
" " << old_dpostht <<
" "
731 if (slope < -990) { newError = sqrt(dpos * dpos + old_dpostht * old_dpostht); }
741 for (
unsigned int iresult = 0; iresult <
results.size(); ++iresult) {
744 new_results.push_back(
res);
750 int station =
m_idHelperSvc->cscIdHelper().stationName(idStrip0) - 49;
753 double pos =
res.position;
761 double dpos =
res.dposition;
763 res.dposition = sqrt(dpos * dpos + dpostht * dpostht);
766 new_results.push_back(
res);