diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 0168f1d368b1f669c76e1132adc6e3dc505a8e47..e57ee75a653da7ffa49e822234029f989c62e34a 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -1,9 +1,9 @@
-#Create a library called "libL1" which includes the source files given in
-#the array.
-#The extension is already found.Any number of sources could be listed here.
+# Create a library called "libL1" which includes the source files given in
+# the array.
+# The extension is already found.Any number of sources could be listed here.
# extra warnings to examine the code
-#ADD_DEFINITIONS(-Wall -Wextra -Wshadow -Weffc++)
+# ADD_DEFINITIONS(-Wall -Wextra -Wshadow -Weffc++)
# L1 defines
ADD_DEFINITIONS(-DDO_TPCCATRACKER_EFF_PERFORMANCE -DNonhomogeneousField -DCBM -DUSE_TIMERS)
@@ -18,31 +18,35 @@ set(INCLUDE_DIRECTORIES
${CMAKE_CURRENT_SOURCE_DIR}/qa
${CMAKE_CURRENT_SOURCE_DIR}/L1Algo/utils
${CMAKE_CURRENT_SOURCE_DIR}/catools
- )
+)
set(SRCS
- #L1Algo / L1CATrackFinder.cxx
- #CbmL1Performance.cxx
- #CbmL1ReadEvent.cxx
- #CbmL1CATrdTrackFinderSA.cxx
- CbmL1TrdTracklet.cxx
- CbmL1TrdTracklet4.cxx
+
+ # L1Algo / L1CATrackFinder.cxx
+ # CbmL1Performance.cxx
+ # CbmL1ReadEvent.cxx
+ # CbmL1CATrdTrackFinderSA.cxx
+ CbmL1TrdTracklet.cxx
+ CbmL1TrdTracklet4.cxx
CbmL1.cxx
- #CbmL1TrdTrackFinderSts.cxx
- CbmL1TrackMerger.cxx
+
+ # CbmL1TrdTrackFinderSts.cxx
+ CbmL1TrackMerger.cxx
CbmL1TofMerger.cxx
- # L1AlgoInputData.cxx
- OffLineInterface/CbmL1StsTrackFinder.cxx
- OffLineInterface/CbmL1GlobalTrackFinder.cxx
+
+ # L1AlgoInputData.cxx
+ OffLineInterface/CbmL1StsTrackFinder.cxx
+ OffLineInterface/CbmL1GlobalTrackFinder.cxx
OffLineInterface/CbmL1GlobalFindTracksEvents.cxx
L1Algo/L1Algo.cxx
+ L1Algo/L1TrackPar.cxx
L1Algo/L1CATrackFinder.cxx
L1Algo/L1BranchExtender.cxx
L1Algo/L1TrackFitter.cxx
L1Algo/L1HitsSortHelper.cxx
L1Algo/L1Grid.cxx
- CbmL1Performance.cxx
+ CbmL1Performance.cxx
CbmL1ReadEvent.cxx
CbmCaPerformance.cxx
L1Algo/L1Station.cxx
@@ -105,41 +109,41 @@ set(HEADERS
qa/CbmCaInputQaBase.h
)
-
-
-
If(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
- ADD_DEFINITIONS(-Wall -Wsign-promo -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wno-parentheses) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
-#-- - Check for compiler flags
+ ADD_DEFINITIONS(-Wall -Wsign-promo -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wno-parentheses) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
+
+ # -- - Check for compiler flags
CHECK_CXX_COMPILER_FLAG("-Werror -Wno-pmf-conversions" HAS_PMF)
+
If(HAS_PMF)
ADD_DEFINITIONS(-Wno-pmf-conversions)
EndIf()
+
CHECK_CXX_COMPILER_FLAG("-Werror -Wstrict-null-sentinel" HAS_SENTINEL)
+
If(HAS_SENTINEL)
ADD_DEFINITIONS(-Wstrict-null-sentinel)
EndIf()
+
CHECK_CXX_COMPILER_FLAG("-Werror -Wno-non-template-friend" HAS_TEMPLATE_FRIEND)
+
If(HAS_TEMPLATE_FRIEND)
ADD_DEFINITIONS(-Wno-non-template-friend)
EndIf()
Else()
- ADD_DEFINITIONS(-Wall -Wsign-promo -Wno-pmf-conversions -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wstrict-null-sentinel -Wno-non-template-friend -Wno-parentheses) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
+ ADD_DEFINITIONS(-Wall -Wsign-promo -Wno-pmf-conversions -Wctor-dtor-privacy -Wreorder -Wno-deprecated -Wstrict-null-sentinel -Wno-non-template-friend -Wno-parentheses) # -Weffc++ -Wnon-virtual-dtor -Woverloaded-virtual -Wold-style-cast : wait for other parts of cbmroot\root.
EndIf()
-
-IF (SSE_FOUND)
+IF(SSE_FOUND)
ADD_DEFINITIONS(-DHAVE_SSE)
- SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS
- "-msse -O3")
+ SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS
+ "-msse -O3")
Message(STATUS "L1 will be compiled with SSE support")
-ELSE (SSE_FOUND)
+ELSE(SSE_FOUND)
Message(STATUS "L1 will be compiled without SSE support")
- SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS
- "-O3")
-ENDIF (SSE_FOUND)
-
-
+ SET_SOURCE_FILES_PROPERTIES(${SRCS} PROPERTIES COMPILE_FLAGS
+ "-O3")
+ENDIF(SSE_FOUND)
set(LIBRARY_NAME L1)
set(LINKDEF ${LIBRARY_NAME}LinkDef.h)
@@ -155,10 +159,10 @@ set(PUBLIC_DEPENDENCIES
ROOT::Graf
ROOT::Hist
ROOT::Physics
- )
+)
if(OpenMP_CXX_FOUND)
- list(APPEND PUBLIC_DEPENDENCIES OpenMP::OpenMP_CXX)
+ list(APPEND PUBLIC_DEPENDENCIES OpenMP::OpenMP_CXX)
endif()
set(PRIVATE_DEPENDENCIES
@@ -181,38 +185,37 @@ set(PRIVATE_DEPENDENCIES
ROOT::Matrix
ROOT::Minuit
ROOT::RIO
- )
-
+)
generate_cbm_library()
add_dependencies(G__L1 KFPARTICLE)
install(FILES CbmL1Counters.h
- L1Algo/L1Assert.h
- L1Algo/L1EventEfficiencies.h
- L1Algo/L1Branch.h
- L1Algo/L1HitPoint.h
- L1Algo/L1Hit.h
- L1Algo/L1Portion.h
- L1Algo/L1Triplet.h
- L1Algo/L1Event.h
- L1Algo/L1EventMatch.h
- L1Algo/L1ObjectInitController.h
- L1Algo/L1Constants.h
- L1Algo/L1Utils.h
- L1Algo/L1NaN.h
- L1Algo/L1SimdSerializer.h
- L1Algo/L1TrackPar.h
- L1Algo/L1Track.h
- qa/CbmCaInputQaBase.h
- DESTINATION include
- )
+ L1Algo/L1Assert.h
+ L1Algo/L1EventEfficiencies.h
+ L1Algo/L1Branch.h
+ L1Algo/L1HitPoint.h
+ L1Algo/L1Hit.h
+ L1Algo/L1Portion.h
+ L1Algo/L1Triplet.h
+ L1Algo/L1Event.h
+ L1Algo/L1EventMatch.h
+ L1Algo/L1ObjectInitController.h
+ L1Algo/L1Constants.h
+ L1Algo/L1Utils.h
+ L1Algo/L1NaN.h
+ L1Algo/L1SimdSerializer.h
+ L1Algo/L1TrackPar.h
+ L1Algo/L1Track.h
+ qa/CbmCaInputQaBase.h
+ DESTINATION include
+)
install(FILES L1Algo/L1Algo.h
- L1Algo/L1Branch.h
- L1Algo/L1Field.h
- L1Algo/L1Hit.h
- L1Algo/L1Vector.h
- DESTINATION include/L1Algo
- )
+ L1Algo/L1Branch.h
+ L1Algo/L1Field.h
+ L1Algo/L1Hit.h
+ L1Algo/L1Vector.h
+ DESTINATION include/L1Algo
+)
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 64beea36593134832fce9bfedaa2cbad142bc1e1..00d60be0de2bf0ec1f6a6dd3cd1e6de48d1348fd 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -70,7 +70,6 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
fvec dzi = fvec(1.) / (z1 - z0);
- const fvec vINF = fvec(.1);
T.x = x0;
T.y = y0;
if (initParams) {
@@ -80,23 +79,15 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
}
fit.SetQp0(qp0);
- T.z = z0;
- T.t = hit0.t;
+ T.z = z0;
+ T.t = hit0.t;
+ T.vi = 0.;
- // T.t[0]=(hit0.t+hit1.t+hit2.t)/3;
- T.chi2 = fvec(0.);
- T.NDF = fvec(2.);
+ T.ResetErrors(1., 1., .1, .1, 1., (sta0.timeInfo ? hit0.dt2 : 1.e6), 1.e6);
+ T.NDF = fvec(2.);
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;
- T.C40 = T.C41 = T.C42 = T.C43 = 0;
- T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = 0;
- T.C22 = T.C33 = vINF;
- T.C44 = 1.;
- T.C55 = sta0.timeInfo ? hit0.dt2 : vINF;
-
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
fvec fldZ0 = sta1.fZ; // suppose field is smoth
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index a35ab55006d2ba094bcef87c882e71c0137c5af2..0a9f7e8efb27f2d8e5a3cc7a54d9556dbd341db3 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -260,12 +260,6 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
L1TrackPar& T = fit.Tr();
- T.chi2 = 0.;
- T.NDF = (fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
-
- // TODO: iteration parameter: "Starting NDF of track parameters"
- // NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
-
T.x = xl;
T.y = yl;
T.z = zl;
@@ -273,16 +267,18 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
T.ty = ty;
T.qp = fvec(0.);
T.t = time;
+ T.vi = fvec(0.);
- T.C20 = T.C21 = fvec(0.);
- T.C30 = T.C31 = T.C32 = fvec(0.);
- T.C40 = T.C41 = T.C42 = T.C43 = fvec(0.);
- T.C50 = T.C51 = T.C52 = T.C53 = T.C54 = fvec(0.);
+ fvec txErr2 = fMaxSlopePV * fMaxSlopePV / fvec(9.);
+ fvec qpErr2 = fMaxInvMom * fMaxInvMom / fvec(9.);
- T.C22 = T.C33 = fMaxSlopePV * fMaxSlopePV / fvec(9.);
+ T.ResetErrors(1., 1., txErr2, txErr2, qpErr2, timeEr2, 1.e2);
+
+ T.chi2 = 0.;
+ T.NDF = (fpCurrentIteration->GetPrimaryFlag()) ? fvec(2.) : fvec(0.);
- T.C44 = fMaxInvMom / fvec(3.) * fMaxInvMom / fvec(3.);
- T.C55 = timeEr2;
+ // TODO: iteration parameter: "Starting NDF of track parameters"
+ // NDF = number of track parameters (6: x, y, tx, ty, qp, time) - number of measured parameters (3: x, y, time) on station or (2: x, y) on target
{ // add the target constraint
@@ -609,7 +605,7 @@ inline void L1Algo::findTripletsStep0( // input
// add the middle hit
- fit.ExtrapolateLine(zPos_2);
+ fit.ExtrapolateLineNoField(zPos_2);
fit.Filter(stam.frontInfo, u_front_2, du2_2);
fit.Filter(stam.backInfo, u_back_2, dv2_2);
@@ -686,14 +682,17 @@ inline void L1Algo::findTripletsStep0( // input
fit3.SetTrack(T2);
L1TrackPar& T_cur = fit3.Tr();
- fit3.ExtrapolateLine(zr);
+ fit3.ExtrapolateLineNoField(zr);
- if ((star.timeInfo) && (stam.timeInfo))
- if (fabs(T_cur.t[i2_4] - hitr.T()) > sqrt(T_cur.C55[i2_4] + hitr.dT2()) * 5) continue;
+ if ((star.timeInfo) && (stam.timeInfo)) {
+ auto dt = T_cur.t[i2_4] - hitr.T();
+ if (dt > (T_cur.C55[i2_4] + hitr.dT2()) * 5) continue;
+ }
// TODO: SG: hardcoded cut of 30 ns
- if ((star.timeInfo) && (stam.timeInfo))
+ if ((star.timeInfo) && (stam.timeInfo)) {
if (fabs(T_cur.t[i2_4] - hitr.T()) > 30) continue;
+ }
// - check whether hit belong to the window ( track position +\- errors ) -
// check lower boundary
@@ -709,6 +708,7 @@ inline void L1Algo::findTripletsStep0( // input
const fscal dY2 = dY * dY;
if (dY2 > dy_est2) continue; // if (yr > y_plus_new [i2_4] ) continue;
+
// check x-boundaries
fvec x, C00;
fit3.ExtrapolateXC00Line(zr, x, C00);
@@ -810,7 +810,7 @@ inline void L1Algo::findTripletsStep1( // input
fit.SetMask(fmask::One());
for (Tindex i3_V = 0; i3_V < n3_V; ++i3_V) {
fit.SetTrack(T_3[i3_V]);
- fit.ExtrapolateLine(z_Pos[i3_V]);
+ fit.ExtrapolateLineNoField(z_Pos[i3_V]);
fit.Filter(star.frontInfo, u_front_[i3_V], du2_3[i3_V]);
fit.Filter(star.backInfo, u_back_[i3_V], dv2_3[i3_V]);
if (kMcbm != fTrackingMode) { fit.FilterTime(t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
@@ -917,6 +917,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.tx = (x[1] - x[0]) * dz01;
T.ty = (y[1] - y[0]) * dz01;
T.qp = 0.; //(fscal) T3.qp[i3_4];
+ T.vi = 0.;
}
// repeat several times in order to improve the precision
@@ -927,21 +928,18 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
fit.Qp0()(fit.Qp0() > GetMaxInvMom()) = GetMaxInvMom();
fit.Qp0()(fit.Qp0() < -GetMaxInvMom()) = -GetMaxInvMom();
- T.ResetErrors(200., 200., 1., 1., 100., 1.e4);
- //T.ResetErrors(200., 200., 10., 10., 100., 1.e4);
- T.NDF = ndf;
T.qp = 0.;
+ T.vi = 0.;
int ih0 = 0;
T.x = x[ih0];
T.y = y[ih0];
T.z = z[ih0];
T.t = t[ih0];
- //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
+ T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
+ std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
- fit.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0]);
- fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0]);
- fit.FilterTime(t[ih0], dt2[ih0], sta[ih0].timeInfo);
+ T.NDF = ndf;
// add the target constraint
fit.AddTargetToLine(fTargX, fTargY, fTargZ, TargetXYInfo, fldTarget);
@@ -965,20 +963,19 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
fit.Qp0()(fit.Qp0() > GetMaxInvMom()) = GetMaxInvMom();
fit.Qp0()(fit.Qp0() < -GetMaxInvMom()) = -GetMaxInvMom();
- T.ResetErrors(200., 200., 1., 1., 100., 1.e4);
T.NDF = ndf;
T.qp = 0.;
+ T.vi = 0.;
int ih0 = NHits - 1;
T.x = x[ih0];
T.y = y[ih0];
T.z = z[ih0];
T.t = t[ih0];
- T.C55 = dt2[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]);
- fit.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0]);
- fit.FilterTime(t[ih0], dt2[ih0], sta[ih0].timeInfo);
+ T.ResetErrors(1., 1., 1., 1., 100., (sta[ih0].timeInfo ? dt2[ih0] : 1.e6), 1.e2);
+ std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du2[ih0], dv2[ih0]);
+
+ T.NDF = ndf;
for (int ih = NHits - 2; ih >= 0; --ih) {
fit.Extrapolate(z[ih], fld);
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index 80a1f062b31b78be6980cd451a39f6059d514fd6..16a469c6d2615c02c57f08c9bfe70a096d718689 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -106,68 +106,12 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
unsigned short stab = firstStation[iTr];
- Tb.x = extTracks[iTr].TFirst[0];
- Tb.y = extTracks[iTr].TFirst[1];
- Tb.tx = extTracks[iTr].TFirst[2];
- Tb.ty = extTracks[iTr].TFirst[3];
- Tb.qp = extTracks[iTr].TFirst[4];
- Tb.z = extTracks[iTr].TFirst[5];
- Tb.t = extTracks[iTr].TFirst[6];
- Tb.C00 = extTracks[iTr].CFirst[0];
- Tb.C10 = extTracks[iTr].CFirst[1];
- Tb.C11 = extTracks[iTr].CFirst[2];
- Tb.C20 = extTracks[iTr].CFirst[3];
- Tb.C21 = extTracks[iTr].CFirst[4];
- Tb.C22 = extTracks[iTr].CFirst[5];
- Tb.C30 = extTracks[iTr].CFirst[6];
- Tb.C31 = extTracks[iTr].CFirst[7];
- Tb.C32 = extTracks[iTr].CFirst[8];
- Tb.C33 = extTracks[iTr].CFirst[9];
- Tb.C40 = extTracks[iTr].CFirst[10];
- Tb.C41 = extTracks[iTr].CFirst[11];
- Tb.C42 = extTracks[iTr].CFirst[12];
- Tb.C43 = extTracks[iTr].CFirst[13];
- Tb.C44 = extTracks[iTr].CFirst[14];
- Tb.C50 = extTracks[iTr].CFirst[15];
- Tb.C51 = extTracks[iTr].CFirst[16];
- Tb.C52 = extTracks[iTr].CFirst[17];
- Tb.C53 = extTracks[iTr].CFirst[18];
- Tb.C54 = extTracks[iTr].CFirst[19];
- Tb.C55 = extTracks[iTr].CFirst[20];
-
+ Tb.copyFromArrays(extTracks[iTr].TFirst, extTracks[iTr].CFirst);
fitB.SetQp0(fitB.Tr().qp);
unsigned short staf = lastStation[jTr];
- Tf.x = extTracks[jTr].TLast[0];
- Tf.y = extTracks[jTr].TLast[1];
- Tf.tx = extTracks[jTr].TLast[2];
- Tf.ty = extTracks[jTr].TLast[3];
- Tf.qp = extTracks[jTr].TLast[4];
- Tf.z = extTracks[jTr].TLast[5];
- Tf.t = extTracks[jTr].TLast[6];
- Tf.C00 = extTracks[jTr].CLast[0];
- Tf.C10 = extTracks[jTr].CLast[1];
- Tf.C11 = extTracks[jTr].CLast[2];
- Tf.C20 = extTracks[jTr].CLast[3];
- Tf.C21 = extTracks[jTr].CLast[4];
- Tf.C22 = extTracks[jTr].CLast[5];
- Tf.C30 = extTracks[jTr].CLast[6];
- Tf.C31 = extTracks[jTr].CLast[7];
- Tf.C32 = extTracks[jTr].CLast[8];
- Tf.C33 = extTracks[jTr].CLast[9];
- Tf.C40 = extTracks[jTr].CLast[10];
- Tf.C41 = extTracks[jTr].CLast[11];
- Tf.C42 = extTracks[jTr].CLast[12];
- Tf.C43 = extTracks[jTr].CLast[13];
- Tf.C44 = extTracks[jTr].CLast[14];
- Tf.C50 = extTracks[jTr].CLast[15];
- Tf.C51 = extTracks[jTr].CLast[16];
- Tf.C52 = extTracks[jTr].CLast[17];
- Tf.C53 = extTracks[jTr].CLast[18];
- Tf.C54 = extTracks[jTr].CLast[19];
- Tf.C55 = extTracks[jTr].CLast[20];
-
+ Tf.copyFromArrays(extTracks[jTr].TLast, extTracks[jTr].CLast);
fitF.SetQp0(fitF.Tr().qp);
if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index 7ea106f6e0904a2398faf80fa1ab8cb4f467b357..1b17d6982fad5f290204f703ea3718fbfe2a76c5 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -9,8 +9,8 @@
void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
{
fvec zeta, HCH;
- fvec F0, F1, F2, F3, F4, F5;
- fvec K1, K2, K3, K4, K5;
+ fvec F0, F1, F2, F3, F4, F5, F6;
+ fvec K1, K2, K3, K4, K5, K6;
zeta = info.cos_phi * fTr.x + info.sin_phi * fTr.y - u;
@@ -24,6 +24,7 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
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;
+ F6 = info.cos_phi * fTr.C60 + info.sin_phi * fTr.C61;
const fmask maskDoFilter = (HCH < sigma2 * 16.f);
//cnst maskDoFilter = _f32vec4_true;
@@ -43,6 +44,7 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
K3 = F3 * wi;
K4 = F4 * wi;
K5 = F5 * wi;
+ K6 = F6 * wi;
fTr.x -= F0 * zetawi;
fTr.y -= F1 * zetawi;
@@ -50,6 +52,7 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
fTr.ty -= F3 * zetawi;
fTr.qp -= F4 * zetawi;
fTr.t -= F5 * zetawi;
+ fTr.vi -= F6 * zetawi;
fTr.C00 -= F0 * F0 * wi;
fTr.C10 -= K1 * F0;
@@ -72,6 +75,13 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
fTr.C53 -= K5 * F3;
fTr.C54 -= K5 * F4;
fTr.C55 -= K5 * F5;
+ fTr.C60 -= K6 * F0;
+ fTr.C61 -= K6 * F1;
+ fTr.C62 -= K6 * F2;
+ fTr.C63 -= K6 * F3;
+ fTr.C64 -= K6 * F4;
+ fTr.C65 -= K6 * F5;
+ fTr.C66 -= K6 * F6;
}
@@ -86,6 +96,7 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
fvec F3 = fTr.C53;
fvec F4 = fTr.C54;
fvec F5 = fTr.C55;
+ fvec F6 = fTr.C65;
fvec HCH = fTr.C55;
@@ -112,6 +123,7 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
fvec K3 = F3 * wi;
fvec K4 = F4 * wi;
fvec K5 = F5 * wi;
+ fvec K6 = F6 * wi;
fTr.x -= F0 * zetawi;
fTr.y -= F1 * zetawi;
@@ -119,6 +131,7 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
fTr.ty -= F3 * zetawi;
fTr.qp -= F4 * zetawi;
fTr.t -= F5 * zetawi;
+ fTr.vi -= F6 * zetawi;
fTr.C00 -= F0 * F0 * wi;
fTr.C10 -= K1 * F0;
@@ -141,6 +154,13 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
fTr.C53 -= K5 * F3;
fTr.C54 -= K5 * F4;
fTr.C55 -= K5 * F5;
+ fTr.C60 -= K6 * F0;
+ fTr.C61 -= K6 * F1;
+ fTr.C62 -= K6 * F2;
+ fTr.C63 -= K6 * F3;
+ fTr.C64 -= K6 * F4;
+ fTr.C65 -= K6 * F5;
+ fTr.C66 -= K6 * F6;
}
@@ -149,8 +169,10 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
cnst TWO(2.);
fvec zeta0, zeta1, S00, S10, S11, si;
- 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;
+ fvec F00, F10, F20, F30, F40, F50, F60;
+ fvec F01, F11, F21, F31, F41, F51, F61;
+ fvec K00, K10, K20, K30, K40, K50, K60;
+ fvec K01, K11, K21, K31, K41, K51, K61;
zeta0 = fTr.x - x;
zeta1 = fTr.y - y;
@@ -162,12 +184,15 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
F30 = fTr.C30;
F40 = fTr.C40;
F50 = fTr.C50;
+ F60 = fTr.C60;
+
F01 = fTr.C10;
F11 = fTr.C11;
F21 = fTr.C21;
F31 = fTr.C31;
F41 = fTr.C41;
F51 = fTr.C51;
+ F61 = fTr.C61;
S00 = F00 + info.C00;
S10 = F10 + info.C10;
@@ -194,34 +219,51 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
K41 = F40 * S10 + F41 * S11;
K50 = F50 * S00 + F51 * S10;
K51 = F50 * S10 + F51 * S11;
+ K60 = F60 * S00 + F61 * S10;
+ K61 = F60 * S10 + F61 * S11;
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.t -= K50 * zeta0 + K51 * zeta1;
+ fTr.vi -= K60 * zeta0 + K61 * 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;
+
+ fTr.C60 -= K60 * F00 + K61 * F01;
+ fTr.C61 -= K60 * F10 + K61 * F11;
+ fTr.C62 -= K60 * F20 + K61 * F21;
+ fTr.C63 -= K60 * F30 + K61 * F31;
+ fTr.C64 -= K60 * F40 + K61 * F41;
+ fTr.C65 -= K60 * F50 + K61 * F51;
+ fTr.C66 -= K60 * F60 + K61 * F61;
}
void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY,
@@ -315,10 +357,313 @@ void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian f
while (!(fMaskF * 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, F);
+ ExtrapolateStepFull(zNew, F);
}
}
+void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
+ (cnst& zOut, // extrapolate to this z position
+ const L1FieldRegion& F)
+{
+ // use Q/p linearisation at fQp0
+ // implementation of the Runge-Kutta method without optimization
+ //
+
+ //
+ // Forth-order Runge-Kutta method for solution of the equation
+ // of motion of a particle with parameter qp = Q /P
+ // in the magnetic field B()
+ //
+ // ( x ) tx
+ // ( y ) ty
+ // ( tx) c_light * qp * L * ( tx*ty * Bx - (1+tx*tx) * By + ty * Bz )
+ // d ( ty) / dz = c_light * qp * L * ( (1+ty*ty) * Bx - tx*ty * By - tx * Bz ) ,
+ // ( qp) 0.
+ // ( t ) L * vi
+ // ( vi) 0.
+ //
+ // where L = sqrt ( 1 + tx*tx + ty*ty ) .
+ // c_light = 0.000299792458 [(GeV/c)/kG/cm]
+ // c_light_ns = 29.9792458 [cm/ns]
+ //
+ // In the following for RK step :
+ // r[7] = {x, y, tx, ty, qp, t, vi}
+ // dr(z)/dz = f(z,r)
+ //
+ //
+ //========================================================================
+ //
+ // NIM A395 (1997) 169-184; NIM A426 (1999) 268-282.
+ //
+ // the routine is based on LHC(b) utility code
+ //
+ //========================================================================
+
+
+ cnst c_light = 0.000299792458;
+
+ //----------------------------------------------------------------
+
+ cnst zMasked = iif(fMask, zOut, fTr.z);
+
+ fvec qp_in = fTr.qp;
+ cnst h = (zMasked - fTr.z);
+
+ cnst a[4] = {0., h * fvec(0.5), h * fvec(0.5), h};
+ cnst c[4] = {h / fvec(6.), h / fvec(3.), h / fvec(3.), h / fvec(6.)};
+
+ fvec f0[4] = {{0.}}; // ( dx/dz )[step]
+ fvec f1[4] = {{0.}}; // ( dy/dz )[step]
+ fvec f2[4] = {{0.}}; // ( dtx/dz )[step]
+ fvec f3[4] = {{0.}}; // ( dty/dz )[step]
+ fvec f4[4] = {{0.}}; // ( dqp/dz )[step]
+ fvec f5[4] = {{0.}}; // ( dt/dz )[step]
+ fvec f6[4] = {{0.}}; // ( dvi/dz )[step]
+
+ //fvec df0[4][7] = {{0.}}; // df0[step] / d {x,y,tx,ty,qp,t,vi}
+ //fvec df1[4][7] = {{0.}}; // df1[step] / d {x,y,tx,ty,qp,t,vi}
+ fvec df2[4][7] = {{0.}}; // df2[step] / d {x,y,tx,ty,qp,t,vi}
+ fvec df3[4][7] = {{0.}}; // df3[step] / d {x,y,tx,ty,qp,t,vi}
+ //fvec df4[4][7] = {{0.}}; // df4[step] / d {x,y,tx,ty,qp,t,vi}
+ fvec df5[4][7] = {{0.}}; // df5[step] / d {x,y,tx,ty,qp,t,vi}
+ //fvec df6[4][7] = {{0.}}; // df6[step] / d {x,y,tx,ty,qp,t,vi}
+
+ fvec k[5][7] = {{0.}}; // [ step-1 ] [ {x,y,tx,ty,qp,t,vi} ]
+
+ // Runge-Kutta steps for track parameters
+ //
+ {
+ fvec r0[7] = {fTr.x, fTr.y, fTr.tx, fTr.ty, fTr.qp, fTr.t, fTr.vi};
+
+ for (int step = 0; step < 4; ++step) {
+ fvec z = fTr.z + a[step];
+ fvec r[7];
+ for (int i = 0; i < 7; ++i) {
+ r[i] = r0[i] + a[step] * k[step][i];
+ }
+
+ fvec B[3];
+ cnst& Bx = B[0];
+ cnst& By = B[1];
+ cnst& Bz = B[2];
+
+ F.Get(r[0], r[1], z, B);
+
+ fvec tx = r[2];
+ fvec ty = r[3];
+ fvec tx2 = tx * tx;
+ fvec ty2 = ty * ty;
+ fvec txty = tx * ty;
+ fvec L2 = fvec(1.) + tx2 + ty2;
+ fvec L2i = fvec(1.) / L2;
+ fvec L = sqrt(L2);
+ fvec cL = c_light * L;
+ fvec cLqp0 = cL * fQp0;
+
+ f0[step] = tx;
+ f1[step] = ty;
+
+ fvec f2tmp = (txty * Bx + ty * Bz) - (fvec(1.) + tx2) * By;
+ f2[step] = cLqp0 * f2tmp;
+ df2[step][2] = cLqp0 * (tx * f2tmp * L2i + ty * Bx - fvec(2.) * tx * By);
+ df2[step][3] = cLqp0 * (ty * f2tmp * L2i + tx * Bx + Bz);
+ df2[step][4] = cL * f2tmp;
+
+ fvec f3tmp = -txty * By - tx * Bz + (fvec(1.) + ty2) * Bx;
+ f3[step] = cLqp0 * f3tmp;
+ df3[step][2] = cLqp0 * (tx * f3tmp * L2i - ty * By - Bz);
+ df3[step][3] = cLqp0 * (ty * f3tmp * L2i + fvec(2.) * ty * Bx - tx * By);
+ df3[step][4] = cL * f3tmp;
+
+ fvec vi = r[6];
+ fvec m2 = fMass2;
+ if (!fDoFitVelocity) { vi = sqrt(fvec(1.) + m2 * fQp0 * fQp0) / fvec(29.9792458f); }
+
+ f5[step] = vi * L;
+ df5[step][2] = vi * tx / L;
+ df5[step][3] = vi * ty / L;
+ df5[step][4] = 0.;
+ df5[step][5] = 0.;
+ df5[step][6] = L;
+
+ if (!fDoFitVelocity) {
+ df5[step][4] = (m2 * fQp0) * (L / sqrt(fvec(1.) + m2 * fQp0 * fQp0) / fvec(29.9792458f));
+ df5[step][6] = 0.;
+ }
+
+ k[step + 1][0] = f0[step];
+ k[step + 1][1] = f1[step];
+ k[step + 1][2] = f2[step];
+ k[step + 1][3] = f3[step];
+ k[step + 1][4] = f4[step]; // == 0
+ k[step + 1][5] = f5[step];
+ k[step + 1][6] = f6[step]; // == 0
+
+ } // end of Runge-Kutta step
+
+ fTr.x = r0[0] + (c[0] * k[1][0] + c[1] * k[2][0] + c[2] * k[3][0] + c[3] * k[4][0]);
+ fTr.y = r0[1] + (c[0] * k[1][1] + c[1] * k[2][1] + c[2] * k[3][1] + c[3] * k[4][1]);
+ fTr.tx = r0[2] + (c[0] * k[1][2] + c[1] * k[2][2] + c[2] * k[3][2] + c[3] * k[4][2]);
+ fTr.ty = r0[3] + (c[0] * k[1][3] + c[1] * k[2][3] + c[2] * k[3][3] + c[3] * k[4][3]);
+ //fTr.qp = r0[4] + (c[0] * k[1][4] + c[1] * k[2][4] + c[2] * k[3][4] + c[3] * k[4][4]);
+ fTr.t = r0[5] + (c[0] * k[1][5] + c[1] * k[2][5] + c[2] * k[3][5] + c[3] * k[4][5]);
+ //fTr.vi = r0[6] + (c[0] * k[1][6] + c[1] * k[2][6] + c[2] * k[3][6] + c[3] * k[4][6]);
+
+ fTr.z = zMasked;
+ }
+
+
+ // Jacobian of extrapolation
+
+ //
+ // derivatives d/dx and d/dy
+ //
+ fvec Jx[7] = {1., 0., 0., 0., 0., 0., 0.}; // D new { x,y,tx,ty,qp,t,vi } / D old x
+ fvec Jy[7] = {0., 1., 0., 0., 0., 0., 0.}; // D new { x,y,tx,ty,qp,t,vi } / D old y
+
+ //
+ // Runge-Kutta steps for derivatives d/dtx
+ //
+ fvec Jtx[7] = {0., 0., 1., 0., 0., 0., 0.}; // D new { x,y,tx,ty,qp,t,vi } / D old tx
+ {
+ for (int step = 0; step < 4; ++step) {
+ fvec r2 = Jtx[2] + a[step] * k[step][2]; // dtx[step]/dtx0
+ fvec r3 = Jtx[3] + a[step] * k[step][3]; // dty[step]/dtx0
+ k[step + 1][0] = r2;
+ k[step + 1][1] = r3;
+ k[step + 1][2] = df2[step][2] * r2 + df2[step][3] * r3;
+ k[step + 1][3] = df3[step][2] * r2 + df3[step][3] * r3;
+ k[step + 1][4] = 0.;
+ k[step + 1][5] = df5[step][2] * r2 + df5[step][3] * r3;
+ k[step + 1][6] = 0.;
+ }
+ for (int i = 0; i < 7; ++i) {
+ Jtx[i] += c[0] * k[1][i] + c[1] * k[2][i] + c[2] * k[3][i] + c[3] * k[4][i];
+ }
+ }
+
+ //
+ // Runge-Kutta steps for derivatives d/dty
+ //
+ fvec Jty[7] = {0., 0., 0., 1., 0., 0., 0.}; // D new { x,y,tx,ty,qp,t,vi } / D old ty
+ {
+ for (int step = 0; step < 4; ++step) {
+ fvec r2 = Jty[2] + a[step] * k[step][2]; // dtx[step]/dty0
+ fvec r3 = Jty[3] + a[step] * k[step][3]; // dty[step]/dty0
+ k[step + 1][0] = r2;
+ k[step + 1][1] = r3;
+ k[step + 1][2] = df2[step][2] * r2 + df2[step][3] * r3;
+ k[step + 1][3] = df3[step][2] * r2 + df3[step][3] * r3;
+ k[step + 1][4] = 0.;
+ k[step + 1][5] = df5[step][2] * r2 + df5[step][3] * r3;
+ k[step + 1][6] = 0.;
+ }
+ for (int i = 0; i < 7; ++i) {
+ Jty[i] += c[0] * k[1][i] + c[1] * k[2][i] + c[2] * k[3][i] + c[3] * k[4][i];
+ }
+ }
+
+ //
+ // Runge-Kutta steps for derivatives d/dqp
+ //
+ fvec Jqp[7] = {0., 0., 0., 0., 1., 0., 0.}; // D new { x,y,tx,ty,qp,t,vi } / D old qp
+ {
+ for (int step = 0; step < 4; ++step) {
+ fvec r2 = Jqp[2] + a[step] * k[step][2]; // dtx/dqp
+ fvec r3 = Jqp[3] + a[step] * k[step][3]; // dty/dqp;
+ fvec r4 = Jqp[4] + a[step] * k[step][4]; // dqp/dqp == 1
+ k[step + 1][0] = r2;
+ k[step + 1][1] = r3;
+ k[step + 1][2] = df2[step][2] * r2 + df2[step][3] * r3 + df2[step][4] * r4;
+ k[step + 1][3] = df3[step][2] * r2 + df3[step][3] * r3 + df3[step][4] * r4;
+ k[step + 1][4] = 0.;
+ k[step + 1][5] = df5[step][2] * r2 + df5[step][3] * r3 + df5[step][4] * r4;
+ k[step + 1][6] = 0.;
+ }
+ for (int i = 0; i < 7; ++i) {
+ Jqp[i] += c[0] * k[1][i] + c[1] * k[2][i] + c[2] * k[3][i] + c[3] * k[4][i];
+ }
+ }
+
+
+ //
+ // derivatives d/dt
+ //
+ fvec Jt[7] = {0., 0., 0., 0., 0., 1., 0.}; // D new { x,y,tx,ty,qp,t,vi } / D old t
+
+ //
+ // derivatives d/dvi
+ //
+ fvec Jvi[7] = {0., 0., 0., 0., 0., 0., 1.}; // D new { x,y,tx,ty,qp,t,vi } / D old vi
+ {
+ for (int step = 0; step < 4; ++step) {
+ fvec r2 = Jvi[2] + a[step] * k[step][2]; // dtx/dvi
+ fvec r3 = Jvi[3] + a[step] * k[step][3]; // dty/dvi;
+ fvec r6 = Jvi[6] + a[step] * k[step][6]; //(dvi/dvi == 1)
+ k[step + 1][0] = r2;
+ k[step + 1][1] = r3;
+ k[step + 1][2] = df2[step][2] * r2 + df2[step][3] * r3;
+ k[step + 1][3] = df3[step][2] * r2 + df3[step][3] * r3;
+ k[step + 1][4] = 0.;
+ k[step + 1][5] = df5[step][2] * r2 + df5[step][3] * r3 + df5[step][6] * r6;
+ k[step + 1][6] = 0.;
+ }
+ for (int i = 0; i < 7; ++i) {
+ Jvi[i] += c[0] * k[1][i] + c[1] * k[2][i] + c[2] * k[3][i] + c[3] * k[4][i];
+ }
+ }
+
+
+ { // update parameters
+ fvec dqp = qp_in - fQp0;
+ fTr.x += Jqp[0] * dqp;
+ fTr.y += Jqp[1] * dqp;
+ fTr.tx += Jqp[2] * dqp;
+ fTr.ty += Jqp[3] * dqp;
+ fTr.t += Jqp[5] * dqp;
+ }
+
+ // covariance matrix transport
+
+ fvec J[7][7];
+ for (int i = 0; i < 7; i++) {
+ J[i][0] = Jx[i]; // d {x,y,tx,ty,qp,t,vi} / dx
+ J[i][1] = Jy[i]; // d * / dy
+ J[i][2] = Jtx[i]; // d * / dtx
+ J[i][3] = Jty[i]; // d * / dty
+ J[i][4] = Jqp[i]; // d * / dqp
+ J[i][5] = Jt[i]; // d * / dt
+ J[i][6] = Jvi[i]; // d * / dvi
+ }
+
+ fvec C[7][7];
+ for (int i = 0; i < 7; i++) {
+ for (int j = 0; j < 7; j++) {
+ C[i][j] = fTr.C(i, j);
+ }
+ }
+
+ fvec JC[7][7];
+ for (int i = 0; i < 7; i++) {
+ for (int j = 0; j < 7; j++) {
+ JC[i][j] = 0.;
+ for (int m = 0; m < 7; m++) {
+ JC[i][j] += J[i][m] * C[m][j];
+ }
+ }
+ }
+ for (int i = 0; i < 7; i++) {
+ for (int j = 0; j < 7; j++) {
+ fvec Cij = 0.;
+ for (int m = 0; m < 7; m++) {
+ Cij += JC[i][m] * J[j][m];
+ }
+ fTr.C(i, j) = Cij;
+ }
+ }
+}
+
+
void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& zOut, // extrapolate to this z position
const L1FieldRegion& F)
@@ -377,7 +722,7 @@ void L1Fit::ExtrapolateStep // extrapolates track parameters and returns jacobi
fvec f2_ty[4], f3_ty[4], f4_ty[4]; // df* / dty
fvec f2_qp[4], f3_qp[4], f4_qp[4]; // df* / dqp
- fvec k[5][5] {{0.}};
+ fvec k[5][5] = {0., 0., 0., 0., 0.};
// Runge-Kutta steps for track parameters
//
@@ -936,13 +1281,14 @@ void L1Fit::ExtrapolateLine(cnst& z_out, const L1FieldRegion& F)
fQp0 = qp0;
}
-void L1Fit::ExtrapolateLine(cnst& z_out)
+void L1Fit::ExtrapolateLineNoField(cnst& z_out)
{
- // extrapolate the track assuming fQp0 == 0
+ // extrapolate the track assuming no field
//
// it is a copy of a sequence two routines
// L1ExtrapolateTime() and L1ExtrapolateLine()
// TODO: make it right
+ //
cnst c_light(29.9792458);
fvec dz = (z_out - fTr.z);
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 4e95f91e7f5e14285390212a55dba257d8f81757..72fe29a978b094b3c60bd9b88987a0018b490d96 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -41,9 +41,9 @@ public:
}
template<typename T>
- void SetTrack(const T* tr, const T* C)
+ void SetTrack(const T p[7], const T c[28])
{
- fTr.Set(tr, C);
+ fTr.copyFromArrays(p, c);
fQp0 = fTr.qp;
}
@@ -59,6 +59,8 @@ public:
fvec& Qp0() { return fQp0; }
+ void SetDoFitVelocity(bool v) { fDoFitVelocity = v; }
+
void SetOneEntry(const int i0, const L1Fit& T1, const int i1);
void Print(int i = -1);
@@ -90,8 +92,9 @@ public:
void Extrapolate(cnst& z_out, const L1FieldRegion& F);
void ExtrapolateStep(cnst& z_out, const L1FieldRegion& F);
+ void ExtrapolateStepFull(cnst& z_out, const L1FieldRegion& F);
void ExtrapolateStepAnalytic(cnst& z_out, const L1FieldRegion& F);
- void ExtrapolateLine(cnst& z_out);
+ void ExtrapolateLineNoField(cnst& z_out);
void ExtrapolateLine(cnst& z_out, const L1FieldRegion& F);
void EnergyLossCorrection(cnst& radThick, cnst& upstreamDirection);
@@ -139,6 +142,9 @@ private:
fvec fPipeRadThick {7.87e-3f}; // 0.7 mm Aluminium // TODO:
fvec fTargetRadThick {3.73e-2f * 2}; // 250 mum Gold // TODO:
+
+ bool fDoFitVelocity {0}; // should the track velocity be fitted as an independent parameter
+
} _fvecalignment;
// =============================================================================================
diff --git a/reco/L1/L1Algo/L1Track.h b/reco/L1/L1Algo/L1Track.h
index 23cbafaf74706b8c5a2ff9861d3f1020b83a0167..e6223f1132fec881708f653cf076fb9c1de9491f 100644
--- a/reco/L1/L1Algo/L1Track.h
+++ b/reco/L1/L1Algo/L1Track.h
@@ -35,12 +35,12 @@ public:
unsigned char n; ///< TODO: ??
float Momentum {kNaN}; ///< TODO: ??
float fTrackTime {kNaN}; ///< Track time
- fscal TFirst[7] {kNaN}; ///< Track parameters on the first station
- fscal CFirst[21] {kNaN}; ///< Track parameter covariation matrix elements on the first station
- fscal TLast[7] {kNaN}; ///< Track parameters on the last station
- fscal CLast[21] {kNaN}; ///< Track parameter covariation matrix elements on the second station
- fscal Tpv[7] {kNaN}; ///< Track parameters in the primary vertex
- fscal Cpv[21] {kNaN}; ///< Track parameter covariation matrix elements in the primary vertex
+ fscal TFirst[8] {kNaN}; ///< Track parameters on the first station
+ fscal CFirst[28] {kNaN}; ///< Track parameter covariation matrix elements on the first station
+ fscal TLast[8] {kNaN}; ///< Track parameters on the last station
+ fscal CLast[28] {kNaN}; ///< Track parameter covariation matrix elements on the second station
+ fscal Tpv[8] {kNaN}; ///< Track parameters in the primary vertex
+ fscal Cpv[28] {kNaN}; ///< Track parameter covariation matrix elements in the primary vertex
fscal chi2 {kNaN}; ///< Track fit chi-square value
short int NDF; ///< Track fit NDF value
@@ -66,11 +66,5 @@ public:
static bool compare(const L1Track& a, const L1Track& b) { return (a.Cpv[20] <= b.Cpv[20]); }
};
-// #include "cmath"
-// bool operator==(const L1Track &other) const {
-// cout<<int(NHits)<<" NHits"<<endl;
-// if ((other.NHits==NHits)&&(fabs(other.Momentum-Momentum)<1.e-6)) return true;
-// else return false;
-// }
#endif // L1Track_H
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 1bc56c0b0d08a321b802526c59287c6341f46ea7..4ed3a4b1dc825420d582555e272ca5a630e63a69 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -17,9 +17,6 @@ using std::vector;
namespace NS_L1TrackFitter
{
const fvec c_light(0.000299792458), c_light_i(fvec(1.) / c_light);
- const fvec ZERO(0.);
- const fvec ONE(1.);
- const fvec vINF(0.1);
} // namespace NS_L1TrackFitter
using namespace NS_L1TrackFitter;
@@ -45,6 +42,7 @@ void L1Algo::L1KFTrackFitter()
L1TrackPar& tr = fit.Tr();
fit.SetParticleMass(GetDefaultParticleMass());
+ fit.SetDoFitVelocity(false);
L1Track* t[fvec::size()] {nullptr};
@@ -228,7 +226,13 @@ void L1Algo::L1KFTrackFitter()
time_last = iif(w_time[ista], time_last, fvec::Zero());
dt2_last = iif(w_time[ista], dt2_last, fvec(1.e6));
- FilterFirst(fit, x_last, y_last, time_last, dt2_last, d_xx_lst, d_yy_lst, d_xy_lst);
+ tr.ResetErrors(d_xx_lst, d_yy_lst, 0.1, 0.1, 1.0, dt2_last, 1.e2);
+ tr.C10 = d_xy_lst;
+ tr.x = x_last;
+ tr.y = y_last;
+ tr.t = time_last;
+ tr.NDF = -3.0;
+
fldZ1 = z[ista];
@@ -301,74 +305,17 @@ void L1Algo::L1KFTrackFitter()
}
}
- const L1TrackPar& Tf = fit.Tr();
+ L1TrackPar Tf = fit.Tr();
+ if (kGlobal == fTrackingMode) { Tf.qp = fitpv.Tr().qp; }
for (int iVec = 0; iVec < nTracks_SIMD; 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.Tr().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];
- t[iVec]->CFirst[2] = Tf.C11[iVec];
- t[iVec]->CFirst[3] = Tf.C20[iVec];
- t[iVec]->CFirst[4] = Tf.C21[iVec];
- t[iVec]->CFirst[5] = Tf.C22[iVec];
- t[iVec]->CFirst[6] = Tf.C30[iVec];
- t[iVec]->CFirst[7] = Tf.C31[iVec];
- t[iVec]->CFirst[8] = Tf.C32[iVec];
- t[iVec]->CFirst[9] = Tf.C33[iVec];
- t[iVec]->CFirst[10] = Tf.C40[iVec];
- t[iVec]->CFirst[11] = Tf.C41[iVec];
- t[iVec]->CFirst[12] = Tf.C42[iVec];
- t[iVec]->CFirst[13] = Tf.C43[iVec];
- t[iVec]->CFirst[14] = Tf.C44[iVec];
- t[iVec]->CFirst[15] = Tf.C50[iVec];
- t[iVec]->CFirst[16] = Tf.C51[iVec];
- t[iVec]->CFirst[17] = Tf.C52[iVec];
- t[iVec]->CFirst[18] = Tf.C53[iVec];
- t[iVec]->CFirst[19] = Tf.C54[iVec];
- t[iVec]->CFirst[20] = Tf.C55[iVec];
-
+ Tf.copyToArrays(iVec, t[iVec]->TFirst, t[iVec]->CFirst);
t[iVec]->chi2 = Tf.chi2[iVec];
t[iVec]->NDF = (int) Tf.NDF[iVec];
}
for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec]->Tpv[0] = fitpv.Tr().x[iVec];
- t[iVec]->Tpv[1] = fitpv.Tr().y[iVec];
- t[iVec]->Tpv[2] = fitpv.Tr().tx[iVec];
- t[iVec]->Tpv[3] = fitpv.Tr().ty[iVec];
- t[iVec]->Tpv[4] = fitpv.Tr().qp[iVec];
- t[iVec]->Tpv[5] = fitpv.Tr().z[iVec];
- t[iVec]->Tpv[6] = fitpv.Tr().t[iVec];
-
- t[iVec]->Cpv[0] = fitpv.Tr().C00[iVec];
- t[iVec]->Cpv[1] = fitpv.Tr().C10[iVec];
- t[iVec]->Cpv[2] = fitpv.Tr().C11[iVec];
- t[iVec]->Cpv[3] = fitpv.Tr().C20[iVec];
- t[iVec]->Cpv[4] = fitpv.Tr().C21[iVec];
- t[iVec]->Cpv[5] = fitpv.Tr().C22[iVec];
- t[iVec]->Cpv[6] = fitpv.Tr().C30[iVec];
- t[iVec]->Cpv[7] = fitpv.Tr().C31[iVec];
- t[iVec]->Cpv[8] = fitpv.Tr().C32[iVec];
- t[iVec]->Cpv[9] = fitpv.Tr().C33[iVec];
- t[iVec]->Cpv[10] = fitpv.Tr().C40[iVec];
- t[iVec]->Cpv[11] = fitpv.Tr().C41[iVec];
- t[iVec]->Cpv[12] = fitpv.Tr().C42[iVec];
- t[iVec]->Cpv[13] = fitpv.Tr().C43[iVec];
- t[iVec]->Cpv[14] = fitpv.Tr().C44[iVec];
- t[iVec]->Cpv[15] = fitpv.Tr().C50[iVec];
- t[iVec]->Cpv[16] = fitpv.Tr().C51[iVec];
- t[iVec]->Cpv[17] = fitpv.Tr().C52[iVec];
- t[iVec]->Cpv[18] = fitpv.Tr().C53[iVec];
- t[iVec]->Cpv[19] = fitpv.Tr().C54[iVec];
- t[iVec]->Cpv[20] = fitpv.Tr().C55[iVec];
+ fitpv.Tr().copyToArrays(iVec, t[iVec]->Tpv, t[iVec]->Cpv);
}
if (iter == 1) continue; // only 1.5 iterations
@@ -415,41 +362,15 @@ void L1Algo::L1KFTrackFitter()
fldZ1 = fldZ0;
}
+ L1TrackPar Tl = fit.Tr();
+ if (kGlobal == fTrackingMode) { Tl.qp = fitpv.Tr().qp; }
+
for (int iVec = 0; iVec < nTracks_SIMD; 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.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];
+ Tl.copyToArrays(iVec, t[iVec]->TLast, t[iVec]->CLast);
+ t[iVec]->chi2 = Tl.chi2[iVec];
+ t[iVec]->NDF = (int) Tl.NDF[iVec];
}
+
} // iter
}
}
@@ -460,22 +381,15 @@ void L1Algo::GuessVecNoField(L1Fit& track, fvec& x_last, fvec& x_2last, fvec& y_
fvec dzi = fvec(1.) / (z_end - z_2last);
- tr.x = x_last;
- tr.y = y_last;
- tr.tx = (x_last - x_2last) * dzi;
- tr.ty = (y_last - y_2last) * dzi;
- tr.t = time_last;
- tr.z = z_end;
- tr.chi2 = fvec(0.);
- tr.NDF = fvec(2.);
-
- tr.C20 = tr.C21 = fvec(0.);
- tr.C30 = tr.C31 = tr.C32 = fvec(0.);
- tr.C40 = tr.C41 = tr.C42 = tr.C43 = fvec(0.);
- tr.C50 = tr.C51 = tr.C52 = tr.C53 = tr.C54 = fvec(0.);
- tr.C22 = tr.C33 = vINF;
- tr.C44 = fvec(1.);
- tr.C55 = fvec(dt2_last);
+ tr.x = x_last;
+ tr.y = y_last;
+ tr.tx = (x_last - x_2last) * dzi;
+ tr.ty = (y_last - y_2last) * dzi;
+ tr.t = time_last;
+ tr.z = z_end;
+
+ tr.ResetErrors(1., 1., 1., 1., 1., dt2_last, 1.e2);
+ tr.NDF = fvec(2.);
}
@@ -483,7 +397,7 @@ void L1Algo::GuessVec(L1TrackPar& tr, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fm
// gives nice initial approximation for x,y,tx,ty - almost same as KF fit. qp - is shifted by 4%, resid_ual - ~3.5% (KF fit resid_ual - 1%).
{
- fvec A0, A1 = ZERO, A2 = ZERO, A3 = ZERO, A4 = ZERO, A5 = ZERO, a0, a1 = ZERO, a2 = ZERO, b0, b1 = ZERO, b2 = ZERO;
+ fvec A0, A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0., a0, a1 = 0., a2 = 0., b0, b1 = 0., b2 = 0.;
fvec z0, x, y, z, S, w, wz, wS;
int i = NHits - 1;
@@ -553,7 +467,7 @@ void L1Algo::GuessVec(L1Fit& track, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fmas
{
L1TrackPar& tr = track.Tr();
- fvec A0, A1 = ZERO, A2 = ZERO, A3 = ZERO, A4 = ZERO, A5 = ZERO, a0, a1 = ZERO, a2 = ZERO, b0, b1 = ZERO, b2 = ZERO;
+ fvec A0, A1 = 0., A2 = 0., A3 = 0., A4 = 0., A5 = 0., a0, a1 = 0., a2 = 0., b0, b1 = 0., b2 = 0.;
fvec z0, x, y, z, S, w, wz, wS, time;
time = fvec(0.);
@@ -639,31 +553,11 @@ void L1Algo::FilterFirst(L1Fit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fve
{
L1TrackPar& tr = track.Tr();
- tr.C00 = d_xx;
+ tr.ResetErrors(d_xx, d_yy, 0.1, 0.1, 1., dt2, 1.e2);
tr.C10 = d_xy;
- tr.C11 = d_yy;
- tr.C20 = ZERO;
- tr.C21 = ZERO;
- tr.C22 = vINF;
- tr.C30 = ZERO;
- tr.C31 = ZERO;
- tr.C32 = ZERO;
- tr.C33 = vINF;
- tr.C40 = ZERO;
- tr.C41 = ZERO;
- tr.C42 = ZERO;
- tr.C43 = ZERO;
- tr.C44 = ONE;
- tr.C50 = ZERO;
- tr.C51 = ZERO;
- tr.C52 = ZERO;
- tr.C53 = ZERO;
- tr.C54 = ZERO;
- tr.C55 = dt2;
-
- tr.x = x;
- tr.y = y;
- tr.t = t;
- tr.NDF = -3.0;
- tr.chi2 = ZERO;
+
+ tr.x = x;
+ tr.y = y;
+ tr.t = t;
+ tr.NDF = -3.0;
}
diff --git a/reco/L1/L1Algo/L1TrackPar.cxx b/reco/L1/L1Algo/L1TrackPar.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..5bc71fb8de151dd5e4288bbae10924761bade891
--- /dev/null
+++ b/reco/L1/L1Algo/L1TrackPar.cxx
@@ -0,0 +1,226 @@
+/* Copyright (C) 2007-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Igor Kulakov [committer], Maksym Zyzak */
+
+#include "L1TrackPar.h"
+
+#include <iomanip>
+#include <iostream>
+
+
+void L1TrackPar::Print(int i) const
+{
+ 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 " << x << std::endl;
+ std::cout << " y " << y << std::endl;
+ std::cout << " tx " << tx << std::endl;
+ std::cout << " ty " << ty << std::endl;
+ std::cout << " qp " << qp << std::endl;
+ std::cout << " t " << t << std::endl;
+ std::cout << " v " << fvec(1.) / vi << std::endl;
+ std::cout << " z " << z << 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;
+ std::cout << " c66 " << C66 << std::endl;
+ }
+ else {
+ std::cout << "T = ";
+ std::cout << " x " << x[i];
+ std::cout << " y " << y[i];
+ std::cout << " tx " << tx[i];
+ std::cout << " ty " << ty[i];
+ std::cout << " qp " << qp[i];
+ std::cout << " t " << t[i];
+ std::cout << " v " << 1. / vi[i];
+
+ std::cout << " z " << z[i] << std::endl;
+ std::cout << "C = ";
+ std::cout << " c00 " << C00[i];
+ std::cout << " c11 " << C11[i];
+ std::cout << " c22 " << C22[i];
+ std::cout << " c33 " << C33[i];
+ std::cout << " c44 " << C44[i] << std::endl;
+ std::cout << " c55 " << C55[i] << std::endl;
+ std::cout << " c66 " << C66[i] << std::endl;
+ }
+ PrintCorrelations(i);
+ std::cout.flags(coutFlags);
+}
+
+void L1TrackPar::PrintCorrelations(int i) const
+{
+ auto flagSafe = std::cout.flags();
+ // std::cout.setf(std::ios::scientific, std::ios::floatfield);
+ std::cout << std::setprecision(6);
+
+ if (i == -1) {
+ fvec s0 = sqrt(C00);
+ fvec s1 = sqrt(C11);
+ fvec s2 = sqrt(C22);
+ fvec s3 = sqrt(C33);
+ fvec s4 = sqrt(C44);
+ fvec s5 = sqrt(C55);
+ fvec s6 = sqrt(C66);
+
+ std::cout << "K = " << std::endl;
+ std::cout << " k10 " << C10 / s1 / s0 << std::endl;
+
+ std::cout << "\n k20 " << C20 / s2 / s0 << std::endl;
+ std::cout << " k21 " << C21 / s2 / s1 << std::endl;
+
+ std::cout << "\n k30 " << C30 / s3 / s0 << std::endl;
+ std::cout << " k31 " << C31 / s3 / s1 << std::endl;
+ std::cout << " k32 " << C32 / s3 / s2 << std::endl;
+
+ std::cout << "\n k40 " << C40 / s4 / s0 << std::endl;
+ std::cout << " k41 " << C41 / s4 / s1 << std::endl;
+ std::cout << " k42 " << C42 / s4 / s2 << std::endl;
+ std::cout << " k43 " << C43 / s4 / s3 << std::endl;
+
+ std::cout << "\n k50 " << C50 / s5 / s0 << std::endl;
+ std::cout << " k51 " << C51 / s5 / s1 << std::endl;
+ std::cout << " k52 " << C52 / s5 / s2 << std::endl;
+ std::cout << " k53 " << C53 / s5 / s3 << std::endl;
+ std::cout << " k54 " << C54 / s5 / s4 << std::endl;
+
+ std::cout << "\n k60 " << C60 / s6 / s0 << std::endl;
+ std::cout << " k61 " << C61 / s6 / s1 << std::endl;
+ std::cout << " k62 " << C62 / s6 / s2 << std::endl;
+ std::cout << " k63 " << C63 / s6 / s3 << std::endl;
+ std::cout << " k64 " << C64 / s6 / s4 << std::endl;
+ std::cout << " k65 " << C65 / s6 / s5 << std::endl;
+ }
+ else {
+ float s0 = sqrt(C00[i]);
+ float s1 = sqrt(C11[i]);
+ float s2 = sqrt(C22[i]);
+ float s3 = sqrt(C33[i]);
+ float s4 = sqrt(C44[i]);
+ float s5 = sqrt(C55[i]);
+ float s6 = sqrt(C66[i]);
+
+ std::cout << "K = " << std::endl;
+ std::cout << " " << C10[i] / s1 / s0 << std::endl;
+ std::cout << " " << C20[i] / s2 / s0 << " " << C21[i] / s2 / s1 << std::endl;
+ std::cout << " " << C30[i] / s3 / s0 << " " << C31[i] / s3 / s1 << " " << C32[i] / s3 / s2 << std::endl;
+ std::cout << " " << C40[i] / s4 / s0 << " " << C41[i] / s4 / s1 << " " << C42[i] / s4 / s2 << " "
+ << C43[i] / s4 / s3 << std::endl;
+ std::cout << " " << C50[i] / s5 / s0 << " " << C51[i] / s5 / s1 << " " << C52[i] / s5 / s2 << " "
+ << C53[i] / s5 / s3 << " " << C54[i] / s5 / s4 << std::endl;
+ std::cout << " " << C60[i] / s6 / s0 << " " << C61[i] / s6 / s1 << " " << C62[i] / s6 / s2 << " "
+ << C63[i] / s6 / s3 << " " << C64[i] / s6 / s4 << " " << C65[i] / s6 / s5 << std::endl;
+ }
+ std::cout.flags(flagSafe);
+}
+
+bool L1TrackPar::IsEntryConsistent(bool printWhenWrong, int k) const
+{
+ bool ok = true;
+
+ // verify that all the numbers in the object are valid floats
+ const fvec* memberFirst = &x;
+ const fvec* memberLast = &NDF;
+ for (int i = 0; i < &memberLast - &memberFirst + 1; i++) {
+ if (!std::isfinite(memberFirst[i][k])) {
+ ok = false;
+ if (printWhenWrong) {
+ std::cout << " L1TrackPar member N " << i << ", vector entry " << k
+ << " is not a number: " << memberFirst[i][k];
+ }
+ }
+ }
+
+ // verify diagonal elements.
+ // Cii is a squared dispersion of i-th parameter, it must be positive
+
+ for (int i = 0; i < 7; i++) {
+ if (C(i, i)[k] <= 0.f) {
+ ok = false;
+ if (printWhenWrong) {
+ std::cout << " L1TrackPar: C[" << i << "," << i << "], vec entry " << k << " is not positive: " << C(i, i)[k]
+ << std::endl;
+ }
+ }
+ }
+
+ // verify non-diagonal elements.
+ // Cij/sqrt(Cii*Cjj) is a correlation between i-th and j-th parameter,
+ // it must belong to [-1,1]
+
+ for (int i = 1; i < 7; i++) {
+ for (int j = 0; j < i; j++) {
+ double tolerance = 1.0;
+ if (C(i, j)[k] * C(i, j)[k] > tolerance * (C(i, i)[k] * C(j, j)[k])) {
+ ok = false;
+ if (printWhenWrong) {
+ std::cout << " L1TrackPar: correlation [" << i << "," << j << "], vec entry " << k
+ << " is too large: " << C(i, i)[k] / sqrt(C(i, i)[k] * C(j, j)[k]) << std::endl;
+ }
+ }
+ }
+ }
+
+ // verify triplets of correlations
+ // Kxy * Kxy + Kxz * Kxz + Kyz * Kyz <= 1 + 2 * Kxy * Kxz * Kyz
+
+ for (int i = 2; i < 7; i++) {
+ for (int j = 1; j < i; j++) {
+ for (int m = 0; m < j; m++) {
+ double tolerance = 1.0;
+ double Cxx = C(i, i)[k];
+ double Cyy = C(j, j)[k];
+ double Czz = C(m, m)[k];
+ double Cxy = C(i, j)[k];
+ double Cxz = C(i, m)[k];
+ double Cyz = C(j, m)[k];
+ if (Cxx * Cyz * Cyz + Cyy * Cxz * Cxz + Czz * Cxy * Cxy
+ > tolerance * (Cxx * Cyy * Czz + 2. * Cxy * Cyz * Cxz)) {
+ ok = false;
+ if (printWhenWrong) {
+ double Kxy = Cxy / sqrt(Cxx * Cyy);
+ double Kxz = Cxz / sqrt(Cxx * Czz);
+ double Kyz = Cyz / sqrt(Cyy * Czz);
+ std::cout << " L1TrackPar: correlations between parametetrs " << i << ", " << j << ", " << m
+ << ", vec entry " << k << " are wrong: " << Kxy << " " << Kxz << " " << Kyz << std::endl;
+ std::cout << " inequation: " << Kxy * Kxy + Kxz * Kxz + Kyz * Kyz << " > " << 1 + 2 * Kxy * Kxz * Kyz
+ << std::endl;
+ }
+ }
+ }
+ }
+ }
+
+ if (!ok && printWhenWrong) {
+ std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
+ Print(k);
+ }
+ return ok;
+}
+
+bool L1TrackPar::IsConsistent(bool printWhenWrong, int nFilled) const
+{
+ assert(nFilled <= (int) fvec::size());
+ bool ok = true;
+ if (nFilled < 0) { nFilled = fvec::size(); }
+ for (int i = 0; i < nFilled; ++i) {
+ ok = ok && IsEntryConsistent(printWhenWrong, i);
+ }
+
+ if (!ok && printWhenWrong) {
+ std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
+ if (nFilled == (int) fvec::size()) { std::cout << " All vector elements are filled " << std::endl; }
+ else {
+ std::cout << " Only first " << nFilled << " vector elements are filled " << std::endl;
+ }
+ Print(-1);
+ }
+ return ok;
+}
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 4f60d6d804055586a3191f48449fd15d9e1cb577..25e4fb43a7d60be9a88519c5683282afde5b3986 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -5,19 +5,31 @@
#ifndef L1TrackPar_h
#define L1TrackPar_h 1
-#include <iomanip>
-#include <iostream>
-
+#include "CaSimd.h"
#include "L1Def.h"
+using namespace cbm::algo::ca;
+
class L1TrackPar {
public:
- fvec x {0.}, y {0.}, tx {0.}, ty {0.}, qp {0.}, z {0.}, t {0.};
+ fvec x {0.}, y {0.}, tx {0.}, ty {0.}, qp {0.}, z {0.}, t {0.}, vi {0.};
+
+ fvec C00 {0.},
+
+ C10 {0.}, C11 {0.},
+
+ C20 {0.}, C21 {0.}, C22 {0.},
+
+ C30 {0.}, C31 {0.}, C32 {0.}, C33 {0.},
- fvec 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.};
+ C40 {0.}, C41 {0.}, C42 {0.}, C43 {0.}, C44 {0.},
+
+ C50 {0.}, C51 {0.}, C52 {0.}, C53 {0.}, C54 {0.}, C55 {0.},
+
+ C60 {0.}, C61 {0.}, C62 {0.}, C63 {0.}, C64 {0.}, C65 {0.}, C66 {0.};
+
+ fvec chi2 {0.}, NDF {0.};
L1TrackPar() = default;
@@ -27,43 +39,8 @@ public:
Set(tr, C);
}
- template<typename T>
- void Set(const T* tr, const T* C)
- {
- x = tr[0];
- y = tr[1];
- tx = tr[2];
- ty = tr[3];
- qp = tr[4];
- z = tr[5];
- t = tr[6];
-
- chi2 = 0.;
- NDF = 0.;
-
- 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];
- }
-
+ //template<typename T>
+ //void Set(const T* tr, const T* C);
void SetOneEntry(const int i0, const L1TrackPar& T1, const int i1);
@@ -81,6 +58,9 @@ public:
return c[ind];
}
+
+ void ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec c44, fvec c55, fvec c66);
+
void Print(int i = -1) const;
void PrintCorrelations(int i = -1) const;
@@ -89,134 +69,91 @@ public:
bool IsConsistent(bool printWhenWrong, int nFilled) const;
- void ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec c44, fvec c55)
- {
- C10 = 0.;
- C20 = C21 = 0.;
- C30 = C31 = C32 = 0.;
- C40 = C41 = C42 = C43 = 0.;
- C50 = C51 = C52 = C53 = C54 = 0.;
-
- C00 = c00;
- C11 = c11;
- C22 = c22;
- C33 = c33;
- C44 = c44;
- C55 = c55;
- chi2 = 0.;
- NDF = 0;
- }
+ template<typename T>
+ void copyToArrays(const int iVec, T p[7], T c[28]) const;
+
+ template<typename T>
+ void copyFromArrays(const int iVec, const T p[7], const T c[28]);
+
+ template<typename T>
+ void copyFromArrays(const T p[7], const T c[28]);
} _fvecalignment;
// =============================================================================================
-inline void L1TrackPar::Print(int i) const
+template<typename T>
+inline void L1TrackPar::copyToArrays(const int iVec, T p[7], T c[28]) const
{
- 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 " << x << std::endl;
- std::cout << " y " << y << std::endl;
- std::cout << " tx " << tx << std::endl;
- std::cout << " ty " << ty << std::endl;
- std::cout << " qp " << qp << std::endl;
- std::cout << " t " << t << std::endl;
- std::cout << " z " << z << 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 << " x " << x[i];
- std::cout << " y " << y[i];
- std::cout << " tx " << tx[i];
- std::cout << " ty " << ty[i];
- std::cout << " qp " << qp[i];
- std::cout << " t " << t[i];
- std::cout << " z " << z[i] << std::endl;
- std::cout << "C = ";
- std::cout << " c00 " << C00[i];
- std::cout << " c11 " << C11[i];
- std::cout << " c22 " << C22[i];
- std::cout << " c33 " << C33[i];
- std::cout << " c44 " << C44[i] << std::endl;
- std::cout << " c55 " << C55[i] << std::endl;
+ p[0] = x[iVec];
+ p[1] = y[iVec];
+ p[2] = tx[iVec];
+ p[3] = ty[iVec];
+ p[4] = qp[iVec];
+ p[5] = z[iVec];
+ p[6] = t[iVec];
+ p[7] = vi[iVec];
+
+ const fvec* tC = &(C00);
+ for (int i = 0; i < 28; i++) {
+ c[i] = tC[i][iVec];
}
- PrintCorrelations(i);
- std::cout.flags(coutFlags);
}
-inline void L1TrackPar::PrintCorrelations(int i) const
+template<typename T>
+inline void L1TrackPar::copyFromArrays(const int iVec, const T p[7], const T c[28])
{
- auto flagSafe = std::cout.flags();
- // std::cout.setf(std::ios::scientific, std::ios::floatfield);
- std::cout << std::setprecision(6);
-
- if (i == -1) {
- fvec s0 = sqrt(C00);
- fvec s1 = sqrt(C11);
- fvec s2 = sqrt(C22);
- fvec s3 = sqrt(C33);
- fvec s4 = sqrt(C44);
- fvec s5 = sqrt(C55);
-
- std::cout << "K = " << std::endl;
- std::cout << " k10 " << C10 / s1 / s0 << std::endl;
-
- std::cout << "\n k20 " << C20 / s2 / s0 << std::endl;
- std::cout << " k21 " << C21 / s2 / s1 << std::endl;
-
- std::cout << "\n k30 " << C30 / s3 / s0 << std::endl;
- std::cout << " k31 " << C31 / s3 / s1 << std::endl;
- std::cout << " k32 " << C32 / s3 / s2 << std::endl;
-
- std::cout << "\n k40 " << C40 / s4 / s0 << std::endl;
- std::cout << " k41 " << C41 / s4 / s1 << std::endl;
- std::cout << " k42 " << C42 / s4 / s2 << std::endl;
- std::cout << " k43 " << C43 / s4 / s3 << std::endl;
-
- std::cout << "\n k50 " << C50 / s5 / s0 << std::endl;
- std::cout << " k51 " << C51 / s5 / s1 << std::endl;
- std::cout << " k52 " << C52 / s5 / s2 << std::endl;
- std::cout << " k53 " << C53 / s5 / s3 << std::endl;
- std::cout << " k54 " << C54 / s5 / s4 << std::endl;
- }
- else {
- float s0 = sqrt(C00[i]);
- float s1 = sqrt(C11[i]);
- float s2 = sqrt(C22[i]);
- float s3 = sqrt(C33[i]);
- float s4 = sqrt(C44[i]);
- float s5 = sqrt(C55[i]);
-
- std::cout << "K = " << std::endl;
- std::cout << " " << C10[i] / s1 / s0 << std::endl;
- std::cout << " " << C20[i] / s2 / s0 << " " << C21[i] / s2 / s1 << std::endl;
- std::cout << " " << C30[i] / s3 / s0 << " " << C31[i] / s3 / s1 << " " << C32[i] / s3 / s2 << std::endl;
- std::cout << " " << C40[i] / s4 / s0 << " " << C41[i] / s4 / s1 << " " << C42[i] / s4 / s2 << " "
- << C43[i] / s4 / s3 << std::endl;
- std::cout << " " << C50[i] / s5 / s0 << " " << C51[i] / s5 / s1 << " " << C52[i] / s5 / s2 << " "
- << C53[i] / s5 / s3 << " " << C54[i] / s5 / s4 << std::endl;
- }
- std::cout.flags(flagSafe);
+ x[iVec] = p[0];
+ y[iVec] = p[1];
+ tx[iVec] = p[2];
+ ty[iVec] = p[3];
+ qp[iVec] = p[4];
+ z[iVec] = p[5];
+ t[iVec] = p[6];
+ vi[iVec] = p[7];
+
+ fvec* tC = &(C00);
+ for (int i = 0; i < 28; i++) {
+ tC[i][iVec] = c[i];
+ }
+
+ chi2 = 0.;
+ NDF = 0.;
+}
+
+template<typename T>
+inline void L1TrackPar::copyFromArrays(const T p[7], const T c[28])
+{
+ x = p[0];
+ y = p[1];
+ tx = p[2];
+ ty = p[3];
+ qp = p[4];
+ z = p[5];
+ t = p[6];
+ vi = p[7];
+
+ fvec* tC = &(C00);
+ for (int i = 0; i < 28; i++) {
+ tC[i] = c[i];
+ }
+
+ chi2 = 0.;
+ NDF = 0.;
}
+
inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const int i1)
{
- x[i0] = T1.x[i1];
- y[i0] = T1.y[i1];
- tx[i0] = T1.tx[i1];
- ty[i0] = T1.ty[i1];
- qp[i0] = T1.qp[i1];
- z[i0] = T1.z[i1];
- t[i0] = T1.t[i1];
+ x[i0] = T1.x[i1];
+ y[i0] = T1.y[i1];
+ tx[i0] = T1.tx[i1];
+ ty[i0] = T1.ty[i1];
+ qp[i0] = T1.qp[i1];
+ z[i0] = T1.z[i1];
+ t[i0] = T1.t[i1];
+ vi[i0] = T1.vi[i1];
+
C00[i0] = T1.C00[i1];
C10[i0] = T1.C10[i1];
C11[i0] = T1.C11[i1];
@@ -238,114 +175,36 @@ inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const in
C53[i0] = T1.C53[i1];
C54[i0] = T1.C54[i1];
C55[i0] = T1.C55[i1];
+ C60[i0] = T1.C60[i1];
+ C61[i0] = T1.C61[i1];
+ C62[i0] = T1.C62[i1];
+ C63[i0] = T1.C63[i1];
+ C64[i0] = T1.C64[i1];
+ C65[i0] = T1.C65[i1];
+ C66[i0] = T1.C66[i1];
chi2[i0] = T1.chi2[i1];
NDF[i0] = T1.NDF[i1];
} // SetOneEntry
-inline bool L1TrackPar::IsEntryConsistent(bool printWhenWrong, int k) const
+inline void L1TrackPar::ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec c44, fvec c55, fvec c66)
{
- bool ok = true;
-
- // verify that all the numbers in the object are valid floats
- const fvec* memberFirst = &x;
- const fvec* memberLast = &NDF;
- for (int i = 0; i < &memberLast - &memberFirst + 1; i++) {
- if (!std::isfinite(memberFirst[i][k])) {
- ok = false;
- if (printWhenWrong) {
- std::cout << " L1TrackPar member N " << i << ", vector entry " << k
- << " is not a number: " << memberFirst[i][k];
- }
- }
- }
-
- // verify diagonal elements.
- // Cii is a squared dispersion of i-th parameter, it must be positive
-
- for (int i = 0; i < 6; i++) {
- if (C(i, i)[k] <= 0.f) {
- ok = false;
- if (printWhenWrong) {
- std::cout << " L1TrackPar: C[" << i << "," << i << "], vec entry " << k << " is not positive: " << C(i, i)[k]
- << std::endl;
- }
- }
- }
-
- // verify non-diagonal elements.
- // Cij/sqrt(Cii*Cjj) is a correlation between i-th and j-th parameter,
- // it must belong to [-1,1]
-
- for (int i = 1; i < 6; i++) {
- for (int j = 0; j < i; j++) {
- double tolerance = 1.0;
- if (C(i, j)[k] * C(i, j)[k] > tolerance * (C(i, i)[k] * C(j, j)[k])) {
- ok = false;
- if (printWhenWrong) {
- std::cout << " L1TrackPar: correlation [" << i << "," << j << "], vec entry " << k
- << " is too large: " << C(i, i)[k] / sqrt(C(i, i)[k] * C(j, j)[k]) << std::endl;
- }
- }
- }
- }
-
- // verify triplets of correlations
- // Kxy * Kxy + Kxz * Kxz + Kyz * Kyz <= 1 + 2 * Kxy * Kxz * Kyz
-
- for (int i = 2; i < 6; i++) {
- for (int j = 1; j < i; j++) {
- for (int m = 0; m < j; m++) {
- double tolerance = 1.0;
- double Cxx = C(i, i)[k];
- double Cyy = C(j, j)[k];
- double Czz = C(m, m)[k];
- double Cxy = C(i, j)[k];
- double Cxz = C(i, m)[k];
- double Cyz = C(j, m)[k];
- if (Cxx * Cyz * Cyz + Cyy * Cxz * Cxz + Czz * Cxy * Cxy
- > tolerance * (Cxx * Cyy * Czz + 2. * Cxy * Cyz * Cxz)) {
- ok = false;
- if (printWhenWrong) {
- double Kxy = Cxy / sqrt(Cxx * Cyy);
- double Kxz = Cxz / sqrt(Cxx * Czz);
- double Kyz = Cyz / sqrt(Cyy * Czz);
- std::cout << " L1TrackPar: correlations between parametetrs " << i << ", " << j << ", " << m
- << ", vec entry " << k << " are wrong: " << Kxy << " " << Kxz << " " << Kyz << std::endl;
- std::cout << " inequation: " << Kxy * Kxy + Kxz * Kxz + Kyz * Kyz << " > " << 1 + 2 * Kxy * Kxz * Kyz
- << std::endl;
- }
- }
- }
- }
- }
-
- if (!ok && printWhenWrong) {
- std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
- Print(k);
- }
- return ok;
+ C10 = 0.;
+ C20 = C21 = 0.;
+ C30 = C31 = C32 = 0.;
+ C40 = C41 = C42 = C43 = 0.;
+ C50 = C51 = C52 = C53 = C54 = 0.;
+ C60 = C61 = C62 = C63 = C64 = C65 = 0.;
+
+ C00 = c00;
+ C11 = c11;
+ C22 = c22;
+ C33 = c33;
+ C44 = c44;
+ C55 = c55;
+ C66 = c66;
+ chi2 = 0.;
+ NDF = -6.;
}
-inline bool L1TrackPar::IsConsistent(bool printWhenWrong, int nFilled) const
-{
- assert(nFilled <= (int) fvec::size());
- bool ok = true;
- if (nFilled < 0) { nFilled = fvec::size(); }
- for (int i = 0; i < nFilled; ++i) {
- ok = ok && IsEntryConsistent(printWhenWrong, i);
- }
-
- if (!ok && printWhenWrong) {
- std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
- if (nFilled == (int) fvec::size()) { std::cout << " All vector elements are filled " << std::endl; }
- else {
- std::cout << " Only first " << nFilled << " vector elements are filled " << std::endl;
- }
- Print(-1);
- }
- return ok;
-}
-
-
#endif
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index ab9b8c0aaddb5723ee1d3fa94397a82f903c3d8e..3264abf8aac9a8744d4deefe9a40c7e4790a5a16 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -100,27 +100,13 @@ inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const L1FieldRegion& fld, int
void FilterFirst(L1Fit& fit, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
{
L1TrackPar& tr = fit.Tr();
-
- tr.C00 = dxx;
+ tr.ResetErrors(dxx, dyy, vINF, vINF, 1., 1.e6, 1.e2);
tr.C10 = dxy;
- tr.C11 = dyy;
- tr.C20 = ZERO;
- tr.C21 = ZERO;
- tr.C22 = vINF;
- tr.C30 = ZERO;
- tr.C31 = ZERO;
- tr.C32 = ZERO;
- tr.C33 = vINF;
- tr.C40 = ZERO;
- tr.C41 = ZERO;
- tr.C42 = ZERO;
- tr.C43 = ZERO;
- tr.C44 = ONE;
-
- tr.x = x;
- tr.y = y;
- tr.NDF = -3.0;
- tr.chi2 = ZERO;
+ tr.x = x;
+ tr.y = y;
+ tr.t = 0.;
+ tr.vi = 0.;
+ tr.NDF = -3.0;
}
void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
@@ -181,6 +167,10 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
+ T.t = 0.;
+ T.vi = 0.;
+ T.ResetErrors(1.e2, 1.e2, 1.e4, 1.e4, 1.e4, 1.e6, 1.e2);
+
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
for (i = 0; i < nTracks_SIMD; i++) {
tr[i] = &Tracks[itrack + i]; // current track
@@ -211,6 +201,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
// mass[i] = TDatabasePDG::Instance()->GetParticle(pid)->Mass();
mass[i] = KFParticleDatabase::Instance()->GetMass(pid);
}
+
fit.SetParticleMass(mass);
// get hits of current track