diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 94dac75987641496a6d198437a26da95d9a3c2d7..6a19a64eeb5e0c673aa93dd6537c31f54283cd37 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -641,8 +641,8 @@ InitStatus CbmL1::Init()
trIterDefault.SetMinNhitsStation0(4);
trIterDefault.SetDoubletChi2Cut(7.56327f); // = 1.3449 * 2.f / 3.f; // prob = 0.1
trIterDefault.SetTripletChi2Cut(23.4450f); // = 7.815 * 3; // prob = 0.05
- trIterDefault.SetTripletFinalChi2Cut(12.5);
- trIterDefault.SetTripletLinkChi2(10.);
+ trIterDefault.SetTripletFinalChi2Cut(23.4450f); //new 7.5);
+ trIterDefault.SetTripletLinkChi2(25.); // new:10
trIterDefault.SetTrackChi2Cut(10.);
trIterDefault.SetMaxSlope(2.748f);
trIterDefault.SetMaxSlopePV(1.1f);
@@ -666,6 +666,7 @@ InitStatus CbmL1::Init()
auto trackingIterFastPrim2 = L1CAIteration(trIterDefault, "FastPrim2Iter");
trackingIterFastPrim2.SetTripletChi2Cut(21.1075f);
+ trackingIterFastPrim2.SetTripletFinalChi2Cut(21.1075f);
trackingIterFastPrim2.SetPickGather(3.0f);
trackingIterFastPrim2.SetMaxInvMom(1.0 / 0.5);
trackingIterFastPrim2.SetMaxDZ(0);
@@ -674,12 +675,15 @@ InitStatus CbmL1::Init()
auto trackingIterAllSec = L1CAIteration(trIterDefault, "AllSecIter");
trackingIterAllSec.SetTripletChi2Cut(18.7560f); // = 6.252 * 3; // prob = 0.1
+ trackingIterAllSec.SetTripletFinalChi2Cut(18.7560f);
trackingIterAllSec.SetMaxInvMom(1.0 / 0.1);
trackingIterAllSec.SetTargetPosSigmaXY(10, 10);
+ trackingIterAllSec.SetMaxSlopePV(1.5);
trackingIterAllSec.SetPrimaryFlag(false);
auto trackingIterFastPrimJump = L1CAIteration(trIterDefault, "FastPrimJumpIter");
trackingIterFastPrimJump.SetTripletChi2Cut(21.1075f); // prob = 0.01
+ trackingIterFastPrimJump.SetTripletFinalChi2Cut(21.1075f);
trackingIterFastPrimJump.SetPickGather(3.0f);
trackingIterFastPrimJump.SetMaxInvMom(1.0 / 0.5);
trackingIterFastPrimJump.SetMaxDZ(0);
@@ -689,6 +693,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllPrimJump = L1CAIteration(trIterDefault, "AllPrimJumpIter");
trackingIterAllPrimJump.SetTripletChi2Cut(18.7560f);
+ trackingIterAllPrimJump.SetTripletFinalChi2Cut(18.7560f);
trackingIterAllPrimJump.SetMaxInvMom(1.0 / 0.1);
trackingIterAllPrimJump.SetTargetPosSigmaXY(5, 5);
trackingIterAllPrimJump.SetPrimaryFlag(true);
@@ -696,6 +701,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllSecJump = L1CAIteration(trIterDefault, "AllSecJumpIter");
trackingIterAllSecJump.SetTripletChi2Cut(21.1075f);
+ trackingIterAllSecJump.SetTripletFinalChi2Cut(21.1075f);
trackingIterAllSecJump.SetMaxInvMom(1.0 / 0.1);
trackingIterAllSecJump.SetMaxSlopePV(1.5f);
trackingIterAllSecJump.SetTargetPosSigmaXY(10, 10);
@@ -711,6 +717,7 @@ InitStatus CbmL1::Init()
auto trackingIterAllSecE = L1CAIteration(trIterDefault, "AllSecEIter");
trackingIterAllSecE.SetTripletChi2Cut(18.7560f);
+ trackingIterAllSecE.SetTripletFinalChi2Cut(18.7560f);
trackingIterAllSecE.SetMaxInvMom(1.0 / 0.05);
trackingIterAllSecE.SetMaxSlopePV(1.5f);
trackingIterAllSecE.SetTargetPosSigmaXY(10, 10);
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 8487f85a7509da9b98fd8f058d7d23439d6d0098..5b855ce4a44726997cf4d025b9f42b7818a12b2c 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -455,7 +455,7 @@ void CbmL1::EfficienciesPerformance()
L1_NTRA.PrintEff(/*ifPrintTableToLog = */ true, false);
cout << "Reconstructible MC tracks/event: "
<< (double(L1_NTRA.mc.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS)) << endl;
- cout << "Reconsted MC tracks/event: "
+ cout << "Reconstructed MC tracks/event: "
<< (double(L1_NTRA.reco.counters[L1_NTRA.indices["total"]]) / double(L1_NEVENTS)) << endl;
cout << endl;
cout << "CA Track Finder: " << L1_CATIME / L1_NEVENTS << (fLegacyEventMode ? " s/ev" : " s/time slice") << endl
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index e65c4cbb152b704261fd1abb41b87be0f60cd395..b565926ac306cd9723cd9e921b6f14bcefd9b97b 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -969,11 +969,14 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
/// Refit Triplets
- ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2");
+ ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2:ndf");
- L1TrackParFit fitNew;
- fitNew.SetParticleMass(GetDefaultParticleMass());
+ L1TrackParFit fit;
+ if (fpCurrentIteration->GetElectronFlag()) { fit.SetParticleMass(L1Constants::phys::kElectronMass); }
+ else {
+ fit.SetParticleMass(GetDefaultParticleMass());
+ }
const int NHits = 3; // triplets
// prepare data
@@ -989,7 +992,10 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
const L1Station& star = fParameters.GetStation(istar);
bool isMomentumFitted = ((stal.fieldStatus != 0) || (stam.fieldStatus != 0) || (star.fieldStatus != 0));
- fvec ndf = isMomentumFitted ? -5 : -4;
+ bool isTimeFitted = ((stal.timeInfo != 0) || (stam.timeInfo != 0) || (star.timeInfo != 0));
+ fvec ndf = -4; // straight line
+ ndf += isMomentumFitted ? -1 : 0;
+ ndf += isTimeFitted ? -1 : 0;
for (int i3 = 0; i3 < n3; ++i3) {
int i3_V = i3 / fvec::size();
@@ -1012,7 +1018,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 du2[NHits], dv2[NHits], dt2[NHits];
+ fscal du[NHits], dv[NHits], dt[NHits], du2[NHits], dv2[NHits], dt2[NHits];
for (int ih = 0; ih < NHits; ++ih) {
const L1Hit& hit = fInputData.GetHit(ihit[ih]);
@@ -1021,15 +1027,20 @@ 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;
- du2[ih] = hit.du * hit.du;
- dv2[ih] = hit.dv * hit.dv;
- dt2[ih] = hit.dt * hit.dt;
+ 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;
};
// find the field along the track
L1FieldValue B[3] _fvecalignment;
L1FieldRegion fld _fvecalignment;
+ L1FieldRegion fldTarget _fvecalignment;
fvec tx[3] = {(x[1] - x[0]) / (z[1] - z[0]), (x[2] - x[0]) / (z[2] - z[0]), (x[2] - x[1]) / (z[2] - z[1])};
fvec ty[3] = {(y[1] - y[0]) / (z[1] - z[0]), (y[2] - y[0]) / (z[2] - z[0]), (y[2] - y[1]) / (z[2] - z[1])};
@@ -1039,9 +1050,9 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
};
fld.Set(B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z);
+ fldTarget.Set(fTargB, fTargZ, B[0], sta[0].z, B[1], sta[1].z);
-
- L1TrackPar& T = fitNew.fTr;
+ L1TrackPar& T = fit.fTr;
// initial parameters
{
@@ -1060,6 +1071,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
if (qp0[0] < -GetMaxInvMom()) { qp0 = -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.;
int ih0 = 0;
@@ -1067,22 +1079,40 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.y = y[ih0];
T.z = z[ih0];
T.t = t[ih0];
- T.C55 = dt2[ih0];
- fitNew.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One());
- fitNew.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
+ //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du[ih0], dv[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);
+
+ { // add the target constraint
+ fvec eX, eY, J04, J14;
+ fvec dz;
+ dz = fTargZ - T.z;
+ L1ExtrapolateJXY0(T.tx, T.ty, dz, fldTarget, eX, eY, J04, J14);
+ fvec J[6];
+ J[0] = dz;
+ J[1] = 0;
+ J[2] = J04;
+ J[3] = 0;
+ J[4] = dz;
+ J[5] = J14;
+ L1FilterVtx(T, fTargX, fTargY, TargetXYInfo, eX, eY, J);
+ }
for (int ih = 1; ih < NHits; ++ih) {
- fitNew.Extrapolate(z[ih], qp0, fld, fvec::One());
+ fit.Extrapolate(z[ih], qp0, fld, fvec::One());
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fitNew.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
}
else {
- fitNew.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
+ fit.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
}
- if (ista[ih] == fNstationsBeforePipe) { fitNew.AddPipeMaterial(qp0, fvec::One()); }
- fitNew.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
- fitNew.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
+ //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);
}
}
@@ -1104,26 +1134,29 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
T.t = t[ih0];
T.C55 = dt2[ih0];
- fitNew.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One());
- fitNew.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
+ //std::tie(T.C00, T.C10, T.C11) = sta[ih0].FormXYCovarianceMatrix(du[ih0], dv[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);
for (int ih = NHits - 2; ih >= 0; --ih) {
- fitNew.Extrapolate(z[ih], qp0, fld, fvec::One());
+ fit.Extrapolate(z[ih], qp0, fld, fvec::One());
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fitNew.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
}
else {
- fitNew.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
+ fit.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
}
- if (ista[ih] == fNstationsBeforePipe - 1) { fitNew.AddPipeMaterial(qp0, fvec::One()); }
- fitNew.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
- fitNew.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
+ //if (ista[ih] == fNstationsBeforePipe - 1) { 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);
}
}
} // for iiter
//cout << " chi2 " << T3.chi2[i3_4] << " " << T.chi2[0] << endl;
- //T.chi2 = (fscal) T3.chi2[i3_4];
+ T.chi2 = (fscal) T3.chi2[i3_4]; //SG!!!
T3.SetOneEntry(i3_4, T, i3_4);
{
@@ -1134,7 +1167,9 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
float ev = 0;
float chi2 = T.chi2[i3_4];
- ca::tools::Debugger::Instance().FillNtuple("triplets", ev, mc1, istal, mctr.p, mctr.x, mctr.y, mctr.z, chi2);
+ float nd = T.NDF[i3_4];
+ ca::tools::Debugger::Instance().FillNtuple("triplets", ev, mc1, istal, mctr.p, mctr.x, mctr.y, mctr.z, chi2,
+ nd);
}
}
} //i3
@@ -1181,7 +1216,7 @@ inline void L1Algo::findTripletsStep3( // input
//if (!T3.IsEntryConsistent(false, i3_4)) continue;
// select
- fscal chi2 = T3.chi2[i3_4];
+ fscal chi2 = T3.chi2[i3_4]; // / T3.NDF[i3_4];
const L1HitIndex_t ihitl = hitsl_3[i3] + HitsUnusedStartIndex[istal];
const L1HitIndex_t ihitm = hitsm_3[i3] + HitsUnusedStartIndex[istam];
@@ -1543,13 +1578,14 @@ inline void L1Algo::CreatePortionOfTriplets(
L1_assert(hitsr_3[i] < HitsUnusedStopIndex[istar] - HitsUnusedStartIndex[istar]);
/// Add the right hits to parameters estimation.
- /*
+
+
Tindex n3_V = (n3 + fvec::size() - 1) / fvec::size();
findTripletsStep1( // input
n3_V, star, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER,
// output
T_3);
-*/
+
/// refit
findTripletsStep2(n3, istal, istam, istar, T_3, hitsl_3, hitsm_3, hitsr_3, 1);