diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx index 8c6e07cc59d037332b08bf6668fcf43ce49f7a04..c4695c66dfbf3f8991c9cf29e427d0a62d5131a1 100644 --- a/reco/L1/CbmL1ReadEvent.cxx +++ b/reco/L1/CbmL1ReadEvent.cxx @@ -1261,17 +1261,18 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int assert(th.iStripB >= 0 || th.iStripB < NStrips); L1Hit h; + h.iSt = th.iStation; h.f = th.iStripF; h.b = th.iStripB; h.ID = th.id; - h.t = th.time; - h.dt = th.dt; - h.du = th.du; - h.dv = th.dv; + h.z = th.z; h.u = th.u; h.v = th.v; - h.z = th.z; - h.iSt = th.iStation; + h.t = th.time; + h.dt2 = th.dt * th.dt; + h.du2 = th.du * th.du; + h.dv2 = th.dv * th.dv; + // save hit fvExternalHits.push_back(CbmL1Hit(iHit, th.ExtIndex, th.Det)); diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx index 71fc7a1a9b7fe5a02795c0d32ba0d17449b95552..4699109dc3a9d6f96b4f7915fb1a42135f693c74 100644 --- a/reco/L1/L1Algo/L1Algo.cxx +++ b/reco/L1/L1Algo/L1Algo.cxx @@ -180,12 +180,12 @@ L1HitPoint L1Algo::CreateHitPoint(const L1Hit& hit) { /// full the hit point by hit information: takes hit as input (2 strips) /// and creates hit_point with all coordinates (x,y,z,u,v,t); - return L1HitPoint(hit.z, hit.u, hit.v, hit.du, hit.dv, hit.t, hit.dt); + return L1HitPoint(hit.z, hit.u, hit.v, hit.t, hit.du2, hit.dv2, hit.dt2); } void L1Algo::CreateHitPoint(const L1Hit& hit, L1HitPoint& point) { - point.Set(hit.z, hit.u, hit.v, hit.du, hit.dv, hit.t, hit.dt); + point.Set(hit.z, hit.u, hit.v, hit.t, hit.du2, hit.dv2, hit.dt2); } int L1Algo::GetMcTrackIdForHit(int iHit) const diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h index f30ce0aaedc63633daab5154847a4b83fd433f0b..f00d23276ebde9e27c52efe530af406563c90af6 100644 --- a/reco/L1/L1Algo/L1Algo.h +++ b/reco/L1/L1Algo/L1Algo.h @@ -264,16 +264,16 @@ public: void findSingletsStep0( // input Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l, // output - fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* HitTime_l, fvec* HitTimeEr, fvec* d_u, - fvec* d_v); + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* t_l, fvec* dt2_l, fvec* du2_l, + fvec* dv2_l); /// Get the field approximation. Add the target to parameters estimation. Propagate to middle station. void findSingletsStep1( // input int istal, int istam, Tindex n1_V, - fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr, + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* t_l, fvec* dt2_l, // output - L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_u, fvec* d_v); + L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* du2_l, fvec* dv2_l); /// Find the doublets. Reformat data in the portion of doublets. void findDoubletsStep0( // input @@ -298,12 +298,12 @@ public: // output Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3, - L1Vector<fvec>& du_, L1Vector<fvec>& dv_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER); + L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3); /// Add the right hits to parameters estimation. void findTripletsStep1( // input Tindex n3_V, const L1Station& star, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3, - L1Vector<fvec>& du_, L1Vector<fvec>& dv_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER, + L1Vector<fvec>& du2_3, L1Vector<fvec>& dv2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3, // output L1Vector<L1TrackPar>& T_3); @@ -346,11 +346,11 @@ public: void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st); void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st); - void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st); + void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st); - void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st, fvec& dx, fvec& dy, + void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& dx2, fvec& dy2, fvec& dxy); - void FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& t_er, L1Station& st, fvec& dx, fvec& dy, + void FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& dx2, fvec& dy2, fvec& dxy); diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx index 73bed78e94f0735bb08478bec8f519cc0d4864cd..3ed23c719478b279dd2e34554a89ecb381ef3a69 100644 --- a/reco/L1/L1Algo/L1BranchExtender.cxx +++ b/reco/L1/L1Algo/L1BranchExtender.cxx @@ -88,7 +88,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, T.C10 = sta0.XYInfo.C10; T.C11 = sta0.XYInfo.C11; - if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du, hit0.dv); } + if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2); } T.C20 = T.C21 = 0; T.C30 = T.C31 = T.C32 = 0; @@ -96,7 +96,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = 0; T.C22 = T.C33 = vINF; T.C44 = 1.; - T.C55 = hit0.dt * hit0.dt; + T.C55 = hit0.dt2; L1FieldValue fldB0, fldB1, fldB2 _fvecalignment; L1FieldRegion fld _fvecalignment; @@ -141,15 +141,15 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir, L1UMeasurementInfo info = sta.frontInfo; - if (fUseHitErrors) { info.sigma2 = hit.du * hit.du; } + if (fUseHitErrors) { info.sigma2 = hit.du2; } L1Filter(T, info, u); info = sta.backInfo; - if (fUseHitErrors) { info.sigma2 = hit.dv * hit.dv; } + if (fUseHitErrors) { info.sigma2 = hit.dv2; } L1Filter(T, info, v); - FilterTime(T, hit.t, hit.dt); + FilterTime(T, hit.t, hit.dt2); fldB0 = fldB1; fldB1 = fldB2; @@ -269,7 +269,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir, ih += HitsUnusedStartIndex[ista]; const L1Hit& hit = (*vHitsUnused)[ih]; //TODO: bug, it should be hit.dt*hit.dt - if (fabs(hit.t - T.t[0]) > sqrt(T.C55[0] + hit.dt * hit.dt) * 5) continue; + if (fabs(hit.t - T.t[0]) > sqrt(T.C55[0] + hit.dt2) * 5) continue; //if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used //L1_SHOW(fvHitKeyFlags.size()); @@ -320,15 +320,15 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir, L1UMeasurementInfo info = sta.frontInfo; - if (fUseHitErrors) { info.sigma2 = hit.du * hit.du; } + if (fUseHitErrors) { info.sigma2 = hit.du2; } L1Filter(T, info, u); info = sta.backInfo; - if (fUseHitErrors) { info.sigma2 = hit.dv * hit.dv; } + if (fUseHitErrors) { info.sigma2 = hit.dv2; } L1Filter(T, info, v); - FilterTime(T, hit.t, hit.dt); + FilterTime(T, hit.t, hit.dt2); fldB0 = fldB1; fldB1 = fldB2; diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx index 5facd7a142a6a8a5ffe82d4655d3a729955f4e6c..3ac5105e2ff00862e705ecc0a590c6603c5c66dc 100644 --- a/reco/L1/L1Algo/L1CATrackFinder.cxx +++ b/reco/L1/L1Algo/L1CATrackFinder.cxx @@ -128,8 +128,7 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc inline void L1Algo::findSingletsStep0( // input Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l, // output - fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* HitTime_l, fvec* HitTimeEr, fvec* d_u, - fvec* d_v) + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* t_l, fvec* dt2_l, fvec* du2_l, fvec* dv2_l) { /// Prepare the portion of data of left hits of a triplet: @@ -141,12 +140,12 @@ inline void L1Algo::findSingletsStep0( // input if (lastV >= 0) { // set some positive errors to unfilled part of vectors in order to avoid nans L1HitPoint& hitl = Hits_l[0]; - d_u[lastV] = hitl.dU(); - d_v[lastV] = hitl.dV(); - HitTimeEr[lastV] = hitl.timeEr; + du2_l[lastV] = hitl.dU2(); + dv2_l[lastV] = hitl.dV2(); + dt2_l[lastV] = hitl.dT2(); u_front_l[lastV] = hitl.U(); u_back_l[lastV] = hitl.V(); - HitTime_l[lastV] = hitl.time; + t_l[lastV] = hitl.T(); zPos_l[lastV] = hitl.Z(); } @@ -156,16 +155,16 @@ inline void L1Algo::findSingletsStep0( // input L1HitPoint& hitl = Hits_l[ilh]; - HitTime_l[i1_V][i1_4] = hitl.time; - HitTimeEr[i1_V][i1_4] = hitl.timeEr; + t_l[i1_V][i1_4] = hitl.T(); + dt2_l[i1_V][i1_4] = hitl.dT2(); hitsl[i1] = ilh; u_front_l[i1_V][i1_4] = hitl.U(); u_back_l[i1_V][i1_4] = hitl.V(); if (fUseHitErrors) { - d_u[i1_V][i1_4] = hitl.dU(); - d_v[i1_V][i1_4] = hitl.dV(); + du2_l[i1_V][i1_4] = hitl.dU2(); + dv2_l[i1_V][i1_4] = hitl.dV2(); } zPos_l[i1_V][i1_4] = hitl.Z(); @@ -177,10 +176,10 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search int istal, int istam, /// indexes of left and middle stations of a triplet Tindex n1_V, /// - fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr, + fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* t_l, fvec* dt2_l, // output // L1TrackPar *T_1, - L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_u, fvec* d_v) + L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* du2_l, fvec* dv2_l) { /// Get the field approximation. Add the target to parameters estimation. @@ -244,8 +243,8 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search fvec& v = u_back_l[i1_V]; auto [xl, yl] = stal.ConvUVtoXY<fvec>(u, v); fvec zl = zPos_l[i1_V]; - fvec& time = HitTime_l[i1_V]; - fvec& timeEr = HitTimeEr[i1_V]; + fvec& time = t_l[i1_V]; + fvec& timeEr2 = dt2_l[i1_V]; const fvec dzli = 1.f / (zl - fTargZ); const fvec tx = (xl - fTargX) * dzli; @@ -283,7 +282,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search T.C44 = fMaxInvMom / fvec(3.) * fMaxInvMom / fvec(3.); - T.C55 = timeEr * timeEr; + T.C55 = timeEr2; { // add the target constraint T.x = xl; @@ -293,7 +292,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search T.C10 = stal.XYInfo.C10; T.C11 = stal.XYInfo.C11; - if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(d_u[i1_V], d_v[i1_V]); } + if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]); } //assert(T.IsConsistent(true, -1)); @@ -406,14 +405,14 @@ inline void L1Algo::findDoubletsStep0( // Pick_m22 is not used, search for mean squared, 2nd version // -- collect possible doublets -- - const fscal iz = 1.f / (T1.z[i1_4] - fParameters.GetTargetPositionZ()[0]); - const fscal timeError = T1.C55[i1_4]; - const fscal time = T1.t[i1_4]; + const fscal iz = 1.f / (T1.z[i1_4] - fParameters.GetTargetPositionZ()[0]); + const fscal timeError2 = T1.C55[i1_4]; + const fscal time = T1.t[i1_4]; L1HitAreaTime areaTime(vGridTime[iStaM], T1.x[i1_4] * iz, T1.y[i1_4] * iz, (sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + fMaxDZ * abs(T1.tx))[i1_4] * iz, (sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + fMaxDZ * abs(T1.ty))[i1_4] * iz, time, - sqrt(timeError) * 5); + sqrt(timeError2) * 5); for (L1HitIndex_t imh = -1; true;) { // loop over the hits in the area if (fParameters.DevIsIgnoreHitSearchAreas()) { @@ -435,8 +434,8 @@ inline void L1Algo::findDoubletsStep0( // check y-boundaries //TODO: move hardcoded cuts to parameters if ((stam.timeInfo) && (stal.timeInfo)) { - if (fabs(time - hitm.time) > sqrt(timeError + hitm.timeEr * hitm.timeEr) * 5) continue; - if (fabs(time - hitm.time) > 30) continue; + if (fabs(time - hitm.T()) > sqrt(timeError2 + hitm.dT2()) * 5) continue; + if (fabs(time - hitm.T()) > 30) continue; } // - check whether hit belong to the window ( track position +\- errors ) - @@ -458,7 +457,7 @@ inline void L1Algo::findDoubletsStep0( fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + stam.XYInfo.C11[i1_4]); /// Covariation matrix of the hit - auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU(), hitm.dV()); + auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU2(), hitm.dV2()); if (fUseHitErrors) { dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + dyyScalMhit); } @@ -487,7 +486,7 @@ inline void L1Algo::findDoubletsStep0( L1UMeasurementInfo info = stam.frontInfo; - if (fUseHitErrors) info.sigma2 = hitm.dU() * hitm.dU(); + if (fUseHitErrors) info.sigma2 = hitm.dU2(); L1FilterChi2XYC00C10C11(info, x, y, C00, C10, C11, chi2, hitm.U()); @@ -498,11 +497,11 @@ inline void L1Algo::findDoubletsStep0( info = stam.backInfo; - if (fUseHitErrors) info.sigma2 = hitm.dV() * hitm.dV(); + if (fUseHitErrors) info.sigma2 = hitm.dV2(); L1FilterChi2(info, x, y, C00, C10, C11, chi2, hitm.V()); - // FilterTime(T1, hitm.time, hitm.timeEr); + // FilterTime(T1, hitm.T(), hitm.dT2()); if (!fpCurrentIteration->GetTrackFromTripletsFlag()) { if (chi2[i1_4] > fDoubletChi2Cut) continue; @@ -543,7 +542,7 @@ inline void L1Algo::findTripletsStep0( // input // output Tindex& n3, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3, L1Vector<fvec>& u_front_3, L1Vector<fvec>& u_back_3, L1Vector<fvec>& z_Pos_3, - L1Vector<fvec>& dv_, L1Vector<fvec>& du_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER) + L1Vector<fvec>& dv2_3, L1Vector<fvec>& du2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3) { int iStaM = &stam - fParameters.GetStations().begin(); int iStaR = &star - fParameters.GetStations().begin(); @@ -572,10 +571,10 @@ inline void L1Algo::findTripletsStep0( // input u_front_3.reset(1, fvec::Zero()); u_back_3.reset(1, fvec::Zero()); z_Pos_3.reset(1, fvec::Zero()); - du_.reset(1, fvec::One()); - dv_.reset(1, fvec::One()); - timeR.reset(1, fvec::Zero()); - timeER.reset(1, fvec::One()); + du2_3.reset(1, fvec::One()); + dv2_3.reset(1, fvec::One()); + t_3.reset(1, fvec::Zero()); + dt2_3.reset(1, fvec::One()); assert(istar < fParameters.GetNstationsActive()); //TODO SG!!! check if it is needed @@ -589,11 +588,11 @@ inline void L1Algo::findTripletsStep0( // input // pack the data fvec u_front_2 = 0.f; fvec u_back_2 = 0.f; - fvec du2 = 1.f; - fvec dv2 = 1.f; + fvec du2_2 = 1.f; + fvec dv2_2 = 1.f; fvec zPos_2 = 0.f; - fvec timeM = 0.f; - fvec timeMEr = 1.f; + fvec t_2 = 0.f; + fvec dt2_2 = 1.f; size_t n2_4 = 0; for (; n2_4 < fvec::size() && i2 < n2; i2++, n2_4++) { @@ -616,11 +615,11 @@ inline void L1Algo::findTripletsStep0( // input u_front_2[n2_4] = hitm.U(); u_back_2[n2_4] = hitm.V(); zPos_2[n2_4] = hitm.Z(); - timeM[n2_4] = hitm.time; - timeMEr[n2_4] = hitm.timeEr; + t_2[n2_4] = hitm.T(); + dt2_2[n2_4] = hitm.dT2(); // num[n2_4] = hitm.track; - du2[n2_4] = hitm.dU(); - dv2[n2_4] = hitm.dV(); + du2_2[n2_4] = hitm.dU2(); + dv2_2[n2_4] = hitm.dV2(); hitsl_2[n2_4] = hitsl_1[i1]; hitsm_2_tmp[n2_4] = hitsm_2[i2]; @@ -641,7 +640,7 @@ inline void L1Algo::findTripletsStep0( // input L1UMeasurementInfo info = stam.frontInfo; - if (fUseHitErrors) info.sigma2 = du2 * du2; + if (fUseHitErrors) info.sigma2 = du2_2; // TODO: SG: L1FilterNoField is wrong. // TODO: If the field was present before, @@ -677,7 +676,7 @@ inline void L1Algo::findTripletsStep0( // input */ info = stam.backInfo; - if (fUseHitErrors) info.sigma2 = dv2 * dv2; + if (fUseHitErrors) info.sigma2 = dv2_2; if (istam < fNfieldStations) { L1Filter(T2, info, u_back_2); } else { @@ -686,7 +685,7 @@ inline void L1Algo::findTripletsStep0( // input // assert(T2.IsConsistent(true, n2_4)); - FilterTime(T2, timeM, timeMEr, stam.timeInfo); + FilterTime(T2, t_2, dt2_2, stam.timeInfo); // assert(T2.IsConsistent(true, n2_4)); @@ -735,9 +734,9 @@ inline void L1Algo::findTripletsStep0( // input if (fabs(T2.tx[i2_4]) > fMaxSlope) continue; if (fabs(T2.ty[i2_4]) > fMaxSlope) continue; - const fvec Pick_r22 = fTripletChi2Cut - T2.chi2 + (!fpCurrentIteration->GetTrackFromTripletsFlag() ? 0 : 1); - const fscal timeError = T2.C55[i2_4]; - const fscal time = T2.t[i2_4]; + const fvec Pick_r22 = fTripletChi2Cut - T2.chi2 + (!fpCurrentIteration->GetTrackFromTripletsFlag() ? 0 : 1); + const fscal timeError2 = T2.C55[i2_4]; + const fscal time = T2.t[i2_4]; // find first possible hit @@ -745,7 +744,7 @@ inline void L1Algo::findTripletsStep0( // input L1HitAreaTime area(vGridTime[&star - fParameters.GetStations().begin()], T2.x[i2_4] * iz, T2.y[i2_4] * iz, (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + fMaxDZ * abs(T2.tx))[i2_4] * iz, (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + fMaxDZ * abs(T2.ty))[i2_4] * iz, time, - sqrt(timeError) * 5); + sqrt(timeError2) * 5); L1HitIndex_t irh = 0; L1HitIndex_t doubletNtriplets = 0; @@ -783,11 +782,11 @@ inline void L1Algo::findTripletsStep0( // input L1ExtrapolateLine(T_cur, zr); if ((star.timeInfo) && (stam.timeInfo)) - if (fabs(T_cur.t[i2_4] - hitr.time) > sqrt(T_cur.C55[i2_4] + hitr.timeEr) * 5) continue; + if (fabs(T_cur.t[i2_4] - hitr.T()) > sqrt(T_cur.C55[i2_4] + sqrt(hitr.dT2())) * 5) continue; // TODO: SG: hardcoded cut of 30 ns if ((star.timeInfo) && (stam.timeInfo)) - if (fabs(T_cur.t[i2_4] - hitr.time) > 30) continue; + if (fabs(T_cur.t[i2_4] - hitr.T()) > 30) continue; // - check whether hit belong to the window ( track position +\- errors ) - // check lower boundary @@ -801,7 +800,7 @@ inline void L1Algo::findTripletsStep0( // input + star.XYInfo.C11[i2_4]))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets /// Covariation matrix of the hit - auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU(), hitr.dV()); + auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU2(), hitr.dV2()); if (fUseHitErrors) { dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + dyyScalRhit))); } @@ -829,16 +828,16 @@ inline void L1Algo::findTripletsStep0( // input info = star.frontInfo; - if (fUseHitErrors) info.sigma2 = hitr.dU() * hitr.dU(); + if (fUseHitErrors) info.sigma2 = hitr.dU2(); L1FilterChi2XYC00C10C11(info, x, y, C00, C10, C11, chi2, hitr.U()); info = star.backInfo; - if (fUseHitErrors) info.sigma2 = hitr.dV() * hitr.dV(); + if (fUseHitErrors) info.sigma2 = hitr.dV2(); L1FilterChi2(info, x, y, C00, C10, C11, chi2, hitr.V()); - FilterTime(T_cur, hitr.time, hitr.timeEr, star.timeInfo); + FilterTime(T_cur, hitr.T(), hitr.dT2(), star.timeInfo); if (!fpCurrentIteration->GetTrackFromTripletsFlag()) { @@ -877,11 +876,11 @@ inline void L1Algo::findTripletsStep0( // input T3.SetOneEntry(n3_4, T2, i2_4); u_front_3[n3_V][n3_4] = hitr.U(); u_back_3[n3_V][n3_4] = hitr.V(); - du_[n3_V][n3_4] = hitr.dU(); - dv_[n3_V][n3_4] = hitr.dV(); + du2_3[n3_V][n3_4] = hitr.dU2(); + dv2_3[n3_V][n3_4] = hitr.dV2(); z_Pos_3[n3_V][n3_4] = zr; - timeR[n3_V][n3_4] = hitr.time; - timeER[n3_V][n3_4] = hitr.timeEr; + t_3[n3_V][n3_4] = hitr.T(); + dt2_3[n3_V][n3_4] = hitr.dT2(); n3++; n3_V = n3 / fvec::size(); @@ -892,10 +891,10 @@ inline void L1Algo::findTripletsStep0( // input u_front_3.push_back(fvec::Zero()); u_back_3.push_back(fvec::Zero()); z_Pos_3.push_back(fvec::Zero()); - du_.push_back(fvec::Zero()); - dv_.push_back(fvec::Zero()); - timeR.push_back(fvec::Zero()); - timeER.push_back(fvec::Zero()); + du2_3.push_back(fvec::One()); + dv2_3.push_back(fvec::One()); + t_3.push_back(fvec::Zero()); + dt2_3.push_back(fvec::One()); } } // search area @@ -906,7 +905,7 @@ inline void L1Algo::findTripletsStep0( // input /// Add the right hits to parameters estimation. inline void L1Algo::findTripletsStep1( // input Tindex n3_V, const L1Station& star, L1Vector<fvec>& u_front_, L1Vector<fvec>& u_back_, L1Vector<fvec>& z_Pos, - L1Vector<fvec>& dv_, L1Vector<fvec>& du_, L1Vector<fvec>& timeR, L1Vector<fvec>& timeER, + L1Vector<fvec>& dv2_3, L1Vector<fvec>& du2_3, L1Vector<fvec>& t_3, L1Vector<fvec>& dt2_3, // output // L1TrackPar *T_3 L1Vector<L1TrackPar>& T_3) @@ -928,7 +927,7 @@ inline void L1Algo::findTripletsStep1( // input L1UMeasurementInfo info = star.frontInfo; - if (fUseHitErrors) info.sigma2 = du_[i3_V] * du_[i3_V]; + if (fUseHitErrors) info.sigma2 = du2_3[i3_V]; bool noField = (&star - fParameters.GetStations().begin() >= fNfieldStations); @@ -941,7 +940,7 @@ inline void L1Algo::findTripletsStep1( // input info = star.backInfo; - if (fUseHitErrors) info.sigma2 = dv_[i3_V] * dv_[i3_V]; + if (fUseHitErrors) info.sigma2 = dv2_3[i3_V]; if (noField) { L1FilterNoField(T3, info, u_back_[i3_V]); } else { @@ -950,9 +949,7 @@ inline void L1Algo::findTripletsStep1( // input // assert(T3.IsConsistent(true, -1)); - if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) { - FilterTime(T3, timeR[i3_V], timeER[i3_V], star.timeInfo); - } + if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) { FilterTime(T3, t_3[i3_V], dt2_3[i3_V], star.timeInfo); } // assert(T3.IsConsistent(true, -1)); // FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V]); @@ -1018,7 +1015,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar } fscal u[NHits], v[NHits], t[NHits], x[NHits], y[NHits], z[NHits]; - fscal du[NHits], dv[NHits], dt[NHits], du2[NHits], dv2[NHits], dt2[NHits]; + fscal du2[NHits], dv2[NHits], dt2[NHits]; for (int ih = 0; ih < NHits; ++ih) { const L1Hit& hit = fInputData.GetHit(ihit[ih]); @@ -1027,13 +1024,10 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar t[ih] = hit.t; std::tie(x[ih], y[ih]) = sta[ih].ConvUVtoXY<fscal>(u[ih], v[ih]); z[ih] = hit.z; - du[ih] = hit.du; - dv[ih] = hit.dv; - dt[ih] = hit.dt; - du2[ih] = hit.du * hit.du; - dv2[ih] = hit.dv * hit.dv; - dt2[ih] = hit.dt * hit.dt; + du2[ih] = hit.du2; + dv2[ih] = hit.dv2; + dt2[ih] = hit.dt2; }; // find the field along the track @@ -1080,11 +1074,11 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar T.z = z[ih0]; T.t = t[ih0]; - //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du[ih0], dv[ih0]); + //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]); fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One()); fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One()); - fit.FilterTime(t[ih0], dt[ih0], fvec::One(), sta[ih0].timeInfo); + fit.FilterTime(t[ih0], dt2[ih0], fvec::One(), sta[ih0].timeInfo); { // add the target constraint fvec eX, eY, J04, J14; @@ -1112,7 +1106,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar //if (ista[ih] == fNstationsBeforePipe) { fit.AddPipeMaterial(qp0, fvec::One()); } fit.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One()); fit.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One()); - fit.FilterTime(t[ih], dt[ih], fvec::One(), sta[ih].timeInfo); + fit.FilterTime(t[ih], dt2[ih], fvec::One(), sta[ih].timeInfo); } } @@ -1134,7 +1128,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar T.t = t[ih0]; T.C55 = dt2[ih0]; - //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du[ih0], dv[ih0]); + //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]); fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One()); fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One()); //fit.FilterTime(t[ih0], dt[ih0], fvec::One(), sta[ih0].timeInfo); @@ -1422,18 +1416,18 @@ inline void L1Algo::CreatePortionOfDoublets( fvec u_front[L1Constants::size::kSingletPortionSizeVec]; fvec u_back[L1Constants::size::kSingletPortionSizeVec]; - fvec dv0[L1Constants::size::kSingletPortionSizeVec]; - fvec du0[L1Constants::size::kSingletPortionSizeVec]; + fvec dv2[L1Constants::size::kSingletPortionSizeVec]; + fvec du2[L1Constants::size::kSingletPortionSizeVec]; fvec zPos[L1Constants::size::kSingletPortionSizeVec]; - fvec HitTime[L1Constants::size::kSingletPortionSizeVec]; - fvec HitTimeEr[L1Constants::size::kSingletPortionSizeVec]; + fvec t[L1Constants::size::kSingletPortionSizeVec]; + fvec dt2[L1Constants::size::kSingletPortionSizeVec]; /// prepare the portion of left hits data findSingletsStep0( // input iSingletPortion * L1Constants::size::kSingletPortionSize, singletPortionSize, vHits_l, // output - u_front, u_back, zPos, hitsl_1, HitTime, HitTimeEr, du0, dv0); + u_front, u_back, zPos, hitsl_1, t, dt2, du2, dv2); for (Tindex i = 0; i < singletPortionSize; ++i) L1_ASSERT(hitsl_1[i] < HitsUnusedStopIndex[istal] - HitsUnusedStartIndex[istal], @@ -1443,9 +1437,9 @@ inline void L1Algo::CreatePortionOfDoublets( /// Get the field approximation. Add the target to parameters estimation. Propagaete to middle station. - findSingletsStep1(istal, istam, portionSize_V, u_front, u_back, zPos, HitTime, HitTimeEr, + findSingletsStep1(istal, istam, portionSize_V, u_front, u_back, zPos, t, dt2, // output - T_1, fld_1, du0, dv0); + T_1, fld_1, du2, dv2); /// Find the doublets. Reformat data in the portion of doublets. @@ -1540,10 +1534,10 @@ inline void L1Algo::CreatePortionOfTriplets( L1Vector<fvec>& u_front3 = fTripletHitR_Ufront[Thread]; L1Vector<fvec>& u_back3 = fTripletHitR_Uback[Thread]; L1Vector<fvec>& z_pos3 = fTripletHitR_Z[Thread]; - L1Vector<fvec>& timeR = fTripletHitR_Time[Thread]; - L1Vector<fvec>& timeER = fTripletHitR_TimeErr[Thread]; - L1Vector<fvec>& du3 = fTripletHitR_dUfront[Thread]; - L1Vector<fvec>& dv3 = fTripletHitR_dUback[Thread]; + L1Vector<fvec>& t_3 = fTripletHitR_Time[Thread]; + L1Vector<fvec>& dt2_3 = fTripletHitR_TimeErr[Thread]; + L1Vector<fvec>& du2_3 = fTripletHitR_dUfront[Thread]; + L1Vector<fvec>& dv2_3 = fTripletHitR_dUback[Thread]; T_3.clear(); hitsl_3.clear(); @@ -1552,10 +1546,10 @@ inline void L1Algo::CreatePortionOfTriplets( u_front3.clear(); u_back3.clear(); z_pos3.clear(); - du3.clear(); - dv3.clear(); - timeR.clear(); - timeER.clear(); + du2_3.clear(); + dv2_3.clear(); + t_3.clear(); + dt2_3.clear(); /// Find the triplets(right hit). Reformat data in the portion of triplets. @@ -1567,7 +1561,7 @@ inline void L1Algo::CreatePortionOfTriplets( n_2, hitsm_2, i1_2, // output - n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER); + n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, du2_3, dv2_3, t_3, dt2_3); for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i) @@ -1582,7 +1576,7 @@ inline void L1Algo::CreatePortionOfTriplets( Tindex n3_V = (n3 + fvec::size() - 1) / fvec::size(); findTripletsStep1( // input - n3_V, star, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER, + n3_V, star, u_front3, u_back3, z_pos3, du2_3, dv2_3, t_3, dt2_3, // output T_3); diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h index f9975ef476f6d321304bcbb816e875d86546906d..9797c173d8514307412d965d3b34ddb29f22b2f2 100644 --- a/reco/L1/L1Algo/L1Filtration.h +++ b/reco/L1/L1Algo/L1Filtration.h @@ -14,7 +14,7 @@ #define cnst const fvec -inline void FilterTime(L1TrackPar& T, fvec t, fvec dt, fvec timeInfo = fvec::One(), fvec w = fvec::One()) +inline void FilterTime(L1TrackPar& T, fvec t, fvec dt2, fvec timeInfo = fvec::One(), fvec w = fvec::One()) { // filter track with a time measurement @@ -30,8 +30,6 @@ inline void FilterTime(L1TrackPar& T, fvec t, fvec dt, fvec timeInfo = fvec::One w.setZero(timeInfo <= fvec::Zero()); - fvec dt2 = dt * dt; - // when dt0 is much smaller than current time error, // set track time exactly to the measurement value without filtering // it helps to keep the initial time errors reasonably small diff --git a/reco/L1/L1Algo/L1Hit.h b/reco/L1/L1Algo/L1Hit.h index acdea6e3c422bc8e15316e71e7a130dc15248cfe..983b0462a980e44b317772f5de7e7b1b32072ef4 100644 --- a/reco/L1/L1Algo/L1Hit.h +++ b/reco/L1/L1Algo/L1Hit.h @@ -38,9 +38,9 @@ public: float v = 0.f; ///< measured V coordinate [cm] float t = 0.f; ///< measured time [ns] float z = 0.f; ///< fixed Z coordinate [cm] - float du = 0.f; ///< measured uncertainty of U coordinate [cm] - float dv = 0.f; ///< measured uncertainty of V coordinate [cm] - float dt = 0.f; ///< measured uncertainty of time [ns] + float du2 = 0.f; ///< measured uncertainty of U coordinate [cm] + float dv2 = 0.f; ///< measured uncertainty of V coordinate [cm] + float dt2 = 0.f; ///< measured uncertainty of time [ns] int ID = 0; ///< index of hit before hits sorting int iSt = 0; ///< index of station in the active stations array // TODO: Test speed penalty of using iSt index @@ -56,9 +56,9 @@ private: ar& v; ar& t; ar& z; - ar& du; - ar& dv; - ar& dt; + ar& du2; + ar& dv2; + ar& dt2; ar& ID; ar& iSt; } diff --git a/reco/L1/L1Algo/L1HitPoint.h b/reco/L1/L1Algo/L1HitPoint.h index ceb27e5c918a7cf2826e2ac525471718a4ba82ca..5fd6b8111fbc615e7b29fcade898978017c811ae 100644 --- a/reco/L1/L1Algo/L1HitPoint.h +++ b/reco/L1/L1Algo/L1HitPoint.h @@ -5,140 +5,50 @@ #ifndef _L1HitPoint_h_ #define _L1HitPoint_h_ -/// contain strips positions and coordinates of hit -#if 1 +/// contain strips positions and coordinates of the hit + struct L1HitPoint { - L1HitPoint() - : // x(0.f) - // , y(0.f) - // , dx(0.f) - // , dy(0.f) - // , dxy(0.f) - u(0.f) - , v(0.f) - , du(0.f) - , dv(0.f) - , z(0.f) - , time(0.f) - , timeEr(2.9f) {}; - L1HitPoint( - // fscal dx_, - // fscal dy_, - // fscal dxy_, - fscal du_, fscal dv_, fscal z_, fscal v_, fscal u_, float time_, float /*timeEv1_*/ = 0, - float timeEr_ = 2.9f) - : // : x(x_) - // , y(y_) - // , dx(dx_) - // , dy(dy_) - // , dxy(dxy_) - u(u_) - , v(v_) - , du(du_) - , dv(dv_) - , z(z_) - , time(time_) - , timeEr(timeEr_) {}; - // L1HitPoint(fscal x_, fscal y_, fscal z_, fscal v_, fscal u_, fscal time_, unsigned short int n_ = 0): - // x(x_), y(y_), z(z_), u(u_), v(v_), time(time_){}; + L1HitPoint() = default; - // fscal Xs() const { return X() / Z(); } - // fscal Ys() const { return Y() / Z(); } // value to sort hits by + L1HitPoint(fscal z_, fscal v_, fscal u_, fscal t_, fscal du2_, fscal dv2_, fscal dt2_) + : z(z_) + , u(u_) + , v(v_) + , t(t_) + , du2(du2_) + , dv2(dv2_) + , dt2(dt2_) {}; - // fscal X() const { return x; } - // fscal Y() const { return y; } - // fscal dX() const { return dx; } - // fscal dY() const { return dy; } - // fscal dXY() const { return dxy; } - fscal dU() const { return du; } - fscal dV() const { return dv; } fscal Z() const { return z; } fscal U() const { return u; } fscal V() const { return v; } - // unsigned short int N() const { return n; } - // int GetSortIndex() const { return SortIndex; } - // void SetX(fscal X1) { x = X1; } - // void SetY(fscal Y1) { y = Y1; } - void SetZ(fscal Z1) { z = Z1; } - void SetU(fscal U1) { u = U1; } - void SetV(fscal V1) { v = V1; } - - void Set( - // const float& dx1, - // const float& dy1, - // const float& xy1, - const float& z1, const fscal& u1, const fscal& v1, const fscal& du1, const fscal& dv1, const float& time1, - float timeEr1) - { - // x = x1; - // y = y1; - // dx = dx1; - // dy = dy1; - du = du1; - // dxy = xy1; - dv = dv1; - z = z1; - u = u1; - v = v1; - time = time1; - timeEr = timeEr1; - } - -private: - // x, y, - float u, v, du, dv; - float z; // TODO: may be we should use iz - // x\u, v - front and back strips; x, y, z - coordinates of hits - -public: - float time, timeEr; - // int track; -}; -#else - -static const float R = 60; -static const float shortPackingConstant = 2 * R / 65535.f; -//TODO: change the Z conversion. Hit Z maybe negative in the new setup. -static const float MZ = 110; -static const float shortPackingConstantZ = MZ / 65535.f; + fscal T() const { return t; } -/// contain strips positions and coordinates of hit -struct L1HitPoint { - L1HitPoint() {}; - L1HitPoint(fscal x_, fscal y_, fscal z_, fscal v_, fscal u_, unsigned short int n_ = 0) - : x(f2s(x_)) - , y(f2s(y_)) - , z(f2sZ(z_)) - , u(f2s(u_)) - , v(f2s(v_)) - , n(n_) {}; + fscal dU2() const { return du2; } + fscal dV2() const { return dv2; } + fscal dT2() const { return dt2; } - fscal Xs() const { return X() / Z(); } - fscal Ys() const { return Y() / Z(); } // value to sort hits by + void SetZ(fscal z_) { z = z_; } + void SetU(fscal u_) { u = u_; } + void SetV(fscal v_) { v = v_; } + void SetT(fscal t_) { t = t_; } - fscal X() const { return s2f(x); } - fscal Y() const { return s2f(y); } - fscal Z() const { return s2fZ(z); } - fscal U() const { return s2f(u); } - fscal V() const { return s2f(v); } - - unsigned short int N() const { return n; } + void Set(const fscal& z_, const fscal& u_, const fscal& v_, const fscal& t_, const fscal& du2_, const fscal& dv2_, + const fscal& dt2_) + { + z = z_; + u = u_; + v = v_; + t = t_; + du2 = du2_; + dv2 = dv2_; + dt2 = dt2_; + } private: - //unsigned short int - unsigned short int f2s(float f) const { return (f + R) / shortPackingConstant; } - float s2f(unsigned short int f) const { return (float(f) + 0.5) * shortPackingConstant - R; } - - unsigned short int f2sZ(float f) const { return (f) / shortPackingConstantZ; } - float s2fZ(unsigned short int f) const { return (float(f) + 0.5) * shortPackingConstantZ; } - - unsigned short int x, y; - unsigned short int z; // TODO: may be we should use iz - unsigned short int u, - v; // x\u, v - front and back strips; x, y, z - coordinates of hits + fscal z {0.}, u {0.}, v {0.}, t {0.}, du2 {0.}, dv2 {0.}, dt2 {0.}; }; -#endif #endif diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h index 7e176b6ef04787d9290b56d334cbe7f181b17b2d..c0103b1b7ffafe5fdbd1900aafa1fb31cd9a6de0 100644 --- a/reco/L1/L1Algo/L1Station.h +++ b/reco/L1/L1Algo/L1Station.h @@ -95,7 +95,7 @@ public: /// \param dv V coordinate uncertainty [length unit] /// \return Covariance matrix of hit position in Cartesian coordinates: [dxx, dxy, dyy] [length unit squared] template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool> = true> - std::tuple<T, T, T> FormXYCovarianceMatrix(T du, T dv) const; + std::tuple<T, T, T> FormXYCovarianceMatrix(T du2, T dv2) const; } _fvecalignment; @@ -133,20 +133,19 @@ std::pair<T, T> L1Station::ConvUVtoXY(T u, T v) const // --------------------------------------------------------------------------------------------------------------------- // template<typename T, std::enable_if_t<std::is_floating_point<T>::value || std::is_same<T, fvec>::value, bool>> -std::tuple<T, T, T> L1Station::FormXYCovarianceMatrix(T du, T dv) const +std::tuple<T, T, T> L1Station::FormXYCovarianceMatrix(T du2, T dv2) const { if constexpr (std::is_same<T, fvec>::value) { - return std::make_tuple( - (xInfo.cos_phi * du) * (xInfo.cos_phi * du) + (xInfo.sin_phi * dv) * (xInfo.sin_phi * dv), // dx2 - (xInfo.cos_phi * du) * (yInfo.cos_phi * du) + (xInfo.sin_phi * dv) * (yInfo.sin_phi * dv), // dxy - (yInfo.cos_phi * du) * (yInfo.cos_phi * du) + (yInfo.sin_phi * dv) * (yInfo.sin_phi * dv) // dy2 + return std::make_tuple(xInfo.cos_phi * xInfo.cos_phi * du2 + xInfo.sin_phi * xInfo.sin_phi * dv2, // dx2 + xInfo.cos_phi * yInfo.cos_phi * du2 + xInfo.sin_phi * yInfo.sin_phi * dv2, // dxy + yInfo.cos_phi * yInfo.cos_phi * du2 + yInfo.sin_phi * yInfo.sin_phi * dv2 // dy2 ); } else { return std::make_tuple( - (xInfo.cos_phi[0] * du) * (xInfo.cos_phi[0] * du) + (xInfo.sin_phi[0] * dv) * (xInfo.sin_phi[0] * dv), // dx2 - (xInfo.cos_phi[0] * du) * (yInfo.cos_phi[0] * du) + (xInfo.sin_phi[0] * dv) * (yInfo.sin_phi[0] * dv), // dxy - (yInfo.cos_phi[0] * du) * (yInfo.cos_phi[0] * du) + (yInfo.sin_phi[0] * dv) * (yInfo.sin_phi[0] * dv) // dy2 + xInfo.cos_phi[0] * xInfo.cos_phi[0] * du2 + xInfo.sin_phi[0] * xInfo.sin_phi[0] * dv2, // dx2 + xInfo.cos_phi[0] * yInfo.cos_phi[0] * du2 + xInfo.sin_phi[0] * yInfo.sin_phi[0] * dv2, // dxy + yInfo.cos_phi[0] * yInfo.cos_phi[0] * du2 + yInfo.sin_phi[0] * yInfo.sin_phi[0] * dv2 // dy2 ); } } diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx index 8a9a07ccf30b9f404879bc5672209458f1bfc2e3..cfeb883c7405729fac8ccd734996c9a580072867 100644 --- a/reco/L1/L1Algo/L1TrackFitter.cxx +++ b/reco/L1/L1Algo/L1TrackFitter.cxx @@ -347,8 +347,8 @@ void L1Algo::L1KFTrackFitter() // NOTE: u- and v-axes are axes, orthogonal to front and back strips of the station, respectively. fvec u[L1Constants::size::kMaxNstations]; // Hit position along the u-axis [cm] fvec v[L1Constants::size::kMaxNstations]; // Hit position along the v-axis [cm] - fvec d_u[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the u-axis [cm] - fvec d_v[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the v-axis [cm] + fvec du2[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the u-axis [cm] + fvec dv2[L1Constants::size::kMaxNstations]; // Hit position uncertainty along the v-axis [cm] fvec x[L1Constants::size::kMaxNstations]; // Hit position along the x-axis [cm] fvec y[L1Constants::size::kMaxNstations]; // Hit position along the y-axis [cm] @@ -358,8 +358,8 @@ void L1Algo::L1KFTrackFitter() fvec z[L1Constants::size::kMaxNstations]; // Hit position along the z-axis (precised) [cm] - fvec time[L1Constants::size::kMaxNstations]; // Hit time [ns] - fvec timeEr[L1Constants::size::kMaxNstations]; // Hit time uncertainty [ns] + fvec time[L1Constants::size::kMaxNstations]; // Hit time [ns] + fvec dt2[L1Constants::size::kMaxNstations]; // Hit time uncertainty [ns] squared fvec x_first; fvec y_first; @@ -368,7 +368,7 @@ void L1Algo::L1KFTrackFitter() fvec d_xy_fst; fvec time_first; - fvec time_er_first; + fvec dt2_first; fvec x_last; fvec y_last; @@ -377,8 +377,8 @@ void L1Algo::L1KFTrackFitter() fvec d_xy_lst; fvec time_last; - fvec time_er_last; - // fvec time_er_lst; /// TODO: Why are there two different variables for the time error on the last station? + fvec dt2_last; + // fvec dt2_lst; /// TODO: Why are there two different variables for the time error on the last station? fvec Sy[L1Constants::size::kMaxNstations]; fvec w[L1Constants::size::kMaxNstations]; @@ -436,16 +436,16 @@ void L1Algo::L1KFTrackFitter() u[ista][iVec] = hit.u; v[ista][iVec] = hit.v; - d_u[ista][iVec] = hit.du; - d_v[ista][iVec] = hit.dv; + du2[ista][iVec] = hit.du2; + dv2[ista][iVec] = hit.dv2; std::tie(x_temp, y_temp) = sta[ista].ConvUVtoXY<fvec>(u[ista], v[ista]); x[ista][iVec] = x_temp[iVec]; y[ista][iVec] = y_temp[iVec]; time[ista][iVec] = hit.t; - timeEr[ista][iVec] = hit.dt; + dt2[ista][iVec] = hit.dt2; z[ista][iVec] = hit.z; sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp); - std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(d_u[ista], d_v[ista]); + std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]); fB[ista].x[iVec] = fB_temp.x[iVec]; fB[ista].y[iVec] = fB_temp.y[iVec]; @@ -455,7 +455,7 @@ void L1Algo::L1KFTrackFitter() x_first[iVec] = x[ista][iVec]; y_first[iVec] = y[ista][iVec]; time_first[iVec] = time[ista][iVec]; - time_er_first[iVec] = timeEr[ista][iVec]; + dt2_first[iVec] = dt2[ista][iVec]; d_xx_fst[iVec] = d_xx[ista][iVec]; d_yy_fst[iVec] = d_yy[ista][iVec]; d_xy_fst[iVec] = d_xy[ista][iVec]; @@ -471,7 +471,7 @@ void L1Algo::L1KFTrackFitter() d_yy_lst[iVec] = d_yy[ista][iVec]; d_xy_lst[iVec] = d_xy[ista][iVec]; time_last[iVec] = time[ista][iVec]; - time_er_last[iVec] = timeEr[ista][iVec]; + dt2_last[iVec] = dt2[ista][iVec]; staLast.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec]; staLast.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec]; staLast.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec]; @@ -505,10 +505,10 @@ void L1Algo::L1KFTrackFitter() int ista = nStations - 1; - 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.)); + time_last = iif(w_time[ista] > fvec::Zero(), time_last, fvec::Zero()); + dt2_last = iif(w_time[ista] > fvec::Zero(), dt2_last, fvec(100.)); - FilterFirst(fit, 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, dt2_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst); fldZ1 = z[ista]; @@ -551,9 +551,9 @@ void L1Algo::L1KFTrackFitter() fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr); } - fit.Filter(sta[ista].frontInfo, u[ista], d_u[ista] * d_u[ista], w1); - fit.Filter(sta[ista].backInfo, v[ista], d_v[ista] * d_v[ista], w1); - fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo); + fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1); + fit.Filter(sta[ista].backInfo, v[ista], dv2[ista], w1); + fit.FilterTime(time[ista], dt2[ista], w1_time, sta[ista].timeInfo); fldB2 = fldB1; fldZ2 = fldZ1; @@ -668,7 +668,7 @@ void L1Algo::L1KFTrackFitter() ista = 0; - FilterFirst(fit, 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, dt2_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst); qp01 = tr.qp; @@ -708,9 +708,9 @@ void L1Algo::L1KFTrackFitter() fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr); } - fit.Filter(sta[ista].frontInfo, u[ista], d_u[ista] * d_u[ista], w1); - fit.Filter(sta[ista].backInfo, v[ista], d_v[ista] * d_v[ista], w1); - fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo); + fit.Filter(sta[ista].frontInfo, u[ista], du2[ista], w1); + fit.Filter(sta[ista].backInfo, v[ista], dv2[ista], w1); + fit.FilterTime(time[ista], dt2[ista], w1_time, sta[ista].timeInfo); fldB2 = fldB1; fldZ2 = fldZ1; @@ -791,8 +791,8 @@ void L1Algo::L1KFTrackFitterMuch() fvec u[L1Constants::size::kMaxNstations]; fvec v[L1Constants::size::kMaxNstations]; - fvec d_u[L1Constants::size::kMaxNstations]; - fvec d_v[L1Constants::size::kMaxNstations]; + fvec du2[L1Constants::size::kMaxNstations]; + fvec dv2[L1Constants::size::kMaxNstations]; fvec x[L1Constants::size::kMaxNstations]; fvec y[L1Constants::size::kMaxNstations]; @@ -803,7 +803,7 @@ void L1Algo::L1KFTrackFitterMuch() fvec z[L1Constants::size::kMaxNstations]; fvec time[L1Constants::size::kMaxNstations]; - fvec timeEr[L1Constants::size::kMaxNstations]; + fvec dt2[L1Constants::size::kMaxNstations]; fvec x_first; fvec y_first; @@ -812,7 +812,7 @@ void L1Algo::L1KFTrackFitterMuch() fvec d_xy_fst; fvec time_first; - fvec time_er_fst; + fvec dt2_fst; fvec x_last; fvec y_last; @@ -821,7 +821,7 @@ void L1Algo::L1KFTrackFitterMuch() fvec d_xy_lst; fvec time_last; - fvec time_er_lst; + fvec dt2_lst; fvec dz; /// !!! fvec Sy[L1Constants::size::kMaxNstations]; @@ -883,10 +883,10 @@ void L1Algo::L1KFTrackFitterMuch() x[ista][iVec] = x_temp[iVec]; y[ista][iVec] = y_temp[iVec]; time[ista][iVec] = hit.t; - timeEr[ista][iVec] = hit.dt; - d_u[ista][iVec] = hit.du; - d_v[ista][iVec] = hit.dv; - std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(d_u[ista], d_v[ista]); + dt2[ista][iVec] = hit.dt2; + du2[ista][iVec] = hit.du2; + dv2[ista][iVec] = hit.dv2; + std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]); // mom[ista][iVec] = hit.p; z[ista][iVec] = hit.z; @@ -899,7 +899,7 @@ void L1Algo::L1KFTrackFitterMuch() x_first[iVec] = x[ista][iVec]; y_first[iVec] = y[ista][iVec]; time_first[iVec] = time[ista][iVec]; - time_er_fst[iVec] = timeEr[ista][iVec]; + dt2_fst[iVec] = dt2[ista][iVec]; d_xx_fst[iVec] = d_xx[ista][iVec]; d_yy_fst[iVec] = d_yy[ista][iVec]; d_xy_fst[iVec] = d_xy[ista][iVec]; @@ -912,7 +912,7 @@ void L1Algo::L1KFTrackFitterMuch() x_last[iVec] = x[ista][iVec]; y_last[iVec] = y[ista][iVec]; time_last[iVec] = time[ista][iVec]; - time_er_lst[iVec] = timeEr[ista][iVec]; + dt2_lst[iVec] = dt2[ista][iVec]; d_xx_lst[iVec] = d_xx[ista][iVec]; d_yy_lst[iVec] = d_yy[ista][iVec]; d_xy_lst[iVec] = d_xy[ista][iVec]; @@ -962,7 +962,7 @@ void L1Algo::L1KFTrackFitterMuch() FilterFirst(T, x_first, y_first, staFirst); - FilterFirst(T1, x_first, y_first, time_first, time_er_fst, staFirst, d_xx_fst, d_yy_fst, d_xy_fst); + FilterFirst(T1, x_first, y_first, time_first, dt2_fst, staFirst, d_xx_fst, d_yy_fst, d_xy_fst); fldZ1 = z[i]; @@ -1029,7 +1029,7 @@ void L1Algo::L1KFTrackFitterMuch() L1Filter(T, info, v[i], w1); T1.Filter(info, v[i], d_yy[i], w1); - T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo); + T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo); } if (i >= 8) { @@ -1107,7 +1107,7 @@ void L1Algo::L1KFTrackFitterMuch() info.sigma2 = d_yy[i]; L1Filter(T, info, v[i], w1); T1.Filter(info, v[i], d_yy[i], w1); - T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo); + T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo); } } // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 ); @@ -1156,7 +1156,7 @@ void L1Algo::L1KFTrackFitterMuch() FilterFirst(T, x_last, y_last, staLast); - FilterFirstL(T1, x_last, y_last, time_last, time_er_lst, staLast, d_xx_lst, d_yy_lst, d_xy_lst); + FilterFirstL(T1, x_last, y_last, time_last, dt2_lst, staLast, d_xx_lst, d_yy_lst, d_xy_lst); qp0 = T.qp; qp01 = tr.qp; @@ -1250,7 +1250,7 @@ void L1Algo::L1KFTrackFitterMuch() L1Filter(T, info, v[i], w1); T1.Filter(info, v[i], d_yy[i], w1); - T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo); + T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo); } if (i < fNfieldStations - 1) { @@ -1294,13 +1294,13 @@ void L1Algo::L1KFTrackFitterMuch() } L1UMeasurementInfo info = sta[i].frontInfo; - // info.sigma2 = d_u[i] * d_u[i]; + info.sigma2 = du2[i]; T1.Filter(info, u[i], info.sigma2, w1); - info = sta[i].backInfo; - // info.sigma2 = d_v[i] * d_v[i]; + info = sta[i].backInfo; + info.sigma2 = dv2[i]; T1.Filter(info, v[i], info.sigma2, w1); - T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo); + T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo); fldB2 = fldB1; @@ -1598,7 +1598,8 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Stat track.fTr.chi2 = ZERO; } -void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt, L1Station& st, fvec& /*d_xx*/, + +void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& /*d_xx*/, fvec& /*d_yy*/, fvec& /*d_xy*/) { track.fTr.C00 = st.XYInfo.C00; @@ -1621,7 +1622,7 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& track.fTr.C52 = ZERO; track.fTr.C53 = ZERO; track.fTr.C54 = ZERO; - track.fTr.C55 = dt * dt; + track.fTr.C55 = dt2; track.fTr.x = x; track.fTr.y = y; @@ -1630,7 +1631,8 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& track.fTr.chi2 = ZERO; } -void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt, L1Station& st) + +void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt2, L1Station& st) { track.fTr.C00 = st.XYInfo.C00; track.fTr.C10 = st.XYInfo.C10; @@ -1652,7 +1654,7 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fv track.fTr.C52 = ZERO; track.fTr.C53 = ZERO; track.fTr.C54 = ZERO; - track.fTr.C55 = dt * dt; + track.fTr.C55 = dt2; track.fTr.x = x; track.fTr.y = y; @@ -1660,8 +1662,7 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fv track.fTr.chi2 = ZERO; } - -void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt, L1Station& /*st*/, fvec& d_xx, +void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt2, L1Station& /*st*/, fvec& d_xx, fvec& d_yy, fvec& d_xy) { // initialize covariance matrix @@ -1688,7 +1689,7 @@ void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, f track.fTr.C52 = ZERO; track.fTr.C53 = ZERO; track.fTr.C54 = ZERO; - track.fTr.C55 = dt * dt; + track.fTr.C55 = dt2; track.fTr.x = x; track.fTr.y = y; diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx index fd310acdd4fbd008a1d90539b0a3aef7a7e769c4..974c6e10fff85db0c2f6425de26ca0246d7aa5ba 100644 --- a/reco/L1/L1Algo/L1TrackParFit.cxx +++ b/reco/L1/L1Algo/L1TrackParFit.cxx @@ -144,7 +144,7 @@ void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w) fTr.C55 -= K5 * F5; } -void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo) +void L1TrackParFit::FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo) { // filter track with a time measurement @@ -160,8 +160,6 @@ void L1TrackParFit::FilterTime(fvec t, fvec dt, fvec w, fvec timeInfo) w.setZero(timeInfo <= fvec::Zero()); - fvec dt2 = dt * dt; - // when dt0 is much smaller than current time error, // set track time exactly to the measurement value without filtering // it helps to keep the initial time errors reasonably small diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h index cb965a9b4712f59a0cb5145ccb9df74f5439f8e7..2a800cf265ecd33d9a5db23ea622ce77723e53c2 100644 --- a/reco/L1/L1Algo/L1TrackParFit.h +++ b/reco/L1/L1Algo/L1TrackParFit.h @@ -49,7 +49,7 @@ public: void Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& sigma2, const fvec& w); void FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y); - void FilterTime(fvec t0, fvec dt0, fvec w = 1., fvec timeInfo = 1.); + void FilterTime(fvec t, fvec dt2, fvec w, fvec timeInfo); void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w = 1.); void Extrapolate(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);