67 PrintL1trk() <<
"=============== FITTING SimpleLR TRACK ====================";
90 double tanLambda = 0.;
100 unsigned int numStubs = 0;
103 for (Stub* stub : l1track3D.stubs()) {
107 const DigitalStub* digiStub = stub->digitalStub();
109 SumRPhi = SumRPhi + digiStub->rt_SF_TF() * digiStub->phiS();
110 SumR = SumR + digiStub->rt_SF_TF();
111 SumPhi = SumPhi + digiStub->phiS();
112 SumR2 = SumR2 + digiStub->rt_SF_TF() * digiStub->rt_SF_TF();
114 PrintL1trk() <<
"Input stub (digi): phiS " << digiStub->iDigi_PhiS() <<
" rT " << digiStub->iDigi_Rt()
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;
183 for (Stub* stub : l1track3D.stubs()) {
184 if (stub->psModule())
189 const DigitalStub* digiStub = stub->digitalStub();
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())
271 const DigitalStub* digiStub = stub->digitalStub();
272 SumRPhi += digiStub->rt_SF_TF() * digiStub->phiS();
273 SumR += digiStub->rt_SF_TF();
274 SumPhi += digiStub->phiS();
275 SumR2 += digiStub->rt_SF_TF() * digiStub->rt_SF_TF();
276 if (stub->psModule()) {
277 SumRZ += digiStub->rt_SF_TF() * digiStub->z();
278 SumZ += digiStub->z();
279 SumR_ps += digiStub->rt_SF_TF();
280 SumR2_ps += digiStub->rt_SF_TF() * digiStub->rt_SF_TF();
283 PrintL1trk() <<
"phiS " << digiStub->iDigi_PhiS() <<
" rT " << digiStub->iDigi_Rt() <<
" z "
284 << digiStub->iDigi_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);
314 denominator = (numStubs * SumR2 - SumR * SumR);
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);
361 PrintL1trk() <<
"HT mbin " << int(l1track3D.cellLocationHT().first) - 16 <<
" cbin "
362 <<
int(l1track3D.cellLocationHT().second) - 32 <<
" iPhi " << l1track3D.iPhiSec() <<
" iEta "
363 << l1track3D.iEtaReg();
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) {
385 const DigitalStub* digiStub = stub->digitalStub();
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;
426 float chi2dof = chi2 /
dof;
431 PrintL1trk() <<
"qOverPt " << qOverPt <<
" phiT " << phiT;
446 const unsigned int hitPattern = 0;
447 L1fittedTrack fitTrk(
448 settings_, &l1track3D, fitStubs, hitPattern, qOverPt, 0., phi0, z0, tanLambda, chi2_phi, chi2_z, nHelixPar);
451 fitTrk.digitizeTrack(
"SimpleLR4");
453 if (
debug_ and digitize_) {
454 PrintL1trk() <<
"Digitized parameters ";
455 PrintL1trk() <<
"HT mbin " << int(l1track3D.cellLocationHT().first) - 16 <<
" cbin "
456 <<
int(l1track3D.cellLocationHT().second) - 32 <<
" iPhi " << l1track3D.iPhiSec() <<
" iEta "
457 << l1track3D.iEtaReg();
458 PrintL1trk() << setw(10) <<
"First Helix parameters: qOverPt = " << fitTrk.qOverPt() <<
" oneOver2r "
459 << fitTrk.digitaltrack()->oneOver2r() <<
" ("
460 << floor(fitTrk.digitaltrack()->oneOver2r() *
qOverPtMult_)
461 <<
"), phi0 = " << fitTrk.digitaltrack()->phi0() <<
" (" << fitTrk.digitaltrack()->iDigi_phi0rel()
462 <<
"), zT = " << zT <<
" (" << floor(zT *
z0Mult_) <<
"), tanLambda = " << tanLambda <<
" ("
467 PrintL1trk() <<
"FitTrack helix parameters " << int(fitTrk.cellLocationFit().first) - 16 <<
", "
468 <<
int(fitTrk.cellLocationFit().second) - 32 <<
" HT parameters "
469 <<
int(fitTrk.cellLocationHT().first) - 16 <<
", " <<
int(fitTrk.cellLocationHT().second) - 32;
471 if (fitTrk.matchedTP() !=
nullptr) {
472 PrintL1trk() <<
"True track: chi2/ndf " << chi2dof;
473 PrintL1trk() <<
"TP qOverPt " << fitTrk.matchedTP()->qOverPt() <<
" phi0 " << fitTrk.matchedTP()->phi0();
475 PrintL1trk() <<
"Track rejected " << chi2 <<
" chi2/ndof " << chi2dof;
477 PrintL1trk() <<
"Fake track!!! " << chi2 <<
" chi2/ndof " << chi2dof;
479 PrintL1trk() <<
"layers in track " << fitTrk.numLayers();
485 L1fittedTrack rejectedTrk;
bool enableDigitize() const
constexpr double deltaPhi(double phi1, double phi2)
unsigned int shiftingBitsz0_
unsigned int shiftingBitsPhi_
unsigned int phiSBits() const
double houghMinPt() const
unsigned int shiftingBitsDenRZ_
float numeratorLambdaMult_
unsigned int minStubLayersRed_
unsigned int shiftingBitsLambda_
static bool pair_compare(std::pair< const Stub *, float > a, std::pair< const Stub *, float > b)
double chosenRofPhi() const
Abs< T >::type abs(const T &t)
const Settings * settings_
unsigned int dividerBitsHelix_
unsigned int dividerBitsHelixZ_
double invPtToDphi() const
unsigned int shiftingBitsPt_
unsigned int shiftingBitsDenRPhi_
double ResidualCut() const
Cut on RPhi Residual (radians)
unsigned int numPhiSectors() const
unsigned int slr_phi0Bits() const
unsigned int numLayerCut(Utility::AlgoStep algo, const Settings *settings, unsigned int iPhiSec, unsigned int iEtaReg, float invPt, float eta=0.)
Power< A, B >::type pow(const A &a, const B &b)