diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx index 94ed00e1601e1cc258a563da52bdbb2235ea07a4..b2b96670a222cf45e7e6c008e9e71ab298b81dc5 100644 --- a/reco/L1/L1Algo/L1TrackFitter.cxx +++ b/reco/L1/L1Algo/L1TrackFitter.cxx @@ -333,8 +333,10 @@ void L1Algo::L1KFTrackFitter() const int nStations = fParameters.GetNstationsActive(); int nTracks_SIMD = fvec::size(); - L1TrackParFit T1; // fitting parametr coresponding to current track - T1.SetParticleMass(GetDefaultParticleMass()); + L1TrackParFit fit; // fitting parametr coresponding to current track + L1TrackPar& tr = fit.fTr; + + fit.SetParticleMass(GetDefaultParticleMass()); L1Track* t[fvec::size()] {nullptr}; @@ -491,15 +493,13 @@ void L1Algo::L1KFTrackFitter() } } - GuessVec(T1, x, y, z, Sy, w, nStations, &z_end, time, w_time); + GuessVec(fit, x, y, z, Sy, w, nStations, &z_end, time, w_time); - if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { - T1.fqp = fvec(0.); - } + if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { tr.qp = fvec(0.); } for (int iter = 0; iter < 2; iter++) { // 1.5 iterations - fvec qp01 = T1.fqp; + fvec qp01 = tr.qp; // fit backward @@ -508,17 +508,17 @@ void L1Algo::L1KFTrackFitter() time_last = iif(w_time[ista] > fvec::Zero(), time_last, fvec::Zero()); time_er_last = iif(w_time[ista] > fvec::Zero(), time_er_last, fvec(100.)); - FilterFirst(T1, x_last, y_last, time_last, time_er_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst); + FilterFirst(fit, x_last, y_last, time_last, time_er_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst); fldZ1 = z[ista]; - sta[ista].fieldSlice.GetFieldValue(T1.fx, T1.fy, fldB1); + sta[ista].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1); fldB1.Combine(fB[ista], w[ista]); fldZ2 = z[ista - 2]; fvec dz = fldZ2 - fldZ1; - sta[ista].fieldSlice.GetFieldValue(T1.fx + T1.ftx * dz, T1.fy + T1.fty * dz, fldB2); + sta[ista].fieldSlice.GetFieldValue(tr.x + tr.tx * dz, tr.y + tr.ty * dz, fldB2); fldB2.Combine(fB[ista - 2], w[ista - 2]); fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0); @@ -526,7 +526,7 @@ void L1Algo::L1KFTrackFitter() fldZ0 = z[ista]; dz = (fldZ1 - fldZ0); - sta[ista].fieldSlice.GetFieldValue(T1.fx - T1.ftx * dz, T1.fy - T1.fty * dz, fldB0); + sta[ista].fieldSlice.GetFieldValue(tr.x - tr.tx * dz, tr.y - tr.ty * dz, fldB0); fldB0.Combine(fB[ista], w[ista]); fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2); @@ -537,30 +537,30 @@ void L1Algo::L1KFTrackFitter() fld1 = fld; - T1.Extrapolate(z[ista], qp01, fld1, wExtr); + fit.Extrapolate(z[ista], qp01, fld1, wExtr); if (ista == fNstationsBeforePipe - 1) { - T1.L1AddPipeMaterial(qp01, wExtr); - T1.EnergyLossCorrection(T1.fPipeRadThick, qp01, fvec(1.f), wExtr); + fit.L1AddPipeMaterial(qp01, wExtr); + fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(1.f), wExtr); } if constexpr (L1Constants::control::kIfUseRadLengthTable) { - T1.L1AddMaterial(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, wExtr); - T1.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, fvec(1.f), wExtr); + fit.L1AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr); + fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(1.f), wExtr); } else { - T1.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr); + fit.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr); } L1UMeasurementInfo info = sta[ista].frontInfo; info.sigma2 = d_u[ista] * d_u[ista]; - T1.Filter(info, u[ista], w1); + fit.Filter(info, u[ista], w1); info = sta[ista].backInfo; info.sigma2 = d_v[ista] * d_v[ista]; - T1.Filter(info, v[ista], w1); - T1.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo); + fit.Filter(info, v[ista], w1); + fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo); fldB2 = fldB1; fldZ2 = fldZ1; @@ -570,7 +570,7 @@ void L1Algo::L1KFTrackFitter() // extrapolate to the PV region - L1TrackParFit T1pv = T1; + L1TrackParFit fitpv = fit; { L1UMeasurementInfo vtxInfoX; vtxInfoX.cos_phi = 1.; @@ -583,34 +583,34 @@ void L1Algo::L1KFTrackFitter() vtxInfoY.sigma2 = fvec(1.e-4); if (kGlobal == fTrackingMode) { - fvec qp00 = T1.fqp; + fvec qp00 = tr.qp; for (int vtxIter = 0; vtxIter < 2; vtxIter++) { - T1pv = T1; - T1pv.fqp = qp00; - T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fqp, fld, fvec(1.)); - T1pv.Filter(vtxInfoX, fParameters.GetTargetPositionX()); - T1pv.Filter(vtxInfoY, fParameters.GetTargetPositionY()); - qp00 = T1pv.fqp; + fitpv = fit; + fitpv.fTr.qp = qp00; + fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fitpv.fTr.qp, fld, fvec(1.)); + fitpv.Filter(vtxInfoX, fParameters.GetTargetPositionX()); + fitpv.Filter(vtxInfoY, fParameters.GetTargetPositionY()); + qp00 = fitpv.fTr.qp; } } else { - T1pv = T1; - T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fqp, fld, fvec(1.)); + fitpv = fit; + fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fitpv.fTr.qp, fld, fvec(1.)); } } - L1TrackParFit Tf = T1; - if (kGlobal == fTrackingMode) { Tf = T1pv; } + L1TrackPar Tf = fit.fTr; + if (kGlobal == fTrackingMode) { Tf = fitpv.fTr; } for (int iVec = 0; iVec < nTracks_SIMD; iVec++) { - t[iVec]->TFirst[0] = Tf.fx[iVec]; - t[iVec]->TFirst[1] = Tf.fy[iVec]; - t[iVec]->TFirst[2] = Tf.ftx[iVec]; - t[iVec]->TFirst[3] = Tf.fty[iVec]; - t[iVec]->TFirst[4] = Tf.fqp[iVec]; - if (kGlobal == fTrackingMode) { t[iVec]->TFirst[4] = T1pv.fqp[iVec]; } - t[iVec]->TFirst[5] = Tf.fz[iVec]; - t[iVec]->TFirst[6] = Tf.ft[iVec]; + t[iVec]->TFirst[0] = Tf.x[iVec]; + t[iVec]->TFirst[1] = Tf.y[iVec]; + t[iVec]->TFirst[2] = Tf.tx[iVec]; + t[iVec]->TFirst[3] = Tf.ty[iVec]; + t[iVec]->TFirst[4] = Tf.qp[iVec]; + if (kGlobal == fTrackingMode) { t[iVec]->TFirst[4] = fitpv.fTr.qp[iVec]; } + t[iVec]->TFirst[5] = Tf.z[iVec]; + t[iVec]->TFirst[6] = Tf.t[iVec]; t[iVec]->CFirst[0] = Tf.C00[iVec]; t[iVec]->CFirst[1] = Tf.C10[iVec]; @@ -639,35 +639,35 @@ void L1Algo::L1KFTrackFitter() } for (int iVec = 0; iVec < nTracks_SIMD; iVec++) { - t[iVec]->Tpv[0] = T1pv.fx[iVec]; - t[iVec]->Tpv[1] = T1pv.fy[iVec]; - t[iVec]->Tpv[2] = T1pv.ftx[iVec]; - t[iVec]->Tpv[3] = T1pv.fty[iVec]; - t[iVec]->Tpv[4] = T1pv.fqp[iVec]; - t[iVec]->Tpv[5] = T1pv.fz[iVec]; - t[iVec]->Tpv[6] = T1pv.ft[iVec]; - - t[iVec]->Cpv[0] = T1pv.C00[iVec]; - t[iVec]->Cpv[1] = T1pv.C10[iVec]; - t[iVec]->Cpv[2] = T1pv.C11[iVec]; - t[iVec]->Cpv[3] = T1pv.C20[iVec]; - t[iVec]->Cpv[4] = T1pv.C21[iVec]; - t[iVec]->Cpv[5] = T1pv.C22[iVec]; - t[iVec]->Cpv[6] = T1pv.C30[iVec]; - t[iVec]->Cpv[7] = T1pv.C31[iVec]; - t[iVec]->Cpv[8] = T1pv.C32[iVec]; - t[iVec]->Cpv[9] = T1pv.C33[iVec]; - t[iVec]->Cpv[10] = T1pv.C40[iVec]; - t[iVec]->Cpv[11] = T1pv.C41[iVec]; - t[iVec]->Cpv[12] = T1pv.C42[iVec]; - t[iVec]->Cpv[13] = T1pv.C43[iVec]; - t[iVec]->Cpv[14] = T1pv.C44[iVec]; - t[iVec]->Cpv[15] = T1pv.C50[iVec]; - t[iVec]->Cpv[16] = T1pv.C51[iVec]; - t[iVec]->Cpv[17] = T1pv.C52[iVec]; - t[iVec]->Cpv[18] = T1pv.C53[iVec]; - t[iVec]->Cpv[19] = T1pv.C54[iVec]; - t[iVec]->Cpv[20] = T1pv.C55[iVec]; + t[iVec]->Tpv[0] = fitpv.fTr.x[iVec]; + t[iVec]->Tpv[1] = fitpv.fTr.y[iVec]; + t[iVec]->Tpv[2] = fitpv.fTr.tx[iVec]; + t[iVec]->Tpv[3] = fitpv.fTr.ty[iVec]; + t[iVec]->Tpv[4] = fitpv.fTr.qp[iVec]; + t[iVec]->Tpv[5] = fitpv.fTr.z[iVec]; + t[iVec]->Tpv[6] = fitpv.fTr.t[iVec]; + + t[iVec]->Cpv[0] = fitpv.fTr.C00[iVec]; + t[iVec]->Cpv[1] = fitpv.fTr.C10[iVec]; + t[iVec]->Cpv[2] = fitpv.fTr.C11[iVec]; + t[iVec]->Cpv[3] = fitpv.fTr.C20[iVec]; + t[iVec]->Cpv[4] = fitpv.fTr.C21[iVec]; + t[iVec]->Cpv[5] = fitpv.fTr.C22[iVec]; + t[iVec]->Cpv[6] = fitpv.fTr.C30[iVec]; + t[iVec]->Cpv[7] = fitpv.fTr.C31[iVec]; + t[iVec]->Cpv[8] = fitpv.fTr.C32[iVec]; + t[iVec]->Cpv[9] = fitpv.fTr.C33[iVec]; + t[iVec]->Cpv[10] = fitpv.fTr.C40[iVec]; + t[iVec]->Cpv[11] = fitpv.fTr.C41[iVec]; + t[iVec]->Cpv[12] = fitpv.fTr.C42[iVec]; + t[iVec]->Cpv[13] = fitpv.fTr.C43[iVec]; + t[iVec]->Cpv[14] = fitpv.fTr.C44[iVec]; + t[iVec]->Cpv[15] = fitpv.fTr.C50[iVec]; + t[iVec]->Cpv[16] = fitpv.fTr.C51[iVec]; + t[iVec]->Cpv[17] = fitpv.fTr.C52[iVec]; + t[iVec]->Cpv[18] = fitpv.fTr.C53[iVec]; + t[iVec]->Cpv[19] = fitpv.fTr.C54[iVec]; + t[iVec]->Cpv[20] = fitpv.fTr.C55[iVec]; } if (iter == 1) continue; // only 1.5 iterations @@ -676,24 +676,24 @@ void L1Algo::L1KFTrackFitter() ista = 0; - FilterFirst(T1, x_first, y_first, time_first, time_er_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst); + FilterFirst(fit, x_first, y_first, time_first, time_er_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst); - qp01 = T1.fqp; + qp01 = tr.qp; fldZ1 = z[ista]; - sta[ista].fieldSlice.GetFieldValue(T1.fx, T1.fy, fldB1); + sta[ista].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1); fldB1.Combine(fB[ista], w[ista]); fldZ2 = z[ista + 2]; dz = fldZ2 - fldZ1; - sta[ista].fieldSlice.GetFieldValue(T1.fx + T1.ftx * dz, T1.fy + T1.fty * dz, fldB2); + sta[ista].fieldSlice.GetFieldValue(tr.x + tr.tx * dz, tr.y + tr.ty * dz, fldB2); fldB2.Combine(fB[ista + 2], w[ista + 2]); fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0); for (++ista; ista < nStations; ista++) { fldZ0 = z[ista]; dz = (fldZ1 - fldZ0); - sta[ista].fieldSlice.GetFieldValue(T1.fx - T1.ftx * dz, T1.fy - T1.fty * dz, fldB0); + sta[ista].fieldSlice.GetFieldValue(tr.x - tr.tx * dz, tr.y - tr.ty * dz, fldB0); fldB0.Combine(fB[ista], w[ista]); fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2); @@ -702,30 +702,30 @@ void L1Algo::L1KFTrackFitter() fvec w1_time = iif(initialised, w_time[ista], fvec::Zero()); fvec wExtr = iif(initialised, fvec::One(), fvec::Zero()); - T1.Extrapolate(z[ista], qp01, fld, w1); + fit.Extrapolate(z[ista], qp01, fld, w1); if (ista == fNstationsBeforePipe) { - T1.L1AddPipeMaterial(qp01, wExtr); - T1.EnergyLossCorrection(T1.fPipeRadThick, qp01, fvec(-1.f), wExtr); + fit.L1AddPipeMaterial(qp01, wExtr); + fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(-1.f), wExtr); } if constexpr (L1Constants::control::kIfUseRadLengthTable) { - T1.L1AddMaterial(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, wExtr); - T1.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, fvec(-1.f), wExtr); + fit.L1AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr); + fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(-1.f), wExtr); } else { - T1.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr); + fit.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr); } L1UMeasurementInfo info = sta[ista].frontInfo; info.sigma2 = d_u[ista] * d_u[ista]; - T1.Filter(info, u[ista], w1); + fit.Filter(info, u[ista], w1); info = sta[ista].backInfo; info.sigma2 = d_v[ista] * d_v[ista]; - T1.Filter(info, v[ista], w1); - T1.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo); + fit.Filter(info, v[ista], w1); + fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo); fldB2 = fldB1; fldZ2 = fldZ1; @@ -734,39 +734,39 @@ void L1Algo::L1KFTrackFitter() } for (int iVec = 0; iVec < nTracks_SIMD; iVec++) { - t[iVec]->TLast[0] = T1.fx[iVec]; - t[iVec]->TLast[1] = T1.fy[iVec]; - t[iVec]->TLast[2] = T1.ftx[iVec]; - t[iVec]->TLast[3] = T1.fty[iVec]; - t[iVec]->TLast[4] = T1.fqp[iVec]; - if (kGlobal == fTrackingMode) { t[iVec]->TLast[4] = T1pv.fqp[iVec]; } - t[iVec]->TLast[5] = T1.fz[iVec]; - t[iVec]->TLast[6] = T1.ft[iVec]; - - t[iVec]->CLast[0] = T1.C00[iVec]; - t[iVec]->CLast[1] = T1.C10[iVec]; - t[iVec]->CLast[2] = T1.C11[iVec]; - t[iVec]->CLast[3] = T1.C20[iVec]; - t[iVec]->CLast[4] = T1.C21[iVec]; - t[iVec]->CLast[5] = T1.C22[iVec]; - t[iVec]->CLast[6] = T1.C30[iVec]; - t[iVec]->CLast[7] = T1.C31[iVec]; - t[iVec]->CLast[8] = T1.C32[iVec]; - t[iVec]->CLast[9] = T1.C33[iVec]; - t[iVec]->CLast[10] = T1.C40[iVec]; - t[iVec]->CLast[11] = T1.C41[iVec]; - t[iVec]->CLast[12] = T1.C42[iVec]; - t[iVec]->CLast[13] = T1.C43[iVec]; - t[iVec]->CLast[14] = T1.C44[iVec]; - t[iVec]->CLast[15] = T1.C50[iVec]; - t[iVec]->CLast[16] = T1.C51[iVec]; - t[iVec]->CLast[17] = T1.C52[iVec]; - t[iVec]->CLast[18] = T1.C53[iVec]; - t[iVec]->CLast[19] = T1.C54[iVec]; - t[iVec]->CLast[20] = T1.C55[iVec]; - - t[iVec]->chi2 = T1.chi2[iVec]; - t[iVec]->NDF = (int) T1.NDF[iVec]; + t[iVec]->TLast[0] = tr.x[iVec]; + t[iVec]->TLast[1] = tr.y[iVec]; + t[iVec]->TLast[2] = tr.tx[iVec]; + t[iVec]->TLast[3] = tr.ty[iVec]; + t[iVec]->TLast[4] = tr.qp[iVec]; + if (kGlobal == fTrackingMode) { t[iVec]->TLast[4] = fitpv.fTr.qp[iVec]; } + t[iVec]->TLast[5] = tr.z[iVec]; + t[iVec]->TLast[6] = tr.t[iVec]; + + t[iVec]->CLast[0] = tr.C00[iVec]; + t[iVec]->CLast[1] = tr.C10[iVec]; + t[iVec]->CLast[2] = tr.C11[iVec]; + t[iVec]->CLast[3] = tr.C20[iVec]; + t[iVec]->CLast[4] = tr.C21[iVec]; + t[iVec]->CLast[5] = tr.C22[iVec]; + t[iVec]->CLast[6] = tr.C30[iVec]; + t[iVec]->CLast[7] = tr.C31[iVec]; + t[iVec]->CLast[8] = tr.C32[iVec]; + t[iVec]->CLast[9] = tr.C33[iVec]; + t[iVec]->CLast[10] = tr.C40[iVec]; + t[iVec]->CLast[11] = tr.C41[iVec]; + t[iVec]->CLast[12] = tr.C42[iVec]; + t[iVec]->CLast[13] = tr.C43[iVec]; + t[iVec]->CLast[14] = tr.C44[iVec]; + t[iVec]->CLast[15] = tr.C50[iVec]; + t[iVec]->CLast[16] = tr.C51[iVec]; + t[iVec]->CLast[17] = tr.C52[iVec]; + t[iVec]->CLast[18] = tr.C53[iVec]; + t[iVec]->CLast[19] = tr.C54[iVec]; + t[iVec]->CLast[20] = tr.C55[iVec]; + + t[iVec]->chi2 = tr.chi2[iVec]; + t[iVec]->NDF = (int) tr.NDF[iVec]; } } // iter } @@ -792,6 +792,7 @@ void L1Algo::L1KFTrackFitterMuch() L1TrackPar T; // fitting parametr coresponding to current track L1TrackParFit T1; // fitting parametr coresponding to current track + L1TrackPar& tr = T1.fTr; T1.SetParticleMass(GetDefaultParticleMass()); L1Fit fit; @@ -970,7 +971,7 @@ void L1Algo::L1KFTrackFitterMuch() fvec qp0 = T.qp; - fvec qp01 = T1.fqp; + fvec qp01 = tr.qp; i = 0; @@ -1024,10 +1025,10 @@ void L1Algo::L1KFTrackFitterMuch() fldZ1 = fldZ0; if constexpr (L1Constants::control::kIfUseRadLengthTable) { - T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, fvec(-1.f), wIn); + T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(-1.f), wIn); - T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, wIn, - sta[i].materialInfo.thick, 1); + T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick, + 1); } else { // L1AddMaterial( T, sta[i].materialInfo, qp0, wIn ); @@ -1050,10 +1051,10 @@ void L1Algo::L1KFTrackFitterMuch() fvec z_last = z[i]; - // T1.ExtrapolateLine( T1.fz + 10, &w1); + // T1.ExtrapolateLine( tr.z + 10, &w1); // L1ExtrapolateLine( T, T.z + 10); - fvec d_z = T1.fz - z_last; + fvec d_z = tr.z - z_last; const fvec st = fvec(10); @@ -1070,14 +1071,14 @@ void L1Algo::L1KFTrackFitterMuch() fvec nofSteps1 = fvec::Zero(); fvec stepSize = iif(nofSteps < fvec::One(), d_z, st) * w1; - fvec z_cur = T1.fz; + fvec z_cur = tr.z; fvec w2 = w1; for (int iStep = 0; iStep < max_steps + 1; iStep++) { const fmask maskLastStep = (nofSteps == nofSteps1); - z_cur = iif(maskLastStep, z_last, T1.fz + stepSize); + z_cur = iif(maskLastStep, z_last, tr.z + stepSize); // fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01)); // T1.ExtrapolateLine1( z, &w2, v_mc); @@ -1091,16 +1092,16 @@ void L1Algo::L1KFTrackFitterMuch() // TODO: Unify energy loss corrections (S.Zharko) if constexpr (L1Constants::control::kIfUseRadLengthTable) { if (i == 11 || i == 14 || i == 17) - T1.EnergyLossCorrectionIron(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + 1), qp01, + T1.EnergyLossCorrectionIron(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01, fvec(-1.f), wIn); if (i == 8) - T1.EnergyLossCorrectionCarbon(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + 1), qp01, + T1.EnergyLossCorrectionCarbon(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01, fvec(-1.f), wIn); if (i == 9 || i == 10 || i == 12 || i == 13 || i == 15 || i == 16 || i == 18 || i == 19) - T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + 1), qp01, + T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01, fvec(-1.f), wIn); - T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + fvec(1)), qp01, wIn, + T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, wIn, sta[i].materialInfo.thick / (nofSteps + fvec(1)), 1); wIn = iif(!maskLastStep, wIn, fvec::Zero()); @@ -1128,39 +1129,39 @@ void L1Algo::L1KFTrackFitterMuch() for (iVec = 0; iVec < nTracks_SIMD; iVec++) { - t[iVec]->TLast[0] = T1.fx[iVec]; - t[iVec]->TLast[1] = T1.fy[iVec]; - t[iVec]->TLast[2] = T1.ftx[iVec]; - t[iVec]->TLast[3] = T1.fty[iVec]; - t[iVec]->TLast[4] = T1.fqp[iVec]; - t[iVec]->TLast[5] = T1.fz[iVec]; - t[iVec]->TLast[6] = T1.ft[iVec]; - - t[iVec]->CLast[0] = T1.C00[iVec]; - t[iVec]->CLast[1] = T1.C10[iVec]; - t[iVec]->CLast[2] = T1.C11[iVec]; - t[iVec]->CLast[3] = T1.C20[iVec]; - t[iVec]->CLast[4] = T1.C21[iVec]; - t[iVec]->CLast[5] = T1.C22[iVec]; - t[iVec]->CLast[6] = T1.C30[iVec]; - t[iVec]->CLast[7] = T1.C31[iVec]; - t[iVec]->CLast[8] = T1.C32[iVec]; - t[iVec]->CLast[9] = T1.C33[iVec]; - t[iVec]->CLast[10] = T1.C40[iVec]; - t[iVec]->CLast[11] = T1.C41[iVec]; - t[iVec]->CLast[12] = T1.C42[iVec]; - t[iVec]->CLast[13] = T1.C43[iVec]; - t[iVec]->CLast[14] = T1.C44[iVec]; - t[iVec]->CLast[15] = T1.C50[iVec]; - t[iVec]->CLast[16] = T1.C51[iVec]; - t[iVec]->CLast[17] = T1.C52[iVec]; - t[iVec]->CLast[18] = T1.C53[iVec]; - t[iVec]->CLast[19] = T1.C54[iVec]; - t[iVec]->CLast[20] = T1.C55[iVec]; - - - t[iVec]->chi2 = T1.chi2[iVec]; - t[iVec]->NDF = (int) T1.NDF[iVec]; + t[iVec]->TLast[0] = tr.x[iVec]; + t[iVec]->TLast[1] = tr.y[iVec]; + t[iVec]->TLast[2] = tr.tx[iVec]; + t[iVec]->TLast[3] = tr.ty[iVec]; + t[iVec]->TLast[4] = tr.qp[iVec]; + t[iVec]->TLast[5] = tr.z[iVec]; + t[iVec]->TLast[6] = tr.t[iVec]; + + t[iVec]->CLast[0] = tr.C00[iVec]; + t[iVec]->CLast[1] = tr.C10[iVec]; + t[iVec]->CLast[2] = tr.C11[iVec]; + t[iVec]->CLast[3] = tr.C20[iVec]; + t[iVec]->CLast[4] = tr.C21[iVec]; + t[iVec]->CLast[5] = tr.C22[iVec]; + t[iVec]->CLast[6] = tr.C30[iVec]; + t[iVec]->CLast[7] = tr.C31[iVec]; + t[iVec]->CLast[8] = tr.C32[iVec]; + t[iVec]->CLast[9] = tr.C33[iVec]; + t[iVec]->CLast[10] = tr.C40[iVec]; + t[iVec]->CLast[11] = tr.C41[iVec]; + t[iVec]->CLast[12] = tr.C42[iVec]; + t[iVec]->CLast[13] = tr.C43[iVec]; + t[iVec]->CLast[14] = tr.C44[iVec]; + t[iVec]->CLast[15] = tr.C50[iVec]; + t[iVec]->CLast[16] = tr.C51[iVec]; + t[iVec]->CLast[17] = tr.C52[iVec]; + t[iVec]->CLast[18] = tr.C53[iVec]; + t[iVec]->CLast[19] = tr.C54[iVec]; + t[iVec]->CLast[20] = tr.C55[iVec]; + + + t[iVec]->chi2 = tr.chi2[iVec]; + t[iVec]->NDF = (int) tr.NDF[iVec]; } if (iter == 1) continue; // only 1.5 iterations @@ -1173,15 +1174,15 @@ void L1Algo::L1KFTrackFitterMuch() FilterFirstL(T1, x_last, y_last, time_last, time_er_lst, staLast, d_xx_lst, d_yy_lst, d_xy_lst); qp0 = T.qp; - qp01 = T1.fqp; + qp01 = tr.qp; // fvec initialised = fvec(z[i] < z_end) & fvec(z_start <= z[i]); // fvec wIn = (ONE & (initialised)); - // T1.EnergyLossCorrectionAl( fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, fvec(1.f), wIn); + // T1.EnergyLossCorrectionAl( fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(1.f), wIn); // - // T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, ONE, sta[i].materialInfo.thick, 0); + // T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, ONE, sta[i].materialInfo.thick, 0); for (--i; i >= 0; i--) { @@ -1193,11 +1194,11 @@ void L1Algo::L1KFTrackFitterMuch() fvec z_last = z[i]; - // T1.ExtrapolateLine( T1.fz + 10, &w1); + // T1.ExtrapolateLine( tr.z + 10, &w1); // L1ExtrapolateLine( T, T.z + 10); - fvec d_z = T1.fz - z_last; + fvec d_z = tr.z - z_last; const fvec st = fvec(10); fvec nofSteps = (abs(d_z) / 10); @@ -1210,13 +1211,13 @@ void L1Algo::L1KFTrackFitterMuch() fvec nofSteps1 = fvec(0); fvec stepSize = wIn * iif((nofSteps < fvec::One()), d_z, st); - fvec z_cur = T1.fz; + fvec z_cur = tr.z; fvec w2 = wIn; for (int iStep = 0; iStep < max_steps + 1; iStep++) { const fmask maskLastStep = (nofSteps == nofSteps1); - z_cur = iif(maskLastStep, z_last, T1.fz - stepSize); + z_cur = iif(maskLastStep, z_last, tr.z - stepSize); // fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01)); // T1.ExtrapolateLine1( z_cur, &w2, v_mc); @@ -1227,16 +1228,16 @@ void L1Algo::L1KFTrackFitterMuch() // TODO: Unify the selection of energy loss correction! (S.Zharko) if constexpr (L1Constants::control::kIfUseRadLengthTable) { if (i == 11 || i == 14 || i == 17) - T1.EnergyLossCorrectionIron(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + 1), qp01, + T1.EnergyLossCorrectionIron(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01, fvec(1.f), w2); if (i == 8) - T1.EnergyLossCorrectionCarbon(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + 1), qp01, + T1.EnergyLossCorrectionCarbon(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01, fvec(1.f), w2); if (i == 9 || i == 10 || i == 12 || i == 13 || i == 15 || i == 16 || i == 18) - T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + 1), qp01, + T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01, fvec(1.f), w2); - T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy) / (nofSteps + fvec(1)), qp01, w2, + T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, w2, sta[i].materialInfo.thick / (nofSteps + fvec(1)), 0); w2 = iif(!maskLastStep, w2, fvec::Zero()); @@ -1274,20 +1275,20 @@ void L1Algo::L1KFTrackFitterMuch() T1.ExtrapolateLine(z[i + 1]); fldZ1 = z[i + 1]; - sta[i + 1].fieldSlice.GetFieldValue(T1.fx, T1.fy, fldB1); + sta[i + 1].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1); fldB1.Combine(fB[i + 1], w[i + 1]); fldZ2 = z[i - 1]; dz = fldZ2 - fldZ1; - sta[i + 1].fieldSlice.GetFieldValue(T1.fx + T1.ftx * dz, T1.fy + T1.fty * dz, fldB2); + sta[i + 1].fieldSlice.GetFieldValue(tr.x + tr.tx * dz, tr.y + tr.ty * dz, fldB2); fldB2.Combine(fB[i - 1], w[i - 1]); } fldZ0 = z[i]; dz = (fldZ1 - fldZ0); - sta[i].fieldSlice.GetFieldValue(T1.fx - T1.ftx * dz, T1.fy - T1.fty * dz, fldB0); + sta[i].fieldSlice.GetFieldValue(tr.x - tr.tx * dz, tr.y - tr.ty * dz, fldB0); fldB0.Combine(fB[i], w[i]); fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2); @@ -1299,9 +1300,9 @@ void L1Algo::L1KFTrackFitterMuch() if constexpr (L1Constants::control::kIfUseRadLengthTable) { - T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, fvec(1.f), wIn); - T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, wIn, - sta[i].materialInfo.thick, 0); + T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(1.f), wIn); + T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick, + 0); } else { fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn); @@ -1327,73 +1328,73 @@ void L1Algo::L1KFTrackFitterMuch() for (iVec = 0; iVec < nTracks_SIMD; iVec++) { - t[iVec]->TFirst[0] = T1.fx[iVec]; - t[iVec]->TFirst[1] = T1.fy[iVec]; - t[iVec]->TFirst[2] = T1.ftx[iVec]; - t[iVec]->TFirst[3] = T1.fty[iVec]; - t[iVec]->TFirst[4] = T1.fqp[iVec]; - t[iVec]->TFirst[5] = T1.fz[iVec]; - t[iVec]->TFirst[6] = T1.ft[iVec]; - - t[iVec]->CFirst[0] = T1.C00[iVec]; - t[iVec]->CFirst[1] = T1.C10[iVec]; - t[iVec]->CFirst[2] = T1.C11[iVec]; - t[iVec]->CFirst[3] = T1.C20[iVec]; - t[iVec]->CFirst[4] = T1.C21[iVec]; - t[iVec]->CFirst[5] = T1.C22[iVec]; - t[iVec]->CFirst[6] = T1.C30[iVec]; - t[iVec]->CFirst[7] = T1.C31[iVec]; - t[iVec]->CFirst[8] = T1.C32[iVec]; - t[iVec]->CFirst[9] = T1.C33[iVec]; - t[iVec]->CFirst[10] = T1.C40[iVec]; - t[iVec]->CFirst[11] = T1.C41[iVec]; - t[iVec]->CFirst[12] = T1.C42[iVec]; - t[iVec]->CFirst[13] = T1.C43[iVec]; - t[iVec]->CFirst[14] = T1.C44[iVec]; - t[iVec]->CFirst[15] = T1.C50[iVec]; - t[iVec]->CFirst[16] = T1.C51[iVec]; - t[iVec]->CFirst[17] = T1.C52[iVec]; - t[iVec]->CFirst[18] = T1.C53[iVec]; - t[iVec]->CFirst[19] = T1.C54[iVec]; - t[iVec]->CFirst[20] = T1.C55[iVec]; - - t[iVec]->chi2 = T1.chi2[iVec]; - t[iVec]->NDF = (int) T1.NDF[iVec]; + t[iVec]->TFirst[0] = tr.x[iVec]; + t[iVec]->TFirst[1] = tr.y[iVec]; + t[iVec]->TFirst[2] = tr.tx[iVec]; + t[iVec]->TFirst[3] = tr.ty[iVec]; + t[iVec]->TFirst[4] = tr.qp[iVec]; + t[iVec]->TFirst[5] = tr.z[iVec]; + t[iVec]->TFirst[6] = tr.t[iVec]; + + t[iVec]->CFirst[0] = tr.C00[iVec]; + t[iVec]->CFirst[1] = tr.C10[iVec]; + t[iVec]->CFirst[2] = tr.C11[iVec]; + t[iVec]->CFirst[3] = tr.C20[iVec]; + t[iVec]->CFirst[4] = tr.C21[iVec]; + t[iVec]->CFirst[5] = tr.C22[iVec]; + t[iVec]->CFirst[6] = tr.C30[iVec]; + t[iVec]->CFirst[7] = tr.C31[iVec]; + t[iVec]->CFirst[8] = tr.C32[iVec]; + t[iVec]->CFirst[9] = tr.C33[iVec]; + t[iVec]->CFirst[10] = tr.C40[iVec]; + t[iVec]->CFirst[11] = tr.C41[iVec]; + t[iVec]->CFirst[12] = tr.C42[iVec]; + t[iVec]->CFirst[13] = tr.C43[iVec]; + t[iVec]->CFirst[14] = tr.C44[iVec]; + t[iVec]->CFirst[15] = tr.C50[iVec]; + t[iVec]->CFirst[16] = tr.C51[iVec]; + t[iVec]->CFirst[17] = tr.C52[iVec]; + t[iVec]->CFirst[18] = tr.C53[iVec]; + t[iVec]->CFirst[19] = tr.C54[iVec]; + t[iVec]->CFirst[20] = tr.C55[iVec]; + + t[iVec]->chi2 = tr.chi2[iVec]; + t[iVec]->NDF = (int) tr.NDF[iVec]; } // extrapolate to the PV region L1TrackParFit T1pv = T1; - T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fqp, fld, fvec(1.)); + T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fTr.qp, fld, fvec(1.)); for (iVec = 0; iVec < nTracks_SIMD; iVec++) { - t[iVec]->Tpv[0] = T1pv.fx[iVec]; - t[iVec]->Tpv[1] = T1pv.fy[iVec]; - t[iVec]->Tpv[2] = T1pv.ftx[iVec]; - t[iVec]->Tpv[3] = T1pv.fty[iVec]; - t[iVec]->Tpv[4] = T1pv.fqp[iVec]; - t[iVec]->Tpv[5] = T1pv.fz[iVec]; - t[iVec]->Tpv[6] = T1pv.ft[iVec]; - - t[iVec]->Cpv[0] = T1pv.C00[iVec]; - t[iVec]->Cpv[1] = T1pv.C10[iVec]; - t[iVec]->Cpv[2] = T1pv.C11[iVec]; - t[iVec]->Cpv[3] = T1pv.C20[iVec]; - t[iVec]->Cpv[4] = T1pv.C21[iVec]; - t[iVec]->Cpv[5] = T1pv.C22[iVec]; - t[iVec]->Cpv[6] = T1pv.C30[iVec]; - t[iVec]->Cpv[7] = T1pv.C31[iVec]; - t[iVec]->Cpv[8] = T1pv.C32[iVec]; - t[iVec]->Cpv[9] = T1pv.C33[iVec]; - t[iVec]->Cpv[10] = T1pv.C40[iVec]; - t[iVec]->Cpv[11] = T1pv.C41[iVec]; - t[iVec]->Cpv[12] = T1pv.C42[iVec]; - t[iVec]->Cpv[13] = T1pv.C43[iVec]; - t[iVec]->Cpv[14] = T1pv.C44[iVec]; - t[iVec]->Cpv[15] = T1pv.C50[iVec]; - t[iVec]->Cpv[16] = T1pv.C51[iVec]; - t[iVec]->Cpv[17] = T1pv.C52[iVec]; - t[iVec]->Cpv[18] = T1pv.C53[iVec]; - t[iVec]->Cpv[19] = T1pv.C54[iVec]; - t[iVec]->Cpv[20] = T1pv.C55[iVec]; + t[iVec]->Tpv[0] = T1pv.fTr.x[iVec]; + t[iVec]->Tpv[1] = T1pv.fTr.y[iVec]; + t[iVec]->Tpv[2] = T1pv.fTr.tx[iVec]; + t[iVec]->Tpv[3] = T1pv.fTr.ty[iVec]; + t[iVec]->Tpv[4] = T1pv.fTr.qp[iVec]; + t[iVec]->Tpv[5] = T1pv.fTr.z[iVec]; + t[iVec]->Tpv[6] = T1pv.fTr.t[iVec]; + + t[iVec]->Cpv[0] = T1pv.fTr.C00[iVec]; + t[iVec]->Cpv[1] = T1pv.fTr.C10[iVec]; + t[iVec]->Cpv[2] = T1pv.fTr.C11[iVec]; + t[iVec]->Cpv[3] = T1pv.fTr.C20[iVec]; + t[iVec]->Cpv[4] = T1pv.fTr.C21[iVec]; + t[iVec]->Cpv[5] = T1pv.fTr.C22[iVec]; + t[iVec]->Cpv[6] = T1pv.fTr.C30[iVec]; + t[iVec]->Cpv[7] = T1pv.fTr.C31[iVec]; + t[iVec]->Cpv[8] = T1pv.fTr.C32[iVec]; + t[iVec]->Cpv[9] = T1pv.fTr.C33[iVec]; + t[iVec]->Cpv[10] = T1pv.fTr.C40[iVec]; + t[iVec]->Cpv[11] = T1pv.fTr.C41[iVec]; + t[iVec]->Cpv[12] = T1pv.fTr.C42[iVec]; + t[iVec]->Cpv[13] = T1pv.fTr.C43[iVec]; + t[iVec]->Cpv[14] = T1pv.fTr.C44[iVec]; + t[iVec]->Cpv[15] = T1pv.fTr.C50[iVec]; + t[iVec]->Cpv[16] = T1pv.fTr.C51[iVec]; + t[iVec]->Cpv[17] = T1pv.fTr.C52[iVec]; + t[iVec]->Cpv[18] = T1pv.fTr.C53[iVec]; + t[iVec]->Cpv[19] = T1pv.fTr.C54[iVec]; + t[iVec]->Cpv[20] = T1pv.fTr.C55[iVec]; } } // iter } @@ -1532,22 +1533,22 @@ void L1Algo::GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec Ai5 = (-A1 * A1 + A0 * A2); fvec L, L1; - t.fx = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det; - t.ftx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det; - fvec txtx1 = fvec(1.) + t.ftx * t.ftx; + t.fTr.x = (Ai0 * a0 + Ai1 * a1 + Ai3 * a2) * det; + t.fTr.tx = (Ai1 * a0 + Ai2 * a1 + Ai4 * a2) * det; + fvec txtx1 = fvec(1.) + t.fTr.tx * t.fTr.tx; L = (Ai3 * a0 + Ai4 * a1 + Ai5 * a2) * det / txtx1; - L1 = L * t.ftx; + L1 = L * t.fTr.tx; A1 = A1 + A3 * L1; A2 = A2 + (A4 + A4 + A5 * L1) * L1; b1 += b2 * L1; det = fvec::One() / (-A1 * A1 + A0 * A2); - t.fy = (A2 * b0 - A1 * b1) * det; - t.fty = (-A1 * b0 + A0 * b1) * det; - t.fqp = -L * c_light_i / sqrt(txtx1 + t.fty * t.fty); - if (timeV) { t.ft = iif(nhits > fvec::Zero(), time / nhits, fvec::Zero()); } + t.fTr.y = (A2 * b0 - A1 * b1) * det; + t.fTr.ty = (-A1 * b0 + A0 * b1) * det; + t.fTr.qp = -L * c_light_i / sqrt(txtx1 + t.fTr.ty * t.fTr.ty); + if (timeV) { t.fTr.t = iif(nhits > fvec::Zero(), time / nhits, fvec::Zero()); } - t.fz = z0; + t.fTr.z = z0; } void L1Algo::FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st) @@ -1583,95 +1584,95 @@ void L1Algo::FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st) void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st) { // initialize covariance matrix - track.C00 = st.XYInfo.C00; - track.C10 = st.XYInfo.C10; - track.C11 = st.XYInfo.C11; - track.C20 = ZERO; - track.C21 = ZERO; - track.C22 = vINF; - track.C30 = ZERO; - track.C31 = ZERO; - track.C32 = ZERO; - track.C33 = vINF; - track.C40 = ZERO; - track.C41 = ZERO; - track.C42 = ZERO; - track.C43 = ZERO; - track.C44 = ONE; - track.C50 = ZERO; - track.C51 = ZERO; - track.C52 = ZERO; - track.C53 = ZERO; - track.C54 = ZERO; - track.C55 = 2.6f * 2.6f; - - track.fx = x; - track.fy = y; - track.ft = t; - track.NDF = -3.0; - track.chi2 = ZERO; + track.fTr.C00 = st.XYInfo.C00; + track.fTr.C10 = st.XYInfo.C10; + track.fTr.C11 = st.XYInfo.C11; + track.fTr.C20 = ZERO; + track.fTr.C21 = ZERO; + track.fTr.C22 = vINF; + track.fTr.C30 = ZERO; + track.fTr.C31 = ZERO; + track.fTr.C32 = ZERO; + track.fTr.C33 = vINF; + track.fTr.C40 = ZERO; + track.fTr.C41 = ZERO; + track.fTr.C42 = ZERO; + track.fTr.C43 = ZERO; + track.fTr.C44 = ONE; + track.fTr.C50 = ZERO; + track.fTr.C51 = ZERO; + track.fTr.C52 = ZERO; + track.fTr.C53 = ZERO; + track.fTr.C54 = ZERO; + track.fTr.C55 = 2.6f * 2.6f; + + track.fTr.x = x; + track.fTr.y = y; + track.fTr.t = t; + track.fTr.NDF = -3.0; + track.fTr.chi2 = ZERO; } void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt, L1Station& st, fvec& /*d_xx*/, fvec& /*d_yy*/, fvec& /*d_xy*/) { - track.C00 = st.XYInfo.C00; - track.C10 = st.XYInfo.C10; - track.C11 = st.XYInfo.C11; - track.C20 = ZERO; - track.C21 = ZERO; - track.C22 = vINF; - track.C30 = ZERO; - track.C31 = ZERO; - track.C32 = ZERO; - track.C33 = vINF; - track.C40 = ZERO; - track.C41 = ZERO; - track.C42 = ZERO; - track.C43 = ZERO; - track.C44 = ONE; - track.C50 = ZERO; - track.C51 = ZERO; - track.C52 = ZERO; - track.C53 = ZERO; - track.C54 = ZERO; - track.C55 = dt * dt; - - track.fx = x; - track.fy = y; - track.ft = t; - track.NDF = -3.0; - track.chi2 = ZERO; + track.fTr.C00 = st.XYInfo.C00; + track.fTr.C10 = st.XYInfo.C10; + track.fTr.C11 = st.XYInfo.C11; + track.fTr.C20 = ZERO; + track.fTr.C21 = ZERO; + track.fTr.C22 = vINF; + track.fTr.C30 = ZERO; + track.fTr.C31 = ZERO; + track.fTr.C32 = ZERO; + track.fTr.C33 = vINF; + track.fTr.C40 = ZERO; + track.fTr.C41 = ZERO; + track.fTr.C42 = ZERO; + track.fTr.C43 = ZERO; + track.fTr.C44 = ONE; + track.fTr.C50 = ZERO; + track.fTr.C51 = ZERO; + track.fTr.C52 = ZERO; + track.fTr.C53 = ZERO; + track.fTr.C54 = ZERO; + track.fTr.C55 = dt * dt; + + track.fTr.x = x; + track.fTr.y = y; + track.fTr.t = t; + track.fTr.NDF = -3.0; + track.fTr.chi2 = ZERO; } void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt, L1Station& st) { - track.C00 = st.XYInfo.C00; - track.C10 = st.XYInfo.C10; - track.C11 = st.XYInfo.C11; - track.C20 = ZERO; - track.C21 = ZERO; - track.C22 = vINF; - track.C30 = ZERO; - track.C31 = ZERO; - track.C32 = ZERO; - track.C33 = vINF; - track.C40 = ZERO; - track.C41 = ZERO; - track.C42 = ZERO; - track.C43 = ZERO; - track.C44 = ONE; - track.C50 = ZERO; - track.C51 = ZERO; - track.C52 = ZERO; - track.C53 = ZERO; - track.C54 = ZERO; - track.C55 = dt * dt; - - track.fx = x; - track.fy = y; - track.NDF = -3.0; - track.chi2 = ZERO; + track.fTr.C00 = st.XYInfo.C00; + track.fTr.C10 = st.XYInfo.C10; + track.fTr.C11 = st.XYInfo.C11; + track.fTr.C20 = ZERO; + track.fTr.C21 = ZERO; + track.fTr.C22 = vINF; + track.fTr.C30 = ZERO; + track.fTr.C31 = ZERO; + track.fTr.C32 = ZERO; + track.fTr.C33 = vINF; + track.fTr.C40 = ZERO; + track.fTr.C41 = ZERO; + track.fTr.C42 = ZERO; + track.fTr.C43 = ZERO; + track.fTr.C44 = ONE; + track.fTr.C50 = ZERO; + track.fTr.C51 = ZERO; + track.fTr.C52 = ZERO; + track.fTr.C53 = ZERO; + track.fTr.C54 = ZERO; + track.fTr.C55 = dt * dt; + + track.fTr.x = x; + track.fTr.y = y; + track.fTr.NDF = -3.0; + track.fTr.chi2 = ZERO; } @@ -1679,34 +1680,34 @@ void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, f fvec& d_yy, fvec& d_xy) { // initialize covariance matrix - track.C00 = d_xx; - track.C10 = d_xy; - track.C11 = d_yy; - // track.C00= st.XYInfo.C00; - // track.C10= st.XYInfo.C10; - // track.C11= st.XYInfo.C11; - track.C20 = ZERO; - track.C21 = ZERO; - track.C22 = vINF; - track.C30 = ZERO; - track.C31 = ZERO; - track.C32 = ZERO; - track.C33 = vINF; - track.C40 = ZERO; - track.C41 = ZERO; - track.C42 = ZERO; - track.C43 = ZERO; - track.C44 = ONE; - track.C50 = ZERO; - track.C51 = ZERO; - track.C52 = ZERO; - track.C53 = ZERO; - track.C54 = ZERO; - track.C55 = dt * dt; - - track.fx = x; - track.fy = y; - // track.ft = t; - track.NDF = -3.0; // TODO: Why -3.0? (S.Zharko) - track.chi2 = ZERO; + track.fTr.C00 = d_xx; + track.fTr.C10 = d_xy; + track.fTr.C11 = d_yy; + // track.fTr.C00= st.XYInfo.C00; + // track.fTr.C10= st.XYInfo.C10; + // track.fTr.C11= st.XYInfo.C11; + track.fTr.C20 = ZERO; + track.fTr.C21 = ZERO; + track.fTr.C22 = vINF; + track.fTr.C30 = ZERO; + track.fTr.C31 = ZERO; + track.fTr.C32 = ZERO; + track.fTr.C33 = vINF; + track.fTr.C40 = ZERO; + track.fTr.C41 = ZERO; + track.fTr.C42 = ZERO; + track.fTr.C43 = ZERO; + track.fTr.C44 = ONE; + track.fTr.C50 = ZERO; + track.fTr.C51 = ZERO; + track.fTr.C52 = ZERO; + track.fTr.C53 = ZERO; + track.fTr.C54 = ZERO; + track.fTr.C55 = dt * dt; + + track.fTr.x = x; + track.fTr.y = y; + // track.fTr.t = t; + track.fTr.NDF = -3.0; // TODO: Why -3.0? (S.Zharko) + track.fTr.chi2 = ZERO; } diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx index cc9291a50f06d682567509e7ea433f345f40a2b1..7c40030a5bc364ef6459a2daceb574e7b38a46cd 100644 --- a/reco/L1/L1Algo/L1TrackParFit.cxx +++ b/reco/L1/L1Algo/L1TrackParFit.cxx @@ -15,18 +15,18 @@ void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w) fvec F0, F1, F2, F3, F4, F5; fvec K1, K2, K3, K4, K5; - zeta = info.cos_phi * fx + info.sin_phi * fy - u; + zeta = info.cos_phi * fTr.x + info.sin_phi * fTr.y - u; // F = CH' - F0 = info.cos_phi * C00 + info.sin_phi * C10; - F1 = info.cos_phi * C10 + info.sin_phi * C11; + F0 = info.cos_phi * fTr.C00 + info.sin_phi * fTr.C10; + F1 = info.cos_phi * fTr.C10 + info.sin_phi * fTr.C11; HCH = (F0 * info.cos_phi + F1 * info.sin_phi); - F2 = info.cos_phi * C20 + info.sin_phi * C21; - F3 = info.cos_phi * C30 + info.sin_phi * C31; - F4 = info.cos_phi * C40 + info.sin_phi * C41; - F5 = info.cos_phi * C50 + info.sin_phi * C51; + F2 = info.cos_phi * fTr.C20 + info.sin_phi * fTr.C21; + F3 = info.cos_phi * fTr.C30 + info.sin_phi * fTr.C31; + F4 = info.cos_phi * fTr.C40 + info.sin_phi * fTr.C41; + F5 = info.cos_phi * fTr.C50 + info.sin_phi * fTr.C51; const fmask maskDoFilter = (HCH < info.sigma2 * 16.f); //const fvec maskDoFilter = _f32vec4_true; @@ -37,9 +37,9 @@ void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w) fvec wi = w / (info.sigma2 + fvec(1.0000001) * HCH); fvec zetawi = w * zeta / (iif(maskDoFilter, info.sigma2, fvec::Zero()) + HCH); - // chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() ); - chi2 += zeta * zeta * wi; - NDF += w; + // fTr.chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() ); + fTr.chi2 += zeta * zeta * wi; + fTr.NDF += w; K1 = F1 * wi; K2 = F2 * wi; @@ -47,34 +47,34 @@ void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w) K4 = F4 * wi; K5 = F5 * wi; - fx -= F0 * zetawi; - fy -= F1 * zetawi; - ftx -= F2 * zetawi; - fty -= F3 * zetawi; - fqp -= F4 * zetawi; - ft -= F5 * zetawi; - - C00 -= F0 * F0 * wi; - C10 -= K1 * F0; - C11 -= K1 * F1; - C20 -= K2 * F0; - C21 -= K2 * F1; - C22 -= K2 * F2; - C30 -= K3 * F0; - C31 -= K3 * F1; - C32 -= K3 * F2; - C33 -= K3 * F3; - C40 -= K4 * F0; - C41 -= K4 * F1; - C42 -= K4 * F2; - C43 -= K4 * F3; - C44 -= K4 * F4; - C50 -= K5 * F0; - C51 -= K5 * F1; - C52 -= K5 * F2; - C53 -= K5 * F3; - C54 -= K5 * F4; - C55 -= K5 * F5; + fTr.x -= F0 * zetawi; + fTr.y -= F1 * zetawi; + fTr.tx -= F2 * zetawi; + fTr.ty -= F3 * zetawi; + fTr.qp -= F4 * zetawi; + fTr.t -= F5 * zetawi; + + fTr.C00 -= F0 * F0 * wi; + fTr.C10 -= K1 * F0; + fTr.C11 -= K1 * F1; + fTr.C20 -= K2 * F0; + fTr.C21 -= K2 * F1; + fTr.C22 -= K2 * F2; + fTr.C30 -= K3 * F0; + fTr.C31 -= K3 * F1; + fTr.C32 -= K3 * F2; + fTr.C33 -= K3 * F3; + fTr.C40 -= K4 * F0; + fTr.C41 -= K4 * F1; + fTr.C42 -= K4 * F2; + fTr.C43 -= K4 * F3; + fTr.C44 -= K4 * F4; + fTr.C50 -= K5 * F0; + fTr.C51 -= K5 * F1; + fTr.C52 -= K5 * F2; + fTr.C53 -= K5 * F3; + fTr.C54 -= K5 * F4; + fTr.C55 -= K5 * F5; } void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w) @@ -83,30 +83,30 @@ void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w) fvec F0, F1, F2, F3, F4, F5; fvec K1, K2, K3, K4, K5; - zeta = info.cos_phi * fx + info.sin_phi * fy - u; + zeta = info.cos_phi * fTr.x + info.sin_phi * fTr.y - u; // F = CH' - F0 = info.cos_phi * C00 + info.sin_phi * C10; - F1 = info.cos_phi * C10 + info.sin_phi * C11; + F0 = info.cos_phi * fTr.C00 + info.sin_phi * fTr.C10; + F1 = info.cos_phi * fTr.C10 + info.sin_phi * fTr.C11; HCH = (F0 * info.cos_phi + F1 * info.sin_phi); - F2 = info.cos_phi * C20 + info.sin_phi * C21; - F3 = info.cos_phi * C30 + info.sin_phi * C31; - F4 = info.cos_phi * C40 + info.sin_phi * C41; - F5 = info.cos_phi * C50 + info.sin_phi * C51; + F2 = info.cos_phi * fTr.C20 + info.sin_phi * fTr.C21; + F3 = info.cos_phi * fTr.C30 + info.sin_phi * fTr.C31; + F4 = info.cos_phi * fTr.C40 + info.sin_phi * fTr.C41; + F5 = info.cos_phi * fTr.C50 + info.sin_phi * fTr.C51; #if 0 // use mask const fmask mask = (HCH < info.sigma2 * 16.); wi = w/( (mask & info.sigma2) +HCH ); zetawi = zeta *wi; - chi2 += mask & (zeta * zetawi); + fTr.chi2 += mask & (zeta * zetawi); #else wi = w / (info.sigma2 + HCH); zetawi = zeta * wi; - chi2 += zeta * zetawi; + fTr.chi2 += zeta * zetawi; #endif // 0 - NDF += w; + fTr.NDF += w; K1 = F1 * wi; K2 = F2 * wi; @@ -114,34 +114,34 @@ void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w) K4 = F4 * wi; K5 = F5 * wi; - fx -= F0 * zetawi; - fy -= F1 * zetawi; - ftx -= F2 * zetawi; - fty -= F3 * zetawi; - // fqp -= F4*zetawi; - ft -= F5 * zetawi; - - C00 -= F0 * F0 * wi; - C10 -= K1 * F0; - C11 -= K1 * F1; - C20 -= K2 * F0; - C21 -= K2 * F1; - C22 -= K2 * F2; - C30 -= K3 * F0; - C31 -= K3 * F1; - C32 -= K3 * F2; - C33 -= K3 * F3; - // C40-= K4*F0; - // C41-= K4*F1; - // C42-= K4*F2; - // C43-= K4*F3; - // C44-= K4*F4; - C50 -= K5 * F0; - C51 -= K5 * F1; - C52 -= K5 * F2; - C53 -= K5 * F3; - C54 -= K5 * F4; - C55 -= K5 * F5; + fTr.x -= F0 * zetawi; + fTr.y -= F1 * zetawi; + fTr.tx -= F2 * zetawi; + fTr.ty -= F3 * zetawi; + // fTr.qp -= F4*zetawi; + fTr.t -= F5 * zetawi; + + fTr.C00 -= F0 * F0 * wi; + fTr.C10 -= K1 * F0; + fTr.C11 -= K1 * F1; + fTr.C20 -= K2 * F0; + fTr.C21 -= K2 * F1; + fTr.C22 -= K2 * F2; + fTr.C30 -= K3 * F0; + fTr.C31 -= K3 * F1; + fTr.C32 -= K3 * F2; + fTr.C33 -= K3 * F3; + // fTr.C40-= K4*F0; + // fTr.C41-= K4*F1; + // fTr.C42-= K4*F2; + // fTr.C43-= K4*F3; + // fTr.C44-= K4*F4; + fTr.C50 -= K5 * F0; + fTr.C51 -= K5 * F1; + fTr.C52 -= K5 * F2; + fTr.C53 -= K5 * F3; + fTr.C54 -= K5 * F4; + fTr.C55 -= K5 * F5; } void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo) @@ -149,14 +149,14 @@ void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo) // filter track with a time measurement // F = CH' - fvec F0 = C50; - fvec F1 = C51; - fvec F2 = C52; - fvec F3 = C53; - fvec F4 = C54; - fvec F5 = C55; + fvec F0 = fTr.C50; + fvec F1 = fTr.C51; + fvec F2 = fTr.C52; + fvec F3 = fTr.C53; + fvec F4 = fTr.C54; + fvec F5 = fTr.C55; - fvec HCH = C55; + fvec HCH = fTr.C55; w.setZero(timeInfo <= fvec::Zero()); @@ -170,11 +170,11 @@ void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo) const fmask maskDoFilter = (HCH < dt2 * 16.f); fvec wi = w / (dt2 + fvec(1.0000001) * HCH); - fvec zeta = ft - t; + fvec zeta = fTr.t - t; fvec zetawi = w * zeta / (iif(maskDoFilter, dt2, fvec::Zero()) + HCH); - chi2(maskDoFilter) += zeta * zeta * wi; - NDF += w; + fTr.chi2(maskDoFilter) += zeta * zeta * wi; + fTr.NDF += w; fvec K1 = F1 * wi; fvec K2 = F2 * wi; @@ -182,34 +182,34 @@ void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo) fvec K4 = F4 * wi; fvec K5 = F5 * wi; - fx -= F0 * zetawi; - fy -= F1 * zetawi; - ftx -= F2 * zetawi; - fty -= F3 * zetawi; - fqp -= F4 * zetawi; - ft -= F5 * zetawi; - - C00 -= F0 * F0 * wi; - C10 -= K1 * F0; - C11 -= K1 * F1; - C20 -= K2 * F0; - C21 -= K2 * F1; - C22 -= K2 * F2; - C30 -= K3 * F0; - C31 -= K3 * F1; - C32 -= K3 * F2; - C33 -= K3 * F3; - C40 -= K4 * F0; - C41 -= K4 * F1; - C42 -= K4 * F2; - C43 -= K4 * F3; - C44 -= K4 * F4; - C50 -= K5 * F0; - C51 -= K5 * F1; - C52 -= K5 * F2; - C53 -= K5 * F3; - C54 -= K5 * F4; - C55 -= K5 * F5; + fTr.x -= F0 * zetawi; + fTr.y -= F1 * zetawi; + fTr.tx -= F2 * zetawi; + fTr.ty -= F3 * zetawi; + fTr.qp -= F4 * zetawi; + fTr.t -= F5 * zetawi; + + fTr.C00 -= F0 * F0 * wi; + fTr.C10 -= K1 * F0; + fTr.C11 -= K1 * F1; + fTr.C20 -= K2 * F0; + fTr.C21 -= K2 * F1; + fTr.C22 -= K2 * F2; + fTr.C30 -= K3 * F0; + fTr.C31 -= K3 * F1; + fTr.C32 -= K3 * F2; + fTr.C33 -= K3 * F3; + fTr.C40 -= K4 * F0; + fTr.C41 -= K4 * F1; + fTr.C42 -= K4 * F2; + fTr.C43 -= K4 * F3; + fTr.C44 -= K4 * F4; + fTr.C50 -= K5 * F0; + fTr.C51 -= K5 * F1; + fTr.C52 -= K5 * F2; + fTr.C53 -= K5 * F3; + fTr.C54 -= K5 * F4; + fTr.C55 -= K5 * F5; } @@ -221,22 +221,22 @@ void L1TrackParFit::FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y) fvec F00, F10, F20, F30, F40, F50, F01, F11, F21, F31, F41, F51; fvec K00, K10, K20, K30, K40, K50, K01, K11, K21, K31, K41, K51; - zeta0 = fx - x; - zeta1 = fy - y; + zeta0 = fTr.x - x; + zeta1 = fTr.y - y; // F = CH' - F00 = C00; - F10 = C10; - F20 = C20; - F30 = C30; - F40 = C40; - F50 = C50; - F01 = C10; - F11 = C11; - F21 = C21; - F31 = C31; - F41 = C41; - F51 = C51; + F00 = fTr.C00; + F10 = fTr.C10; + F20 = fTr.C20; + F30 = fTr.C30; + F40 = fTr.C40; + F50 = fTr.C50; + F01 = fTr.C10; + F11 = fTr.C11; + F21 = fTr.C21; + F31 = fTr.C31; + F41 = fTr.C41; + F51 = fTr.C51; S00 = F00 + info.C00; S10 = F10 + info.C10; @@ -248,8 +248,8 @@ void L1TrackParFit::FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y) S10 = -si * S10; S11 = si * S00tmp; - chi2 += zeta0 * zeta0 * S00 + 2.f * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11; - NDF += TWO; + fTr.chi2 += zeta0 * zeta0 * S00 + 2.f * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11; + fTr.NDF += TWO; K00 = F00 * S00 + F01 * S10; K01 = F00 * S10 + F01 * S11; @@ -264,33 +264,33 @@ void L1TrackParFit::FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y) K50 = F50 * S00 + F51 * S10; K51 = F50 * S10 + F51 * S11; - fx -= K00 * zeta0 + K01 * zeta1; - fy -= K10 * zeta0 + K11 * zeta1; - ftx -= K20 * zeta0 + K21 * zeta1; - fty -= K30 * zeta0 + K31 * zeta1; - fqp -= K40 * zeta0 + K41 * zeta1; - - C00 -= K00 * F00 + K01 * F01; - C10 -= K10 * F00 + K11 * F01; - C11 -= K10 * F10 + K11 * F11; - C20 -= K20 * F00 + K21 * F01; - C21 -= K20 * F10 + K21 * F11; - C22 -= K20 * F20 + K21 * F21; - C30 -= K30 * F00 + K31 * F01; - C31 -= K30 * F10 + K31 * F11; - C32 -= K30 * F20 + K31 * F21; - C33 -= K30 * F30 + K31 * F31; - C40 -= K40 * F00 + K41 * F01; - C41 -= K40 * F10 + K41 * F11; - C42 -= K40 * F20 + K41 * F21; - C43 -= K40 * F30 + K41 * F31; - C44 -= K40 * F40 + K41 * F41; - C50 -= K50 * F00 + K51 * F01; - C51 -= K50 * F10 + K51 * F11; - C52 -= K50 * F20 + K51 * F21; - C53 -= K50 * F30 + K51 * F31; - C54 -= K50 * F40 + K51 * F41; - C55 -= K50 * F50 + K51 * F51; + fTr.x -= K00 * zeta0 + K01 * zeta1; + fTr.y -= K10 * zeta0 + K11 * zeta1; + fTr.tx -= K20 * zeta0 + K21 * zeta1; + fTr.ty -= K30 * zeta0 + K31 * zeta1; + fTr.qp -= K40 * zeta0 + K41 * zeta1; + + fTr.C00 -= K00 * F00 + K01 * F01; + fTr.C10 -= K10 * F00 + K11 * F01; + fTr.C11 -= K10 * F10 + K11 * F11; + fTr.C20 -= K20 * F00 + K21 * F01; + fTr.C21 -= K20 * F10 + K21 * F11; + fTr.C22 -= K20 * F20 + K21 * F21; + fTr.C30 -= K30 * F00 + K31 * F01; + fTr.C31 -= K30 * F10 + K31 * F11; + fTr.C32 -= K30 * F20 + K31 * F21; + fTr.C33 -= K30 * F30 + K31 * F31; + fTr.C40 -= K40 * F00 + K41 * F01; + fTr.C41 -= K40 * F10 + K41 * F11; + fTr.C42 -= K40 * F20 + K41 * F21; + fTr.C43 -= K40 * F30 + K41 * F31; + fTr.C44 -= K40 * F40 + K41 * F41; + fTr.C50 -= K50 * F00 + K51 * F01; + fTr.C51 -= K50 * F10 + K51 * F11; + fTr.C52 -= K50 * F20 + K51 * F21; + fTr.C53 -= K50 * F30 + K51 * F31; + fTr.C54 -= K50 * F40 + K51 * F41; + fTr.C55 -= K50 * F50 + K51 * F51; } @@ -299,45 +299,45 @@ void L1TrackParFit::ExtrapolateLine(fvec z_out, fvec* w) cnst c_light = 29.9792458; - fvec dz = (z_out - fz); + fvec dz = (z_out - fTr.z); if (w) { dz.setZero(*w <= fvec::Zero()); } - fx += dz * ftx; - fy += dz * fty; - fz += dz; - ft += dz * sqrt(fvec(1.) + ftx * ftx + fty * fty) / c_light; + fTr.x += dz * fTr.tx; + fTr.y += dz * fTr.ty; + fTr.z += dz; + fTr.t += dz * sqrt(fvec(1.) + fTr.tx * fTr.tx + fTr.ty * fTr.ty) / c_light; - const fvec k1 = dz * ftx / (c_light * sqrt(ftx * ftx + fty * fty + fvec(1.))); - const fvec k2 = dz * fty / (c_light * sqrt(ftx * ftx + fty * fty + fvec(1.))); + const fvec k1 = dz * fTr.tx / (c_light * sqrt(fTr.tx * fTr.tx + fTr.ty * fTr.ty + fvec(1.))); + const fvec k2 = dz * fTr.ty / (c_light * sqrt(fTr.tx * fTr.tx + fTr.ty * fTr.ty + fvec(1.))); - const fvec dzC32_in = dz * C32; + const fvec dzC32_in = dz * fTr.C32; - C21 += dzC32_in; - C10 += dz * (C21 + C30); + fTr.C21 += dzC32_in; + fTr.C10 += dz * (fTr.C21 + fTr.C30); - const fvec C20_in = C20; + const fvec C20_in = fTr.C20; - C20 += dz * C22; - C00 += dz * (C20 + C20_in); + fTr.C20 += dz * fTr.C22; + fTr.C00 += dz * (fTr.C20 + C20_in); - const fvec C31_in = C31; + const fvec C31_in = fTr.C31; - C31 += dz * C33; - C11 += dz * (C31 + C31_in); - C30 += dzC32_in; + fTr.C31 += dz * fTr.C33; + fTr.C11 += dz * (fTr.C31 + C31_in); + fTr.C30 += dzC32_in; - C40 += dz * C42; - C41 += dz * C43; + fTr.C40 += dz * fTr.C42; + fTr.C41 += dz * fTr.C43; - fvec c52 = C52; - fvec c53 = C53; + fvec c52 = fTr.C52; + fvec c53 = fTr.C53; - C50 += k1 * C20 + k2 * C30; - C51 += k1 * C21 + k2 * C31; - C52 += k1 * C22 + k2 * C32; - C53 += k1 * C32 + k2 * C33; - C54 += k1 * C42 + k2 * C43; - C55 += k1 * (C52 + c52) + k2 * (C53 + c53); + fTr.C50 += k1 * fTr.C20 + k2 * fTr.C30; + fTr.C51 += k1 * fTr.C21 + k2 * fTr.C31; + fTr.C52 += k1 * fTr.C22 + k2 * fTr.C32; + fTr.C53 += k1 * fTr.C32 + k2 * fTr.C33; + fTr.C54 += k1 * fTr.C42 + k2 * fTr.C43; + fTr.C55 += k1 * (fTr.C52 + c52) + k2 * (fTr.C53 + c53); } void L1TrackParFit::ExtrapolateLine1(fvec z_out, fvec* w, fvec v) @@ -345,46 +345,46 @@ void L1TrackParFit::ExtrapolateLine1(fvec z_out, fvec* w, fvec v) cnst c_light(29.9792458); - fvec dz = (z_out - fz); + fvec dz = (z_out - fTr.z); if (w) { dz.setZero(*w <= fvec::Zero()); } - fx += dz * ftx; - fy += dz * fty; - fz += dz; + fTr.x += dz * fTr.tx; + fTr.y += dz * fTr.ty; + fTr.z += dz; - ft += dz * sqrt(fvec(1.) + ftx * ftx + fty * fty) / (v * c_light); + fTr.t += dz * sqrt(fvec(1.) + fTr.tx * fTr.tx + fTr.ty * fTr.ty) / (v * c_light); - const fvec k1 = dz * ftx / ((v * c_light) * sqrt(ftx * ftx + fty * fty + fvec(1.))); - const fvec k2 = dz * fty / ((v * c_light) * sqrt(ftx * ftx + fty * fty + fvec(1.))); + const fvec k1 = dz * fTr.tx / ((v * c_light) * sqrt(fTr.tx * fTr.tx + fTr.ty * fTr.ty + fvec(1.))); + const fvec k2 = dz * fTr.ty / ((v * c_light) * sqrt(fTr.tx * fTr.tx + fTr.ty * fTr.ty + fvec(1.))); - const fvec dzC32_in = dz * C32; + const fvec dzC32_in = dz * fTr.C32; - C21 += dzC32_in; - C10 += dz * (C21 + C30); + fTr.C21 += dzC32_in; + fTr.C10 += dz * (fTr.C21 + fTr.C30); - const fvec C20_in = C20; + const fvec C20_in = fTr.C20; - C20 += dz * C22; - C00 += dz * (C20 + C20_in); + fTr.C20 += dz * fTr.C22; + fTr.C00 += dz * (fTr.C20 + C20_in); - const fvec C31_in = C31; + const fvec C31_in = fTr.C31; - C31 += dz * C33; - C11 += dz * (C31 + C31_in); - C30 += dzC32_in; + fTr.C31 += dz * fTr.C33; + fTr.C11 += dz * (fTr.C31 + C31_in); + fTr.C30 += dzC32_in; - C40 += dz * C42; - C41 += dz * C43; + fTr.C40 += dz * fTr.C42; + fTr.C41 += dz * fTr.C43; - fvec c52 = C52; - fvec c53 = C53; + fvec c52 = fTr.C52; + fvec c53 = fTr.C53; - C50 += k1 * C20 + k2 * C30; - C51 += k1 * C21 + k2 * C31; - C52 += k1 * C22 + k2 * C32; - C53 += k1 * C32 + k2 * C33; - C54 += k1 * C42 + k2 * C43; - C55 += k1 * (C52 + c52) + k2 * (C53 + c53); + fTr.C50 += k1 * fTr.C20 + k2 * fTr.C30; + fTr.C51 += k1 * fTr.C21 + k2 * fTr.C31; + fTr.C52 += k1 * fTr.C22 + k2 * fTr.C32; + fTr.C53 += k1 * fTr.C32 + k2 * fTr.C33; + fTr.C54 += k1 * fTr.C42 + k2 * fTr.C43; + fTr.C55 += k1 * (fTr.C52 + c52) + k2 * (fTr.C53 + c53); } void L1TrackParFit::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix @@ -392,9 +392,9 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja fvec qp0, // use Q/p linearisation at this value const L1FieldRegion& F, const fvec& w) { - fvec sgn = iif(fz < z_out, fvec(1.), fvec(-1.)); - while (!(w * abs(z_out - fz) <= fvec(1.e-6)).isFull()) { - fvec zNew = fz + sgn * fvec(50.); // max. 50 cm step + fvec sgn = iif(fTr.z < z_out, fvec(1.), fvec(-1.)); + while (!(w * abs(z_out - fTr.z) <= fvec(1.e-6)).isFull()) { + fvec zNew = fTr.z + sgn * fvec(50.); // max. 50 cm step zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out; ExtrapolateStep(zNew, qp0, F, w); } @@ -413,10 +413,10 @@ void // // | x | tx // | y | ty - // d/dz = | tx| = ft * ( ty * ( B(3)+tx*b(1) ) - (1+tx**2)*B(2) ) - // | ty| ft * (-tx * ( B(3)+ty+b(2) - (1+ty**2)*B(1) ) , + // d/dz = | tx| = fTr.t * ( ty * ( B(3)+tx*b(1) ) - (1+tx**2)*B(2) ) + // | ty| fTr.t * (-tx * ( B(3)+ty+b(2) - (1+ty**2)*B(1) ) , // - // where ft = c_light*qp*sqrt ( 1 + tx**2 + ty**2 ) . + // where fTr.t = c_light*qp*sqrt ( 1 + tx**2 + ty**2 ) . // // In the following for RK step : // @@ -442,23 +442,23 @@ void //---------------------------------------------------------------- - z_out = iif((fvec(0.f) < w), z_out, fz); + z_out = iif((fvec(0.f) < w), z_out, fTr.z); - fvec qp_in = fqp; - const fvec h = (z_out - fz); + fvec qp_in = fTr.qp; + const fvec h = (z_out - fTr.z); // cout<<h<<" h"<<endl; - // cout<<ftx<<" ftx"<<endl; - // cout<<fty<<" fty"<<endl; + // cout<<fTr.tx<<" fTr.tx"<<endl; + // cout<<fTr.ty<<" fTr.ty"<<endl; fvec hC = h * c_light; - //std::cout << "fx " << fx << std::endl; + //std::cout << "fTr.x " << fTr.x << std::endl; fvec x0[5]; - x0[0] = fx; - x0[1] = fy; - x0[2] = ftx; - x0[3] = fty; - x0[4] = ft; + x0[0] = fTr.x; + x0[1] = fTr.y; + x0[2] = fTr.tx; + x0[3] = fTr.ty; + x0[4] = fTr.t; // // Runge-Kutta step @@ -475,7 +475,7 @@ void fvec B[3]; - F.Get(x[0], x[1], fz + a[step] * h, B); + F.Get(x[0], x[1], fTr.z + a[step] * h, B); //std::cout << "extrapolation step " << step << " z " << z_in + a[step] * h << " field " << B[0] << " " << B[1] << " " // << B[2] << std::endl; @@ -526,21 +526,21 @@ void // cout << "w = " << *w << "; "; // cout << "initialised = " << initialised << "; "; - // cout << "fx = " << fx; + // cout << "fTr.x = " << fTr.x; //std::cout << " x : x0[0] " << x0[0] << " k[0] " << k[0] << " k[5] " << k[5] << " k[10] " << k[10] << " k[15] " // << k[15] << std::endl; - fx = x0[0] + c[0] * k[0] + c[1] * k[5 + 0] + c[2] * k[10 + 0] + c[3] * k[15 + 0]; - fy = x0[1] + c[0] * k[1] + c[1] * k[5 + 1] + c[2] * k[10 + 1] + c[3] * k[15 + 1]; - ftx = x0[2] + c[0] * k[2] + c[1] * k[5 + 2] + c[2] * k[10 + 2] + c[3] * k[15 + 2]; - fty = x0[3] + c[0] * k[3] + c[1] * k[5 + 3] + c[2] * k[10 + 3] + c[3] * k[15 + 3]; - ft = x0[4] + c[0] * k[4] + c[1] * k[5 + 4] + c[2] * k[10 + 4] + c[3] * k[15 + 4]; - fz = z_out; + fTr.x = x0[0] + c[0] * k[0] + c[1] * k[5 + 0] + c[2] * k[10 + 0] + c[3] * k[15 + 0]; + fTr.y = x0[1] + c[0] * k[1] + c[1] * k[5 + 1] + c[2] * k[10 + 1] + c[3] * k[15 + 1]; + fTr.tx = x0[2] + c[0] * k[2] + c[1] * k[5 + 2] + c[2] * k[10 + 2] + c[3] * k[15 + 2]; + fTr.ty = x0[3] + c[0] * k[3] + c[1] * k[5 + 3] + c[2] * k[10 + 3] + c[3] * k[15 + 3]; + fTr.t = x0[4] + c[0] * k[4] + c[1] * k[5 + 4] + c[2] * k[10 + 4] + c[3] * k[15 + 4]; + fTr.z = z_out; - // cout << "; fx = " << fx << endl; + // cout << "; fTr.x = " << fTr.x << endl; } - // cout<<fx<<" fx"<<endl; + // cout<<fTr.x<<" fTr.x"<<endl; fvec J[36]; // Jacobian of extrapolation @@ -688,16 +688,16 @@ void fvec dqp = qp_in - qp0; { // update parameters - fx += J[6 * 4 + 0] * dqp; - fy += J[6 * 4 + 1] * dqp; - ftx += J[6 * 4 + 2] * dqp; - fty += J[6 * 4 + 3] * dqp; - ft += J[6 * 4 + 5] * dqp; + fTr.x += J[6 * 4 + 0] * dqp; + fTr.y += J[6 * 4 + 1] * dqp; + fTr.tx += J[6 * 4 + 2] * dqp; + fTr.ty += J[6 * 4 + 3] * dqp; + fTr.t += J[6 * 4 + 5] * dqp; } - // cout<<fx<<" fx"<<endl; + // cout<<fTr.x<<" fTr.x"<<endl; // covariance matrix transport - //cout<< (ft - ft_old)<<" ft dt "<<endl; + //cout<< (fTr.t - fTr.t_old)<<" fTr.t dt "<<endl; // // if ( C_in&&C_out ) CbmKFMath::multQtSQ( 5, J, C_in, C_out); // TODO // j(0,2) = J[6*2 + 0]; @@ -721,75 +721,75 @@ void // j(4,4) = J[6*4 + 4]; // j(5,4) = J[6*4 + 5]; - const fvec c42 = C42, c43 = C43; - - const fvec cj00 = C00 + C20 * J[6 * 2 + 0] + C30 * J[6 * 3 + 0] + C40 * J[6 * 4 + 0]; - const fvec cj10 = C10 + C21 * J[6 * 2 + 0] + C31 * J[6 * 3 + 0] + C41 * J[6 * 4 + 0]; - const fvec cj20 = C20 + C22 * J[6 * 2 + 0] + C32 * J[6 * 3 + 0] + c42 * J[6 * 4 + 0]; - const fvec cj30 = C30 + C32 * J[6 * 2 + 0] + C33 * J[6 * 3 + 0] + c43 * J[6 * 4 + 0]; - const fvec cj40 = C40 + C42 * J[6 * 2 + 0] + C43 * J[6 * 3 + 0] + C44 * J[6 * 4 + 0]; - const fvec cj50 = C50 + C52 * J[6 * 2 + 0] + C53 * J[6 * 3 + 0] + C54 * J[6 * 4 + 0]; - - // const fvec cj01 = C10 + C20*J[6*2 + 1] + C30*J[6*3 + 1] + C40*J[6*4 + 1]; - const fvec cj11 = C11 + C21 * J[6 * 2 + 1] + C31 * J[6 * 3 + 1] + C41 * J[6 * 4 + 1]; - const fvec cj21 = C21 + C22 * J[6 * 2 + 1] + C32 * J[6 * 3 + 1] + c42 * J[6 * 4 + 1]; - const fvec cj31 = C31 + C32 * J[6 * 2 + 1] + C33 * J[6 * 3 + 1] + c43 * J[6 * 4 + 1]; - const fvec cj41 = C41 + C42 * J[6 * 2 + 1] + C43 * J[6 * 3 + 1] + C44 * J[6 * 4 + 1]; - const fvec cj51 = C51 + C52 * J[6 * 2 + 1] + C53 * J[6 * 3 + 1] + C54 * J[6 * 4 + 1]; - - // const fvec cj02 = C20*J[6*2 + 2] + C30*J[6*3 + 2] + C40*J[6*4 + 2]; - // const fvec cj12 = C21*J[6*2 + 2] + C31*J[6*3 + 2] + C41*J[6*4 + 2]; - const fvec cj22 = C22 * J[6 * 2 + 2] + C32 * J[6 * 3 + 2] + c42 * J[6 * 4 + 2]; - const fvec cj32 = C32 * J[6 * 2 + 2] + C33 * J[6 * 3 + 2] + c43 * J[6 * 4 + 2]; - const fvec cj42 = C42 * J[6 * 2 + 2] + C43 * J[6 * 3 + 2] + C44 * J[6 * 4 + 2]; - const fvec cj52 = C52 * J[6 * 2 + 2] + C53 * J[6 * 3 + 2] + C54 * J[6 * 4 + 2]; - - // const fvec cj03 = C20*J[6*2 + 3] + C30*J[6*3 + 3] + C40*J[6*4 + 3]; - // const fvec cj13 = C21*J[6*2 + 3] + C31*J[6*3 + 3] + C41*J[6*4 + 3]; - const fvec cj23 = C22 * J[6 * 2 + 3] + C32 * J[6 * 3 + 3] + c42 * J[6 * 4 + 3]; - const fvec cj33 = C32 * J[6 * 2 + 3] + C33 * J[6 * 3 + 3] + c43 * J[6 * 4 + 3]; - const fvec cj43 = C42 * J[6 * 2 + 3] + C43 * J[6 * 3 + 3] + C44 * J[6 * 4 + 3]; - const fvec cj53 = C52 * J[6 * 2 + 3] + C53 * J[6 * 3 + 3] + C54 * J[6 * 4 + 3]; - - const fvec cj24 = C42; - const fvec cj34 = C43; - const fvec cj44 = C44; - const fvec cj54 = C54; - - // const fvec cj05 = C50 + C20*J[17] + C30*J[23] + C40*J[29]; - // const fvec cj15 = C51 + C21*J[17] + C31*J[23] + C41*J[29]; - const fvec cj25 = C52 + C22 * J[17] + C32 * J[23] + C42 * J[29]; - const fvec cj35 = C53 + C32 * J[17] + C33 * J[23] + C43 * J[29]; - const fvec cj45 = C54 + C42 * J[17] + C43 * J[23] + C44 * J[29]; - const fvec cj55 = C55 + C52 * J[17] + C53 * J[23] + C54 * J[29]; - - - C00 = cj00 + cj20 * J[12] + cj30 * J[18] + cj40 * J[24]; - - C10 = cj10 + cj20 * J[13] + cj30 * J[19] + cj40 * J[25]; - C11 = cj11 + cj21 * J[13] + cj31 * J[19] + cj41 * J[25]; - - C20 = cj20 + cj30 * J[20] + cj40 * J[26]; - C21 = cj21 + cj31 * J[20] + cj41 * J[26]; - C22 = cj22 + cj32 * J[20] + cj42 * J[26]; - - C30 = cj30 + cj20 * J[15] + cj40 * J[27]; - C31 = cj31 + cj21 * J[15] + cj41 * J[27]; - C32 = cj32 + cj22 * J[15] + cj42 * J[27]; - C33 = cj33 + cj23 * J[15] + cj43 * J[27]; - - C40 = cj40; - C41 = cj41; - C42 = cj42; - C43 = cj43; - C44 = cj44; - - C50 = cj50 + cj20 * J[17] + cj30 * J[23] + cj40 * J[29]; - C51 = cj51 + cj21 * J[17] + cj31 * J[23] + cj41 * J[29]; - C52 = cj52 + cj22 * J[17] + cj32 * J[23] + cj42 * J[29]; - C53 = cj53 + cj23 * J[17] + cj33 * J[23] + cj43 * J[29]; - C54 = cj54 + cj24 * J[17] + cj34 * J[23] + cj44 * J[29]; - C55 = cj55 + cj25 * J[17] + cj35 * J[23] + cj45 * J[29]; + const fvec c42 = fTr.C42, c43 = fTr.C43; + + const fvec cj00 = fTr.C00 + fTr.C20 * J[6 * 2 + 0] + fTr.C30 * J[6 * 3 + 0] + fTr.C40 * J[6 * 4 + 0]; + const fvec cj10 = fTr.C10 + fTr.C21 * J[6 * 2 + 0] + fTr.C31 * J[6 * 3 + 0] + fTr.C41 * J[6 * 4 + 0]; + const fvec cj20 = fTr.C20 + fTr.C22 * J[6 * 2 + 0] + fTr.C32 * J[6 * 3 + 0] + c42 * J[6 * 4 + 0]; + const fvec cj30 = fTr.C30 + fTr.C32 * J[6 * 2 + 0] + fTr.C33 * J[6 * 3 + 0] + c43 * J[6 * 4 + 0]; + const fvec cj40 = fTr.C40 + fTr.C42 * J[6 * 2 + 0] + fTr.C43 * J[6 * 3 + 0] + fTr.C44 * J[6 * 4 + 0]; + const fvec cj50 = fTr.C50 + fTr.C52 * J[6 * 2 + 0] + fTr.C53 * J[6 * 3 + 0] + fTr.C54 * J[6 * 4 + 0]; + + // const fvec cj01 = fTr.C10 + fTr.C20*J[6*2 + 1] + fTr.C30*J[6*3 + 1] + fTr.C40*J[6*4 + 1]; + const fvec cj11 = fTr.C11 + fTr.C21 * J[6 * 2 + 1] + fTr.C31 * J[6 * 3 + 1] + fTr.C41 * J[6 * 4 + 1]; + const fvec cj21 = fTr.C21 + fTr.C22 * J[6 * 2 + 1] + fTr.C32 * J[6 * 3 + 1] + c42 * J[6 * 4 + 1]; + const fvec cj31 = fTr.C31 + fTr.C32 * J[6 * 2 + 1] + fTr.C33 * J[6 * 3 + 1] + c43 * J[6 * 4 + 1]; + const fvec cj41 = fTr.C41 + fTr.C42 * J[6 * 2 + 1] + fTr.C43 * J[6 * 3 + 1] + fTr.C44 * J[6 * 4 + 1]; + const fvec cj51 = fTr.C51 + fTr.C52 * J[6 * 2 + 1] + fTr.C53 * J[6 * 3 + 1] + fTr.C54 * J[6 * 4 + 1]; + + // const fvec cj02 = fTr.C20*J[6*2 + 2] + fTr.C30*J[6*3 + 2] + fTr.C40*J[6*4 + 2]; + // const fvec cj12 = fTr.C21*J[6*2 + 2] + fTr.C31*J[6*3 + 2] + fTr.C41*J[6*4 + 2]; + const fvec cj22 = fTr.C22 * J[6 * 2 + 2] + fTr.C32 * J[6 * 3 + 2] + c42 * J[6 * 4 + 2]; + const fvec cj32 = fTr.C32 * J[6 * 2 + 2] + fTr.C33 * J[6 * 3 + 2] + c43 * J[6 * 4 + 2]; + const fvec cj42 = fTr.C42 * J[6 * 2 + 2] + fTr.C43 * J[6 * 3 + 2] + fTr.C44 * J[6 * 4 + 2]; + const fvec cj52 = fTr.C52 * J[6 * 2 + 2] + fTr.C53 * J[6 * 3 + 2] + fTr.C54 * J[6 * 4 + 2]; + + // const fvec cj03 = fTr.C20*J[6*2 + 3] + fTr.C30*J[6*3 + 3] + fTr.C40*J[6*4 + 3]; + // const fvec cj13 = fTr.C21*J[6*2 + 3] + fTr.C31*J[6*3 + 3] + fTr.C41*J[6*4 + 3]; + const fvec cj23 = fTr.C22 * J[6 * 2 + 3] + fTr.C32 * J[6 * 3 + 3] + c42 * J[6 * 4 + 3]; + const fvec cj33 = fTr.C32 * J[6 * 2 + 3] + fTr.C33 * J[6 * 3 + 3] + c43 * J[6 * 4 + 3]; + const fvec cj43 = fTr.C42 * J[6 * 2 + 3] + fTr.C43 * J[6 * 3 + 3] + fTr.C44 * J[6 * 4 + 3]; + const fvec cj53 = fTr.C52 * J[6 * 2 + 3] + fTr.C53 * J[6 * 3 + 3] + fTr.C54 * J[6 * 4 + 3]; + + const fvec cj24 = fTr.C42; + const fvec cj34 = fTr.C43; + const fvec cj44 = fTr.C44; + const fvec cj54 = fTr.C54; + + // const fvec cj05 = fTr.C50 + fTr.C20*J[17] + fTr.C30*J[23] + fTr.C40*J[29]; + // const fvec cj15 = fTr.C51 + fTr.C21*J[17] + fTr.C31*J[23] + fTr.C41*J[29]; + const fvec cj25 = fTr.C52 + fTr.C22 * J[17] + fTr.C32 * J[23] + fTr.C42 * J[29]; + const fvec cj35 = fTr.C53 + fTr.C32 * J[17] + fTr.C33 * J[23] + fTr.C43 * J[29]; + const fvec cj45 = fTr.C54 + fTr.C42 * J[17] + fTr.C43 * J[23] + fTr.C44 * J[29]; + const fvec cj55 = fTr.C55 + fTr.C52 * J[17] + fTr.C53 * J[23] + fTr.C54 * J[29]; + + + fTr.C00 = cj00 + cj20 * J[12] + cj30 * J[18] + cj40 * J[24]; + + fTr.C10 = cj10 + cj20 * J[13] + cj30 * J[19] + cj40 * J[25]; + fTr.C11 = cj11 + cj21 * J[13] + cj31 * J[19] + cj41 * J[25]; + + fTr.C20 = cj20 + cj30 * J[20] + cj40 * J[26]; + fTr.C21 = cj21 + cj31 * J[20] + cj41 * J[26]; + fTr.C22 = cj22 + cj32 * J[20] + cj42 * J[26]; + + fTr.C30 = cj30 + cj20 * J[15] + cj40 * J[27]; + fTr.C31 = cj31 + cj21 * J[15] + cj41 * J[27]; + fTr.C32 = cj32 + cj22 * J[15] + cj42 * J[27]; + fTr.C33 = cj33 + cj23 * J[15] + cj43 * J[27]; + + fTr.C40 = cj40; + fTr.C41 = cj41; + fTr.C42 = cj42; + fTr.C43 = cj43; + fTr.C44 = cj44; + + fTr.C50 = cj50 + cj20 * J[17] + cj30 * J[23] + cj40 * J[29]; + fTr.C51 = cj51 + cj21 * J[17] + cj31 * J[23] + cj41 * J[29]; + fTr.C52 = cj52 + cj22 * J[17] + cj32 * J[23] + cj42 * J[29]; + fTr.C53 = cj53 + cj23 * J[17] + cj33 * J[23] + cj43 * J[29]; + fTr.C54 = cj54 + cj24 * J[17] + cj34 * J[23] + cj44 * J[29]; + fTr.C55 = cj55 + cj25 * J[17] + cj35 * J[23] + cj45 * J[29]; } void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w) @@ -801,10 +801,10 @@ void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w) //const fscal RadThick=0.0009f;//0.5/18.76; const fscal logRadThick = log(fPipeRadThick[0]); - fvec tx = ftx; - fvec ty = fty; - fvec txtx = ftx * ftx; - fvec tyty = fty * fty; + fvec tx = fTr.tx; + fvec ty = fTr.ty; + fvec txtx = fTr.tx * fTr.tx; + fvec tyty = fTr.ty * fTr.ty; fvec txtx1 = txtx + ONE; fvec h = txtx + tyty; fvec t = sqrt(txtx1 + tyty); @@ -816,9 +816,9 @@ void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w) //fvec a = ( (ONE+mass2*qp0*qp0t)*RadThick*s0*s0 ); fvec a = ((t + fMass2 * qp0 * qp0t) * fPipeRadThick * s0 * s0); - C22 += w * txtx1 * a; - C32 += w * tx * ty * a; - C33 += w * (ONE + tyty) * a; + fTr.C22 += w * txtx1 * a; + fTr.C32 += w * tx * ty * a; + fTr.C33 += w * (ONE + tyty) * a; } @@ -826,8 +826,8 @@ void L1TrackParFit::L1AddMaterial(const fvec& radThick, fvec qp0, fvec w) { cnst ONE = 1.; - fvec tx = ftx; - fvec ty = fty; + fvec tx = fTr.tx; + fvec ty = fTr.ty; fvec txtx = tx * tx; fvec tyty = ty * ty; fvec txtx1 = txtx + ONE; @@ -842,17 +842,17 @@ void L1TrackParFit::L1AddMaterial(const fvec& radThick, fvec qp0, fvec w) //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 ); fvec a = ((t + fMass2 * qp0 * qp0t) * radThick * s0 * s0); - C22 += w * txtx1 * a; - C32 += w * tx * ty * a; - C33 += w * (ONE + tyty) * a; + fTr.C22 += w * txtx1 * a; + fTr.C32 += w * tx * ty * a; + fTr.C33 += w * (ONE + tyty) * a; } void L1TrackParFit::L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream) { cnst ONE = 1.; - fvec tx = ftx; - fvec ty = fty; + fvec tx = fTr.tx; + fvec ty = fTr.ty; fvec txtx = tx * tx; fvec tyty = ty * ty; fvec txtx1 = txtx + ONE; @@ -872,18 +872,18 @@ void L1TrackParFit::L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thi fvec T23 = (thickness * thickness) / fvec(3.0); fvec T2 = thickness / fvec(2.0); - C00 += w * txtx1 * a * T23; - C10 += w * tx * ty * a * T23; - C20 += w * txtx1 * a * D * T2; - C30 += w * tx * ty * a * D * T2; + fTr.C00 += w * txtx1 * a * T23; + fTr.C10 += w * tx * ty * a * T23; + fTr.C20 += w * txtx1 * a * D * T2; + fTr.C30 += w * tx * ty * a * D * T2; - C11 += w * (ONE + tyty) * a * T23; - C21 += w * tx * ty * a * D * T2; - C31 += w * (ONE + tyty) * a * D * T2; + fTr.C11 += w * (ONE + tyty) * a * T23; + fTr.C21 += w * tx * ty * a * D * T2; + fTr.C31 += w * (ONE + tyty) * a * D * T2; - C22 += w * txtx1 * a; - C32 += w * tx * ty * a; - C33 += w * (ONE + tyty) * a; + fTr.C22 += w * txtx1 * a; + fTr.C32 += w * tx * ty * a; + fTr.C33 += w * (ONE + tyty) * a; } @@ -891,9 +891,9 @@ void L1TrackParFit::L1AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w) { cnst ONE = 1.f; - fvec tx = ftx; - fvec ty = fty; - // fvec time = ft; + fvec tx = fTr.tx; + fvec ty = fTr.ty; + // fvec time = fTr.t; fvec txtx = tx * tx; fvec tyty = ty * ty; fvec txtx1 = txtx + ONE; @@ -908,9 +908,9 @@ void L1TrackParFit::L1AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w) //fvec a = ( (ONE+mass2*qp0*qp0t)*info.RadThick*s0*s0 ); fvec a = ((t + fMass2 * qp0 * qp0t) * info.RadThick * s0 * s0); - C22 += w * txtx1 * a; - C32 += w * tx * ty * a; - C33 += w * (ONE + tyty) * a; + fTr.C22 += w * txtx1 * a; + fTr.C32 += w * tx * ty * a; + fTr.C33 += w * (ONE + tyty) * a; } @@ -923,7 +923,7 @@ void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec d const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2); - fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty); + fvec tr = sqrt(fvec(1.f) + fTr.tx * fTr.tx + fTr.ty * fTr.ty); const fvec dE = bethe * radThick * tr * 2.33f * 9.34961f; @@ -933,13 +933,13 @@ void L1TrackParFit::EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec d corr = iif(ok, corr, fvec::One()); qp0 *= corr; - fqp *= corr; - C40 *= corr; - C41 *= corr; - C42 *= corr; - C43 *= corr; - C44 *= corr * corr; - C54 *= corr; + fTr.qp *= corr; + fTr.C40 *= corr; + fTr.C41 *= corr; + fTr.C42 *= corr; + fTr.C43 *= corr; + fTr.C44 *= corr * corr; + fTr.C54 *= corr; } template<int atomicZ> @@ -958,7 +958,7 @@ void L1TrackParFit::EnergyLossCorrection(float atomicA, float rho, float radLen, const fvec bethe = L1Fit::ApproximateBetheBloch(p2 / fMass2, rho, 0.20, 3.00, i, atomicZ / atomicA); - fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty); + fvec tr = sqrt(fvec(1.f) + fTr.tx * fTr.tx + fTr.ty * fTr.ty); fvec dE = bethe * radThick * tr * radLen * rho; @@ -968,7 +968,7 @@ void L1TrackParFit::EnergyLossCorrection(float atomicA, float rho, float radLen, corr = iif(ok, corr, fvec::One()); qp0 *= corr; - fqp *= corr; + fTr.qp *= corr; fvec P(sqrt(p2)); // GeV @@ -998,12 +998,12 @@ void L1TrackParFit::EnergyLossCorrection(float atomicA, float rho, float radLen, fvec P2 = P * P; fvec SDEDX = (E2 * DEDX2) / (P2 * P2 * P2); - // T.C40 *= corr; - // T.C41 *= corr; - // T.C42 *= corr; - // T.C43 *= corr; - // T.C44 *= corr*corr; - C44 += abs(SDEDX); + // T.fTr.C40 *= corr; + // T.fTr.C41 *= corr; + // T.fTr.C42 *= corr; + // T.fTr.C43 *= corr; + // T.fTr.C44 *= corr*corr; + fTr.C44 += abs(SDEDX); } diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h index 904bf44ad3a40bb9097a253e8cc89602ded86be2..e66c7eb0b060c9004b9b12543e24777e2a582b07 100644 --- a/reco/L1/L1Algo/L1TrackParFit.h +++ b/reco/L1/L1Algo/L1TrackParFit.h @@ -15,10 +15,8 @@ class L1TrackParFit { public: - fvec fx {0.}, fy {0.}, ftx {0.}, fty {0.}, fqp {0.}, fz {0.}, ft {0.}, C00 {0.}, C10 {0.}, C11 {0.}, C20 {0.}, - C21 {0.}, C22 {0.}, C30 {0.}, C31 {0.}, C32 {0.}, C33 {0.}, C40 {0.}, C41 {0.}, C42 {0.}, C43 {0.}, C44 {0.}, - C50 {0.}, C51 {0.}, C52 {0.}, C53 {0.}, C54 {0.}, C55 {0.}, chi2 {0.}, NDF {0.}; - // fvec n; + L1TrackPar fTr {}; + fvec fQp0; fvec fMass = 0.10565800; // muon mass fvec fMass2 = fMass * fMass; // mass squared @@ -28,37 +26,7 @@ public: L1TrackParFit() = default; - L1TrackParFit(double* T, double* C) - : fx(T[0]) - , fy(T[1]) - , ftx(T[3]) - , fty(T[4]) - , fqp(T[5]) - , fz(T[6]) - , ft(T[7]) - , C00(C[0]) - , C10(C[1]) - , C11(C[2]) - , C20(C[3]) - , C21(C[4]) - , C22(C[5]) - , C30(C[6]) - , C31(C[7]) - , C32(C[8]) - , C33(C[9]) - , C40(C[10]) - , C41(C[11]) - , C42(C[12]) - , C43(C[13]) - , C44(C[14]) - , C50(C[15]) - , C51(C[16]) - , C52(C[17]) - , C53(C[18]) - , C54(C[19]) - , C55(C[20]) - , chi2(0.) - , NDF(0.) {}; + L1TrackParFit(double* T, double* C) : fTr(T, C) {} void SetOneEntry(const int i0, const L1TrackParFit& T1, const int i1); @@ -121,144 +89,13 @@ public: // ============================================================================================= -inline void L1TrackParFit::Print(int i) -{ - std::ios_base::fmtflags coutFlags(std::cout.flags()); - std::cout.setf(std::ios::scientific, std::ios::floatfield); - - if (i == -1) { - std::cout << "T = " << std::endl; - std::cout << " x " << fx << std::endl; - std::cout << " y " << fy << std::endl; - std::cout << " tx " << ftx << std::endl; - std::cout << " ty " << fty << std::endl; - std::cout << " qp " << fqp << std::endl; - std::cout << " t " << ft << std::endl; - std::cout << " z " << fz << std::endl; - std::cout << "C = " << std::endl; - std::cout << " c00 " << C00 << std::endl; - std::cout << " c11 " << C11 << std::endl; - std::cout << " c22 " << C22 << std::endl; - std::cout << " c33 " << C33 << std::endl; - std::cout << " c44 " << C44 << std::endl; - std::cout << " c55 " << C55 << std::endl; - } - else { - std::cout << "T = "; - std::cout << fx[i] << " "; - std::cout << fy[i] << " "; - std::cout << ftx[i] << " "; - std::cout << fty[i] << " "; - std::cout << fqp[i] << " "; - std::cout << fz[i] << std::endl; - std::cout << "C = "; - std::cout << C00[i] << " "; - std::cout << C11[i] << " "; - std::cout << C22[i] << " "; - std::cout << C33[i] << " "; - std::cout << C44[i] << " "; - std::cout << C55[i] << std::endl; - } - std::cout.flags(coutFlags); -} - - -inline void L1TrackParFit::Compare(L1TrackPar& T) -{ - std::cout.precision(8); - - // if (fabs(T.x[0]-fx[0])/T.x[0] > 1.e-7) std::cout << fx <<" x "<< T.x << std::endl; - // - // if (fabs(T.y[0]-fy[0])/T.y[0] > 1.e-7) std::cout << fy <<" y "<< T.y << std::endl; - // - // if (fabs(T.tx[0]-ftx[0])/T.tx[0] > 1.e-7) std::cout << ftx <<" tx "<< T.tx << std::endl; - // - // if (fabs(T.ty[0]-fty[0])/T.ty[0] > 1.e-7) std::cout << fty <<" ty "<< T.ty << std::endl; - // - // if (fabs(T.qp[0]-fqp[0])/T.qp[0] > 1.e-7) std::cout << fqp <<" qp "<< T.qp << std::endl; - // - // if (fabs(T.z[0]-fz[0])/T.z[0] > 1.e-7) std::cout << fz <<" z "<< T.z << std::endl; - // - // if (fabs(T.C00[0]-C00[0])/T.C00[0] > 1.e-7) std::cout << C00 <<" C00 "<< T.C00 << std::endl; - // if (fabs(T.C11[0]-C11[0])/T.C11[0] > 1.e-7) std::cout << C11 <<" C11 "<< T.C11 << std::endl; - // if (fabs(T.C22[0]-C22[0])/T.C22[0] > 1.e-7) std::cout << C22 <<" C22 "<< T.C22 << std::endl; - // if (fabs(T.C33[0]-C33[0])/T.C33[0] > 1.e-7) std::cout << C33 <<" C33 "<< T.C33 << std::endl; - // if (fabs(T.C44[0]-C44[0])/T.C44[0] > 1.e-7) std::cout << C44 <<" C44 "<< T.C44 << std::endl; - +inline void L1TrackParFit::Print(int i) { fTr.Print(i); } - // if (!(T.x==fx)[0]) std::cout << fx <<" x "<< T.x << std::endl; - // if (!(T.y==fy)[0]) std::cout << fy << " y "<< T.y << std::endl; - // if (!(T.tx==ftx)[0]) std::cout << ftx << " ty "<< T.tx << std::endl; - // if (!(T.ty==fty)[0]) std::cout << fty << " ty "<< T.ty << std::endl; - // if (!(T.qp==fqp)[0]) std::cout << fqp << " qp "<< T.qp << std::endl; - // if (!(T.z==fz)[0]) std::cout << fz << " z "<< T.z << std::endl; - // // if (T.t!=ft) std::cout << ft << " fx "<< T.x << std::endl; - - - // if (!(T.C00==C00)[0]) std::cout << C00 << " C00 "<< T.C00 << std::endl; - // if (!(T.C11==C11)[0]) std::cout << C11 << " C11 "<< T.C11 << std::endl; - // if (!(T.C22==C22)[0]) std::cout << C22 << " C22 "<< T.C22 << std::endl; - // if (!(T.C33==C33)[0]) std::cout << C33 << " C33 "<< T.C33 << std::endl; - // if (!(T.C44==C44)[0]) std::cout << C44 << " C44 "<< T.C44 << std::endl; - // if (T.C55!=C55) std::cout << C55 << " fx "<< T.C55 << std::endl; - - - std::cout << "parameters:" << std::endl; - std::cout << T.x[0] << " " << T.y[0] << " " << T.z[0] << " " << T.tx[0] << " " << T.ty[0] << " " << T.qp[0] - << std::endl; - std::cout << fx[0] << " " << fy[0] << " " << fz[0] << " " << ftx[0] << " " << fty[0] << " " << fqp[0] << " " << ft[0] - << std::endl; - - std::cout << "Covariance matrix:" << std::endl; - std::cout << T.C00[0] << " " << T.C10[0] << " " << T.C11[0] << " " << T.C20[0] << " " << T.C21[0] << " " << T.C22[0] - << " " << T.C30[0] << " " << T.C31[0] << " " << T.C32[0] << " " << T.C33[0] << " " << T.C40[0] << " " - << T.C41[0] << " " << T.C42[0] << " " << T.C43[0] << " " << T.C44[0] << std::endl; - - std::cout << C00[0] << " " << C10[0] << " " << C11[0] << " " << C20[0] << " " << C21[0] << " " << C22[0] << " " - << C30[0] << " " << C31[0] << " " << C32[0] << " " << C33[0] << " " << C40[0] << " " << C41[0] << " " - << C42[0] << " " << C43[0] << " " << C44[0] << std::endl; - std::cout << " Time covariance:" << std::endl; - std::cout << " " << C50[0] << " " << C51[0] << " " << C52[0] << " " << C53[0] << " " << C54[0] << " " << C55[0] - << " " << std::endl; - - std::cout << std::endl; -} inline void L1TrackParFit::SetOneEntry(const int i0, const L1TrackParFit& T1, const int i1) { - fx[i0] = T1.fx[i1]; - fy[i0] = T1.fy[i1]; - ftx[i0] = T1.ftx[i1]; - fty[i0] = T1.fty[i1]; - fqp[i0] = T1.fqp[i1]; - fz[i0] = T1.fz[i1]; - ft[i0] = T1.ft[i1]; - C00[i0] = T1.C00[i1]; - C10[i0] = T1.C10[i1]; - C11[i0] = T1.C11[i1]; - C20[i0] = T1.C20[i1]; - C21[i0] = T1.C21[i1]; - C22[i0] = T1.C22[i1]; - C30[i0] = T1.C30[i1]; - C31[i0] = T1.C31[i1]; - C32[i0] = T1.C32[i1]; - C33[i0] = T1.C33[i1]; - C40[i0] = T1.C40[i1]; - C41[i0] = T1.C41[i1]; - C42[i0] = T1.C42[i1]; - C43[i0] = T1.C43[i1]; - C44[i0] = T1.C44[i1]; - C50[i0] = T1.C50[i1]; - C51[i0] = T1.C51[i1]; - C52[i0] = T1.C52[i1]; - C53[i0] = T1.C53[i1]; - C54[i0] = T1.C54[i1]; - C55[i0] = T1.C55[i1]; - - chi2[i0] = T1.chi2[i1]; - NDF[i0] = T1.NDF[i1]; - // time[i0] = T1.time[i1]; - // n[i0] = T1.n[i1]; -} // SetOneEntry + fTr.SetOneEntry(i0, T1.fTr, i1); + fQp0[i0] = T1.fQp0[i1]; +} #endif