61 static bool pair_compare(std::pair<const Stub*, float>
a, std::pair<const Stub*, float>
b) {
62 return (
a.second <
b.second);
67 PrintL1trk() <<
"=============== FITTING SimpleLR TRACK ====================";
90 double tanLambda = 0.;
100 unsigned int numStubs = 0;
109 SumRPhi = SumRPhi + digiStub->
rt_SF_TF() * digiStub->
phiS();
111 SumPhi = SumPhi + digiStub->
phiS();
115 <<
" z " << digiStub->
iDigi_Z();
118 if (l1track3D.
iPhiSec() == 0 and stub->phi() > 0) {
119 phi = stub->phi() - 2 *
M_PI;
121 phi = stub->phi() + 2 *
M_PI;
125 SumRPhi = SumRPhi + stub->r() * phi;
126 SumR = SumR + stub->r();
127 SumPhi = SumPhi + phi;
128 SumR2 = SumR2 + stub->r() * stub->r();
130 PrintL1trk() <<
"InputStub (float): phi " << phi <<
" r " << stub->r() <<
" z " << stub->z();
134 double numeratorPt, digiNumeratorPt;
136 double numeratorPhi, digiNumeratorPhi;
137 double reciprocal, digiReciprocal;
139 double numeratorLambda;
141 digiNumeratorPt = (numStubs * SumRPhi - SumR * SumPhi);
142 digiDenominator = (numStubs * SumR2 - SumR * SumR);
143 digiNumeratorPhi = (SumR2 * SumPhi - SumR * SumRPhi);
146 qOverPt = (numStubs * SumRPhi - SumR * SumPhi) / (numStubs * SumR2 - SumR * SumR);
147 phi0 = (SumR2 * SumPhi - SumR * SumRPhi) / (numStubs * SumR2 - SumR * SumR);
173 PrintL1trk() << setw(10) <<
"Input helix (digi): qOverPt = " << qOverPt <<
" (" << floor(qOverPt *
qOverPtMult_)
174 <<
"), phiT = " << phiT <<
" (" << floor(phiT *
phiTMult_) <<
") ";
176 PrintL1trk() <<
"Input Helix (float): qOverPt = " << qOverPt <<
" phi0 " << phi0;
181 std::vector<std::pair<Stub*, double> >
vRes;
182 unsigned int psStubs = 0;
184 if (stub->psModule())
206 std::pair<Stub*, double> ResStubPair(stub, Res);
207 vRes.push_back(ResStubPair);
209 if (stub->assocTP() !=
nullptr)
210 PrintL1trk() <<
" Stub rphi residual " << Res <<
" TP " << stub->assocTP()->index();
212 PrintL1trk() <<
" Stub rphi residual " << Res <<
" TP nullptr";
216 double largestResidual = 9999.;
219 std::vector<std::pair<Stub*, double> >::iterator maxResIt = max_element(
vRes.begin(),
vRes.end(),
pair_compare);
220 largestResidual = (*maxResIt).second;
222 PrintL1trk() <<
"Largest Residual " << largestResidual;
225 if ((*maxResIt).first->psModule()) {
228 PrintL1trk() <<
"removing PS residual " << (*maxResIt).second;
229 vRes.erase(maxResIt);
233 PrintL1trk() <<
"residual " << (*maxResIt).second <<
" set to -1. ";
234 (*maxResIt).second = -1.;
237 vRes.erase(maxResIt);
239 PrintL1trk() <<
"removing residual " << (*maxResIt).second;
244 std::vector<Stub*> fitStubs;
245 fitStubs.reserve(
vRes.size());
246 for (std::pair<Stub*, double> ResStubPair :
vRes) {
247 fitStubs.push_back(ResStubPair.first);
260 double SumR2_ps = 0.;
265 for (
const Stub* stub : fitStubs) {
266 if (stub->psModule())
274 SumPhi += digiStub->
phiS();
276 if (stub->psModule()) {
277 SumRZ += digiStub->
rt_SF_TF() * digiStub->
z();
278 SumZ += digiStub->
z();
288 if (l1track3D.
iPhiSec() == 0 and stub->phi() > 0) {
289 phi = stub->phi() - 2 *
M_PI;
291 phi = stub->phi() + 2 *
M_PI;
296 SumRPhi += stub->r() * phi;
299 SumR2 += stub->r() * stub->r();
300 if (stub->psModule()) {
301 SumRZ += stub->r() * stub->z();
303 SumR_ps += stub->r();
304 SumR2_ps += stub->r() * stub->r();
307 PrintL1trk() <<
"phi " << phi <<
" r " << stub->r() <<
" z " << stub->z();
311 numeratorZ0 = (SumR2_ps * SumZ - SumR_ps * SumRZ);
312 numeratorLambda = (psStubs * SumRZ - SumR_ps * SumZ);
313 numeratorPt = (numStubs * SumRPhi - SumR * SumPhi);
315 double denominatorZ = (psStubs * SumR2_ps - SumR_ps * SumR_ps);
316 numeratorPhi = (SumR2 * SumPhi - SumR * SumRPhi);
319 z0 = numeratorZ0 / denominatorZ;
320 tanLambda = numeratorLambda / denominatorZ;
321 qOverPt = (numStubs * SumRPhi - SumR * SumPhi) / (numStubs * SumR2 - SumR * SumR);
322 phi0 = (SumR2 * SumPhi - SumR * SumRPhi) / (numStubs * SumR2 - SumR * SumR);
364 PrintL1trk() <<
"Second Helix variables: numeratorPt = " << numeratorPt <<
", numeratorPhi = " << numeratorPhi
365 <<
", numeratorZ0 = " << numeratorZ0 <<
" numeratorLambda = " << numeratorLambda
366 <<
" denominator = " <<
denominator <<
" reciprocal = " << reciprocal
367 <<
" denominatorZ = " << denominatorZ <<
" reciprocalZ = " << reciprocalZ;
368 PrintL1trk() << setw(10) <<
"Final Helix parameters: qOverPt = " << qOverPt <<
" ("
370 <<
"), zT = " << zT <<
" (" << floor(zT *
z0Mult_) <<
"), tanLambda = " << tanLambda <<
" ("
374 PrintL1trk() << setw(10) <<
"Final Helix parameters: qOverPt = " << qOverPt <<
", phi0 = " << phi0
375 <<
", z0 = " <<
z0 <<
", tanLambda = " << tanLambda;
378 double chi2_phi = 0.;
381 for (
const Stub* stub : fitStubs) {
386 ResPhi = digiStub->
phiS() - phiT - qOverPt * digiStub->
rt_SF_TF();
387 ResZ = digiStub->
z() - zT - tanLambda * digiStub->
rt_SF_TF();
390 ResZ = stub->z() -
z0 - tanLambda * stub->r();
393 double RPhiSigma = 0.0002;
394 float RZSigma = stub->sigmaZ() +
std::abs(tanLambda) * stub->sigmaR();
396 if (not stub->barrel())
406 chi2_phi +=
std::abs(ResPhi * ResPhi);
409 PrintL1trk() <<
"Stub ResPhi " << ResPhi * RPhiSigma <<
" ResSigma " << RPhiSigma <<
" Res " << ResPhi
410 <<
" chi2 " << chi2_phi;
411 PrintL1trk() <<
"Stub ResZ " << ResZ * RZSigma <<
" ResSigma " << RZSigma <<
" Res " << ResZ <<
" chi2 "
417 bool accepted =
false;
420 double chi2 = chi2_phi + chi2_z;
424 constexpr
unsigned int nHelixPar = 4;
425 float dof = 2 * fitStubs.size() - nHelixPar;
431 PrintL1trk() <<
"qOverPt " << qOverPt <<
" phiT " << phiT;
446 const unsigned int hitPattern = 0;
448 settings_, &l1track3D, fitStubs, hitPattern, qOverPt, 0., phi0,
z0, tanLambda, chi2_phi, chi2_z, nHelixPar);
458 PrintL1trk() << setw(10) <<
"First Helix parameters: qOverPt = " << fitTrk.
qOverPt() <<
" oneOver2r "
462 <<
"), zT = " << zT <<
" (" << floor(zT *
z0Mult_) <<
"), tanLambda = " << tanLambda <<
" ("
472 PrintL1trk() <<
"True track: chi2/ndf " << chi2dof;
475 PrintL1trk() <<
"Track rejected " <<
chi2 <<
" chi2/ndof " << chi2dof;
477 PrintL1trk() <<
"Fake track!!! " <<
chi2 <<
" chi2/ndof " << chi2dof;