30 #include "gsl/gsl_statistics.h"
39 debug =
pset.getParameter<
bool>(
"DebugMuonSeed");
40 sysError =
pset.getParameter<
double>(
"SeedPtSystematics");
42 DT12 =
pset.getParameter<std::vector<double> >(
"DT_12");
43 DT13 =
pset.getParameter<std::vector<double> >(
"DT_13");
44 DT14 =
pset.getParameter<std::vector<double> >(
"DT_14");
45 DT23 =
pset.getParameter<std::vector<double> >(
"DT_23");
46 DT24 =
pset.getParameter<std::vector<double> >(
"DT_24");
47 DT34 =
pset.getParameter<std::vector<double> >(
"DT_34");
49 CSC01 =
pset.getParameter<std::vector<double> >(
"CSC_01");
50 CSC12 =
pset.getParameter<std::vector<double> >(
"CSC_12");
51 CSC02 =
pset.getParameter<std::vector<double> >(
"CSC_02");
52 CSC13 =
pset.getParameter<std::vector<double> >(
"CSC_13");
53 CSC03 =
pset.getParameter<std::vector<double> >(
"CSC_03");
54 CSC14 =
pset.getParameter<std::vector<double> >(
"CSC_14");
55 CSC23 =
pset.getParameter<std::vector<double> >(
"CSC_23");
56 CSC24 =
pset.getParameter<std::vector<double> >(
"CSC_24");
57 CSC34 =
pset.getParameter<std::vector<double> >(
"CSC_34");
59 OL1213 =
pset.getParameter<std::vector<double> >(
"OL_1213");
60 OL1222 =
pset.getParameter<std::vector<double> >(
"OL_1222");
61 OL1232 =
pset.getParameter<std::vector<double> >(
"OL_1232");
62 OL1213 =
pset.getParameter<std::vector<double> >(
"OL_1213");
63 OL2222 =
pset.getParameter<std::vector<double> >(
"OL_1222");
65 SME11 =
pset.getParameter<std::vector<double> >(
"SME_11");
66 SME12 =
pset.getParameter<std::vector<double> >(
"SME_12");
67 SME13 =
pset.getParameter<std::vector<double> >(
"SME_13");
68 SME21 =
pset.getParameter<std::vector<double> >(
"SME_21");
69 SME22 =
pset.getParameter<std::vector<double> >(
"SME_22");
70 SME31 =
pset.getParameter<std::vector<double> >(
"SME_31");
71 SME32 =
pset.getParameter<std::vector<double> >(
"SME_32");
72 SME41 =
pset.getParameter<std::vector<double> >(
"SME_41");
74 SMB10 =
pset.getParameter<std::vector<double> >(
"SMB_10");
75 SMB11 =
pset.getParameter<std::vector<double> >(
"SMB_11");
76 SMB12 =
pset.getParameter<std::vector<double> >(
"SMB_12");
77 SMB20 =
pset.getParameter<std::vector<double> >(
"SMB_20");
78 SMB21 =
pset.getParameter<std::vector<double> >(
"SMB_21");
79 SMB22 =
pset.getParameter<std::vector<double> >(
"SMB_22");
80 SMB30 =
pset.getParameter<std::vector<double> >(
"SMB_30");
81 SMB31 =
pset.getParameter<std::vector<double> >(
"SMB_31");
82 SMB32 =
pset.getParameter<std::vector<double> >(
"SMB_32");
85 CSC01_1 =
pset.getParameter<std::vector<double> >(
"CSC_01_1_scale");
86 CSC12_1 =
pset.getParameter<std::vector<double> >(
"CSC_12_1_scale");
87 CSC12_2 =
pset.getParameter<std::vector<double> >(
"CSC_12_2_scale");
88 CSC12_3 =
pset.getParameter<std::vector<double> >(
"CSC_12_3_scale");
89 CSC13_2 =
pset.getParameter<std::vector<double> >(
"CSC_13_2_scale");
90 CSC13_3 =
pset.getParameter<std::vector<double> >(
"CSC_13_3_scale");
91 CSC14_3 =
pset.getParameter<std::vector<double> >(
"CSC_14_3_scale");
92 CSC23_1 =
pset.getParameter<std::vector<double> >(
"CSC_23_1_scale");
93 CSC23_2 =
pset.getParameter<std::vector<double> >(
"CSC_23_2_scale");
94 CSC24_1 =
pset.getParameter<std::vector<double> >(
"CSC_24_1_scale");
95 CSC34_1 =
pset.getParameter<std::vector<double> >(
"CSC_34_1_scale");
97 DT12_1 =
pset.getParameter<std::vector<double> >(
"DT_12_1_scale");
98 DT12_2 =
pset.getParameter<std::vector<double> >(
"DT_12_2_scale");
99 DT13_1 =
pset.getParameter<std::vector<double> >(
"DT_13_1_scale");
100 DT13_2 =
pset.getParameter<std::vector<double> >(
"DT_13_2_scale");
101 DT14_1 =
pset.getParameter<std::vector<double> >(
"DT_14_1_scale");
102 DT14_2 =
pset.getParameter<std::vector<double> >(
"DT_14_2_scale");
103 DT23_1 =
pset.getParameter<std::vector<double> >(
"DT_23_1_scale");
104 DT23_2 =
pset.getParameter<std::vector<double> >(
"DT_23_2_scale");
105 DT24_1 =
pset.getParameter<std::vector<double> >(
"DT_24_1_scale");
106 DT24_2 =
pset.getParameter<std::vector<double> >(
"DT_24_2_scale");
107 DT34_1 =
pset.getParameter<std::vector<double> >(
"DT_34_1_scale");
108 DT34_2 =
pset.getParameter<std::vector<double> >(
"DT_34_2_scale");
110 OL_1213 =
pset.getParameter<std::vector<double> >(
"OL_1213_0_scale");
111 OL_1222 =
pset.getParameter<std::vector<double> >(
"OL_1222_0_scale");
112 OL_1232 =
pset.getParameter<std::vector<double> >(
"OL_1232_0_scale");
113 OL_2213 =
pset.getParameter<std::vector<double> >(
"OL_2213_0_scale");
114 OL_2222 =
pset.getParameter<std::vector<double> >(
"OL_2222_0_scale");
116 SMB_10S =
pset.getParameter<std::vector<double> >(
"SMB_10_0_scale");
117 SMB_11S =
pset.getParameter<std::vector<double> >(
"SMB_11_0_scale");
118 SMB_12S =
pset.getParameter<std::vector<double> >(
"SMB_12_0_scale");
119 SMB_20S =
pset.getParameter<std::vector<double> >(
"SMB_20_0_scale");
120 SMB_21S =
pset.getParameter<std::vector<double> >(
"SMB_21_0_scale");
121 SMB_22S =
pset.getParameter<std::vector<double> >(
"SMB_22_0_scale");
122 SMB_30S =
pset.getParameter<std::vector<double> >(
"SMB_30_0_scale");
123 SMB_31S =
pset.getParameter<std::vector<double> >(
"SMB_31_0_scale");
124 SMB_32S =
pset.getParameter<std::vector<double> >(
"SMB_32_0_scale");
126 SME_11S =
pset.getParameter<std::vector<double> >(
"SME_11_0_scale");
127 SME_12S =
pset.getParameter<std::vector<double> >(
"SME_12_0_scale");
128 SME_13S =
pset.getParameter<std::vector<double> >(
"SME_13_0_scale");
129 SME_21S =
pset.getParameter<std::vector<double> >(
"SME_21_0_scale");
130 SME_22S =
pset.getParameter<std::vector<double> >(
"SME_22_0_scale");
192 double chi2_dof = 9999.0;
193 unsigned int ini_seg = 0;
197 for (
size_t i = ini_seg;
i < seg.size();
i++) {
198 double dof = static_cast<double>(seg[
i]->degreesOfFreedom());
199 if (chi2_dof < (seg[
i]->
chi2() /
dof))
201 chi2_dof = seg[
i]->chi2() /
dof;
202 best_seg = static_cast<int>(
i);
211 segPos = seg[
last]->localPosition();
217 polar *= fabs(ptmean) / polar.
perp();
220 int chargeI = static_cast<int>(
charge);
223 p_err = (sptmean * sptmean) / (polar.mag() * polar.mag() * ptmean * ptmean);
224 mat = seg[
last]->parametersError().similarityT(seg[
last]->projectionMatrix());
227 mat[0][0] = mat[0][0] / fabs(
tan(mom.
theta()));
228 mat[1][1] = mat[1][1] / fabs(
tan(mom.
theta()));
229 mat[3][3] = 2.25 * mat[3][3];
230 mat[4][4] = 2.25 * mat[4][4];
233 mat[0][0] = mat[0][0] / fabs(
tan(mom.
theta()));
234 mat[1][1] = mat[1][1] / fabs(
tan(mom.
theta()));
235 mat[2][2] = 2.25 * mat[2][2];
236 mat[3][3] = 2.25 * mat[3][3];
237 mat[4][4] = 2.25 * mat[4][4];
239 double dh = fabs(seg[
last]->globalPosition().
eta()) - 1.6;
240 if (fabs(
dh) < 0.1 &&
type == 1) {
241 mat[1][1] = 4. * mat[1][1];
242 mat[2][2] = 4. * mat[2][2];
243 mat[3][3] = 9. * mat[3][3];
244 mat[4][4] = 9. * mat[4][4];
256 segPos = seg[
last]->localPosition();
267 polar *= fabs(ptmean) / polar.
perp();
270 int chargeI = static_cast<int>(
charge);
273 p_err = (sptmean * sptmean) / (polar.mag() * polar.mag() * ptmean * ptmean);
274 mat = seg[
last]->parametersError().similarityT(seg[
last]->projectionMatrix());
281 float Geta =
gp.eta();
284 std::cout <<
"pt " << ptmean <<
" errpt " << sptmean <<
" eta " << Geta <<
" charge " <<
charge << std::endl;
306 for (
unsigned l = 0;
l < seg.size();
l++) {
325 const std::vector<int>&
layers,
328 unsigned size = seg.size();
338 std::vector<double> ptEstimate;
339 std::vector<double> sptEstimate;
349 segPos[0] = seg[0]->globalPosition();
350 float eta = fabs(segPos[0].
eta());
366 unsigned idx1 =
size;
370 int layer1 =
layers[idx1];
371 if (layer0 == layer1)
373 segPos[1] = seg[idx1]->globalPosition();
375 double dphi = segPos[0].
phi() - segPos[1].
phi();
377 double temp_dphi = dphi;
380 if (temp_dphi < 0.) {
381 temp_dphi = -1.0 * temp_dphi;
386 if (temp_dphi < 0.0001) {
390 ptEstimate.push_back(
pt *
sign);
391 sptEstimate.push_back(spt);
394 if (layer0 == 0 && temp_dphi >= 0.0001) {
402 else if (layer1 == 2) {
408 else if (layer1 == 3) {
419 ptEstimate.push_back(
pt *
sign);
420 sptEstimate.push_back(spt);
433 else if (layer1 == 3) {
444 ptEstimate.push_back(
pt *
sign);
445 sptEstimate.push_back(spt);
449 if (layer0 == 2 && temp_dphi > 0.0001) {
468 ptEstimate.push_back(
pt *
sign);
469 sptEstimate.push_back(spt);
473 if (layer0 == 3 && temp_dphi > 0.0001) {
477 ptEstimate.push_back(
pt *
sign);
478 sptEstimate.push_back(spt);
484 if (!ptEstimate.empty())
485 weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
495 const std::vector<int>&
layers,
498 unsigned size = seg.size();
502 std::vector<double> ptEstimate;
503 std::vector<double> sptEstimate;
513 segPos[0] = seg[0]->globalPosition();
514 float eta = fabs(segPos[0].
eta());
523 for (
unsigned idx1 = 1; idx1 <
size; ++idx1) {
524 int layer1 =
layers[idx1];
525 segPos[1] = seg[idx1]->globalPosition();
530 double dphi = segPos[0].
phi() - segPos[1].
phi();
531 double temp_dphi = dphi;
536 if (temp_dphi < 0.) {
537 temp_dphi = -temp_dphi;
541 if (temp_dphi < 0.0001) {
545 ptEstimate.push_back(
pt *
sign);
546 sptEstimate.push_back(spt);
550 if (layer0 == -1 && temp_dphi > 0.0001) {
563 else if (layer1 == -3) {
584 ptEstimate.push_back(
pt *
sign);
585 sptEstimate.push_back(spt);
589 if (layer0 == -2 && temp_dphi > 0.0001) {
612 ptEstimate.push_back(
pt *
sign);
613 sptEstimate.push_back(spt);
617 if (layer0 == -3 && temp_dphi > 0.0001) {
628 ptEstimate.push_back(
pt *
sign);
629 sptEstimate.push_back(spt);
635 if (!ptEstimate.empty())
636 weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
644 const std::vector<int>&
layers,
653 std::vector<int> layersCSC;
655 std::vector<int> layersDT;
658 for (
unsigned j = 0;
j <
layers.size(); ++
j) {
660 segCSC.push_back(seg[
j]);
661 layersCSC.push_back(
layers[
j]);
663 segDT.push_back(seg[
j]);
668 std::vector<double> ptEstimate;
669 std::vector<double> sptEstimate;
673 segPos[0] = seg[0]->globalPosition();
674 float eta = fabs(segPos[0].
eta());
677 if (!segDT.empty() && !segCSC.empty()) {
679 segPos[1] = seg[
size - 1]->globalPosition();
681 double dphi = segPos[0].
phi() - segPos[1].
phi();
682 double temp_dphi = dphi;
687 if (temp_dphi < 0.) {
688 temp_dphi = -temp_dphi;
692 if (temp_dphi < 0.0001) {
696 ptEstimate.push_back(thePt *
sign);
697 sptEstimate.push_back(theSpt);
701 if (layer0 == -1 && temp_dphi > 0.0001) {
709 else if (layer1 == 2) {
720 ptEstimate.push_back(thePt *
sign);
721 sptEstimate.push_back(theSpt);
724 if (layer0 == -2 && temp_dphi > 0.0001) {
730 ptEstimate.push_back(thePt *
sign);
731 sptEstimate.push_back(theSpt);
742 if (segDT.size() > 1) {
744 ptEstimate.push_back(thePt);
745 sptEstimate.push_back(theSpt);
765 if (!ptEstimate.empty())
766 weightedPt(ptEstimate, sptEstimate, thePt, theSpt);
774 const std::vector<int>&
layers,
781 double eta = segPos.
eta();
786 double cosDpsi = (gv.
x() * segPos.
x() + gv.
y() * segPos.
y());
787 cosDpsi /=
sqrt(segPos.
x() * segPos.
x() + segPos.
y() * segPos.
y());
788 cosDpsi /=
sqrt(gv.
x() * gv.
x() + gv.
y() * gv.
y());
790 double axb = (segPos.
x() * gv.
y()) - (segPos.
y() * gv.
x());
791 double sign = (axb < 0.) ? 1.0 : -1.0;
793 double dpsi = fabs(acos(cosDpsi));
794 if (dpsi > 1.570796) {
795 dpsi = 3.141592 - dpsi;
798 if (fabs(dpsi) < 0.00005) {
805 if (fabs(
eta) < 0.3) {
811 if (fabs(
eta) >= 0.3 && fabs(
eta) < 0.82) {
817 if (fabs(
eta) >= 0.82 && fabs(
eta) < 1.2) {
825 if (fabs(
eta) > 0.92 && fabs(
eta) < 1.16) {
831 if (fabs(
eta) >= 1.16 && fabs(
eta) <= 1.6) {
839 if (fabs(
eta) > 1.6) {
848 if (fabs(
eta) < 0.25) {
854 if (fabs(
eta) >= 0.25 && fabs(
eta) < 0.72) {
860 if (fabs(
eta) >= 0.72 && fabs(
eta) < 1.04) {
868 if (fabs(
eta) > 0.95 && fabs(
eta) <= 1.6) {
874 if (fabs(
eta) > 1.6 && fabs(
eta) < 2.45) {
884 if (fabs(
eta) <= 0.22) {
890 if (fabs(
eta) > 0.22 && fabs(
eta) <= 0.6) {
896 if (fabs(
eta) > 0.6 && fabs(
eta) < 0.95) {
902 thePt = fabs(thePt) *
sign;
903 theSpt = fabs(theSpt);
910 if (NShowers > 2 && thePt < 300.) {
914 if (NShowers == 2 && NShowerSegments > 11 && thePt < 150.) {
918 if (NShowers == 2 && NShowerSegments <= 11 && thePt < 50.) {
922 if (NShowers == 1 && NShowerSegments <= 5 && thePt < 10.) {
935 const std::vector<double>& sptEstimate,
938 int size = ptEstimate.size();
942 thePt = ptEstimate[0];
943 theSpt = sptEstimate[0];
951 for (
unsigned j = 0;
j < ptEstimate.size();
j++) {
953 if (ptEstimate[
j] < 0.) {
956 charge -= 1. * (ptEstimate[
j] * ptEstimate[
j]) / (sptEstimate[
j] * sptEstimate[
j]);
959 charge += 1. * (ptEstimate[
j] * ptEstimate[
j]) / (sptEstimate[
j] * sptEstimate[
j]);
971 double weightPtSum = 0.;
972 double sigmaSqr_sum = 0.;
976 for (
unsigned j = 0;
j < ptEstimate.size(); ++
j) {
980 sigmaSqr_sum += 1.0 / (sptEstimate[
j] * sptEstimate[
j]);
981 weightPtSum += fabs(ptEstimate[
j]) / (sptEstimate[
j] * sptEstimate[
j]);
993 thePt = (
charge * weightPtSum) / sigmaSqr_sum;
994 theSpt =
sqrt(1.0 / sigmaSqr_sum);
1002 double h = fabs(
eta);
1003 double estPt = (vPara[0] + vPara[1] *
h + vPara[2] *
h *
h) /
dPhi;
1004 double estSPt = (vPara[3] + vPara[4] *
h + vPara[5] *
h *
h) * estPt;
1005 std::vector<double> paraPt;
1006 paraPt.push_back(estPt);
1007 paraPt.push_back(estSPt);
1015 double oPhi = 1. / dphi;
1016 dphi = dphi / (1. +
t1 / (oPhi + 10.));