diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 6e430593537c0e1da6540e4d45312d5e9795f03f..4aeed635f919d5788c94d0d4c2133e16270ae9c5 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -223,6 +223,7 @@ InitStatus CbmL1::Init()
fUseMUCH = 0;
fUseTRD = 1;
fUseTOF = 0;
+ L1FieldRegion::gkUseOriginalField = true;
fInitManager.DevSetIgnoreHitSearchAreas(true);
//fInitManager.DevSetFitSingletsFromTarget(true);
//fInitManager.DevSetIsMatchDoubletsViaMc(true);
diff --git a/reco/L1/L1Algo/L1Field.cxx b/reco/L1/L1Algo/L1Field.cxx
index 6119b5b159db5c6b75b09052699a19dbbe29276d..9f4091592bba77558919696749d5830619e11102 100644
--- a/reco/L1/L1Algo/L1Field.cxx
+++ b/reco/L1/L1Algo/L1Field.cxx
@@ -4,6 +4,8 @@
#include "L1Field.h"
+#include "CbmKF.h"
+
#include <iomanip>
#include <iostream>
#include <sstream>
@@ -14,6 +16,8 @@
// L1FieldValue methods
//
+bool L1FieldRegion::gkUseOriginalField = false;
+
//----------------------------------------------------------------------------------------------------------------------
//
void L1FieldValue::CheckConsistency() const
@@ -182,15 +186,46 @@ L1FieldValue L1FieldRegion::Get(const fvec z)
return B;
}
+//----------------------------------------------------------------------------------------------------------------------
+// TODO: Should it be inline? (S.Zharko)
+void L1FieldRegion::Get(const fvec x, const fvec y, const fvec z, fvec* B) const
+{
+ if (gkUseOriginalField) {
+ for (size_t i = 0; i < fvec::size(); i++) {
+ double inPos[3] = {x[i], y[i], z[i]};
+ double outB[3];
+ CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB);
+ B[0][i] = outB[0];
+ B[1][i] = outB[1];
+ B[2][i] = outB[2];
+ }
+ }
+ else {
+ Get(z, B);
+ }
+}
+
//----------------------------------------------------------------------------------------------------------------------
// TODO: Should it be inline? (S.Zharko)
void L1FieldRegion::Get(const fvec z_, fvec* B) const
{
- fvec dz = z_ - z0;
- fvec dz2 = dz * dz;
- B[0] = cx0 + cx1 * dz + cx2 * dz2;
- B[1] = cy0 + cy1 * dz + cy2 * dz2;
- B[2] = cz0 + cz1 * dz + cz2 * dz2;
+ if (gkUseOriginalField) {
+ for (size_t i = 0; i < fvec::size(); i++) {
+ double inPos[3] = {0., 0., z_[i]};
+ double outB[3];
+ CbmKF::Instance()->GetMagneticField()->GetFieldValue(inPos, outB);
+ B[0][i] = outB[0];
+ B[1][i] = outB[1];
+ B[2][i] = outB[2];
+ }
+ }
+ else {
+ fvec dz = z_ - z0;
+ fvec dz2 = dz * dz;
+ B[0] = cx0 + cx1 * dz + cx2 * dz2;
+ B[1] = cy0 + cy1 * dz + cy2 * dz2;
+ B[2] = cz0 + cz1 * dz + cz2 * dz2;
+ }
}
//----------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1Field.h b/reco/L1/L1Algo/L1Field.h
index 2b39ae1bb32a7ed879260fcd99ea9e445ccb7142..4d51fbc7823ad05953baf46d6c0a93ab31f067ee 100644
--- a/reco/L1/L1Algo/L1Field.h
+++ b/reco/L1/L1Algo/L1Field.h
@@ -108,6 +108,13 @@ public:
// TODO: Probably we need a const specifier here, because the method does not change the fields
L1FieldValue Get(const fvec z);
+ /// Gets the field vector and writes it into B pointer
+ /// \param x x-coordinate of the point to calculate the field
+ /// \param y y-coordinate of the point to calculate the field
+ /// \param z z-coordinate of the point to calculate the field
+ /// \param B pointer to the output fvec array of the magnetic field
+ void Get(const fvec x, const fvec y, const fvec z, fvec* B) const;
+
/// Gets the field vector and writes it into B pointer
/// \param z_ z-coordinate of the point to calculate the field
/// \param B pointer to the output fvec array of the magnetic field
@@ -164,6 +171,8 @@ public:
std::string ToString(int indentLevel = 0) const;
public:
+ static bool gkUseOriginalField;
+
// TODO: Probably it's better to have arrays instead of separate fvec values? (S.Zharko)
// Bx(z) = cx0 + cx1*(z-z0) + cx2*(z-z0)^2
fvec cx0 {0.f};
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 6fa007c17a350a44d1ed1a21ef17146ea5154754..5cf80ecedb6a21cdc8453d0641756c3704f605f4 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -113,7 +113,6 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
fvec fldZ1 = sta2.z;
fvec fldZ2 = sta0.z;
-
sta1.fieldSlice.GetFieldValue(x1, y1, fldB0);
sta2.fieldSlice.GetFieldValue(x2, y2, fldB1);
sta0.fieldSlice.GetFieldValue(x0, y0, fldB2);
@@ -331,9 +330,8 @@ void L1Algo::L1KFTrackFitter()
L1FieldValue fldB01, fldB11, fldB21 _fvecalignment;
L1FieldRegion fld1 _fvecalignment;
- const int nHits = fParameters.GetNstationsActive();
- int iVec = 0, i = 0;
- int nTracks_SIMD = fvec::size();
+ const int nStations = fParameters.GetNstationsActive();
+ int nTracks_SIMD = fvec::size();
L1TrackPar T; // fitting parametr coresponding to current track
L1TrackParFit T1; // fitting parametr coresponding to current track
@@ -399,34 +397,44 @@ void L1Algo::L1KFTrackFitter()
fvec ZSta[L1Constants::size::kMaxNstations];
- for (int iHit = 0; iHit < nHits; iHit++) {
- ZSta[iHit] = sta[iHit].z;
+ for (int ista = 0; ista < nStations; ista++) {
+ ZSta[ista] = sta[ista].z;
}
unsigned short N_vTracks = fTracks.size(); // number of tracks processed per one SSE register
for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
+
if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
- for (i = 0; i < nTracks_SIMD; i++)
- t[i] = &fTracks[itrack + i]; // current track
+ for (int i = 0; i < nTracks_SIMD; i++) {
+ t[i] = &fTracks[itrack + i];
+ }
+
+ // fill the rest of the SIMD vectors with something reasonable
+ for (uint i = nTracks_SIMD; i < fvec::size(); i++) {
+ t[i] = &fTracks[itrack];
+ }
// get hits of current track
- for (i = 0; i < nHits; i++) {
- w[i] = ZERO;
- w_time[i] = ZERO;
- z[i] = ZSta[i];
+ for (int ista = 0; ista < nStations; ista++) {
+ w[ista] = ZERO;
+ w_time[ista] = ZERO;
+ z[ista] = ZSta[ista];
}
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
+ for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
+
int nHitsTrack = t[iVec]->NHits;
int iSta[L1Constants::size::kMaxNstations];
- for (i = 0; i < nHitsTrack; i++) {
+
+ for (int ih = 0; ih < nHitsTrack; ih++) {
+
const L1Hit& hit = fInputData.GetHit(fRecoHits[start_hit++]);
const int ista = hit.iSt;
- iSta[i] = ista;
+ iSta[ih] = ista;
w[ista][iVec] = 1.;
- if (ista > fNstationsBeforePipe) w_time[ista][iVec] = 1.;
+ if (sta[ista].timeInfo) { w_time[ista][iVec] = 1.; }
u[ista][iVec] = hit.u;
v[ista][iVec] = hit.v;
@@ -444,7 +452,7 @@ void L1Algo::L1KFTrackFitter()
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
- if (i == 0) {
+ if (ih == 0) {
z_start[iVec] = z[ista][iVec];
x_first[iVec] = x[ista][iVec];
y_first[iVec] = y[ista][iVec];
@@ -457,7 +465,7 @@ void L1Algo::L1KFTrackFitter()
staFirst.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
staFirst.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec];
}
- else if (i == nHitsTrack - 1) {
+ else if (ih == nHitsTrack - 1) {
z_end[iVec] = z[ista][iVec];
x_last[iVec] = x[ista][iVec];
y_last[iVec] = y[ista][iVec];
@@ -474,8 +482,8 @@ void L1Algo::L1KFTrackFitter()
fscal prevZ = z_end[iVec];
fscal cursy = 0., curSy = 0.;
- for (i = nHitsTrack - 1; i >= 0; i--) {
- const int ista = iSta[i];
+ for (int ih = nHitsTrack - 1; ih >= 0; ih--) {
+ const int ista = iSta[ih];
const L1Station& st = fParameters.GetStation(ista);
const fscal curZ = z[ista][iVec];
fscal dZ = curZ - prevZ;
@@ -488,95 +496,99 @@ void L1Algo::L1KFTrackFitter()
}
- //fit backward
+ GuessVec(T, x, y, z, Sy, w, nStations, &z_end);
+ GuessVec(T1, x, y, z, Sy, w, nStations, &z_end, time, w_time);
- GuessVec(T, x, y, z, Sy, w, nHits, &z_end);
- GuessVec(T1, x, y, z, Sy, w, nHits, &z_end, time, w_time);
+ if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) {
+ T1.fqp = fvec(0.);
+ T.qp = fvec(0.);
+ }
for (int iter = 0; iter < 2; iter++) { // 1.5 iterations
- fvec qp0 = T.qp;
-
+ fvec qp0 = T.qp;
fvec qp01 = T1.fqp;
- i = nHits - 1;
+ // fit backward
+
+ 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.));
FilterFirst(T, x_last, y_last, staLast);
// FilterFirst( T1, x_last, y_last, time_last, time_er_last, staLast );
FilterFirst(T1, x_last, y_last, time_last, time_er_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
- T1.FilterTime(time[i], timeEr[i], w_time[i], sta[i].timeInfo);
+ T1.FilterTime(time[ista], timeEr[ista], w_time[ista], sta[ista].timeInfo);
// fit.L1AddMaterial( T, sta[i].materialInfo, qp0, 1 );
- fldZ1 = z[i];
-
- sta[i].fieldSlice.GetFieldValue(T.x, T.y, fldB1);
+ fldZ1 = z[ista];
+ sta[ista].fieldSlice.GetFieldValue(T.x, T.y, fldB1);
- fldB1.Combine(fB[i], w[i]);
+ fldB1.Combine(fB[ista], w[ista]);
- fldZ2 = z[i - 2];
+ fldZ2 = z[ista - 2];
fvec dz = fldZ2 - fldZ1;
- sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fldB2);
- fldB2.Combine(fB[i - 2], w[i - 2]);
+ sta[ista].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fldB2);
+ fldB2.Combine(fB[ista - 2], w[ista - 2]);
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- for (--i; i >= 0; i--) {
- fldZ0 = z[i];
+ for (--ista; ista >= 0; ista--) {
+
+ fldZ0 = z[ista];
dz = (fldZ1 - fldZ0);
- sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fldB0);
- fldB0.Combine(fB[i], w[i]);
+ sta[ista].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fldB0);
+ fldB0.Combine(fB[ista], w[ista]);
fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
- fmask initialised = (z[i] < z_end) & (z_start <= z[i]);
- fvec w1 = iif(initialised, w[i], fvec::Zero());
- fvec w1_time = iif(initialised, w_time[i], fvec::Zero());
- fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
+ fmask initialised = (z[ista] < z_end) & (z_start <= z[ista]);
+ fvec w1 = iif(initialised, w[ista], fvec::Zero());
+ fvec wExtr = iif(initialised, fvec::One(), fvec::Zero());
+ fvec w1_time = iif(initialised, w_time[ista], fvec::Zero());
fld1 = fld;
+ L1Extrapolate(T, z[ista], qp0, fld, &wExtr);
- T1.Extrapolate(z[i], qp01, fld1, &w1);
+ T1.Extrapolate(z[ista], qp01, fld1, wExtr);
// T1.ExtrapolateLine(z[i]);
- L1Extrapolate(T, z[i], qp0, fld, &w1);
-
- if (i == fNstationsBeforePipe - 1) {
-
- fit.L1AddPipeMaterial(T, qp0, wIn);
- fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wIn);
-
+ if (ista == fNstationsBeforePipe - 1) {
+ fit.L1AddPipeMaterial(T, qp0, wExtr);
+ fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(1.f), wExtr);
- T1.L1AddPipeMaterial(qp01, wIn);
- T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(1.f), wIn);
+ T1.L1AddPipeMaterial(qp01, wExtr);
+ T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(1.f), wExtr);
}
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, fvec(1.f), wIn);
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, wExtr);
+ fit.EnergyLossCorrection(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec(1.f), wExtr);
- T1.L1AddMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, wIn);
- T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, fvec(1.f), wIn);
+ T1.L1AddMaterial(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, wExtr);
+ T1.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, fvec(1.f), wExtr);
}
else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
- T1.L1AddMaterial(sta[i].materialInfo, qp01, wIn);
+ fit.L1AddMaterial(T, sta[ista].materialInfo, qp0, wExtr);
+ T1.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr);
}
- L1UMeasurementInfo info = sta[i].frontInfo;
- info.sigma2 = d_u[i] * d_u[i];
- L1Filter(T, sta[i].frontInfo, u[i], w1);
- T1.Filter(info, u[i], w1);
+ L1UMeasurementInfo info = sta[ista].frontInfo;
+ info.sigma2 = d_u[ista] * d_u[ista];
- info = sta[i].backInfo;
- info.sigma2 = d_v[i] * d_v[i];
+ L1Filter(T, info, u[ista], w1);
+ T1.Filter(info, u[ista], w1);
- L1Filter(T, sta[i].backInfo, v[i], w1);
- T1.Filter(info, v[i], w1);
+ info = sta[ista].backInfo;
+ info.sigma2 = d_v[ista] * d_v[ista];
- T1.FilterTime(time[i], timeEr[i], w1_time, sta[i].timeInfo);
+ L1Filter(T, info, v[ista], w1);
+ T1.Filter(info, v[ista], w1);
+ T1.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo);
fldB2 = fldB1;
fldZ2 = fldZ1;
@@ -585,127 +597,156 @@ void L1Algo::L1KFTrackFitter()
}
// fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec]->TFirst[0] = T1.fx[iVec];
- t[iVec]->TFirst[1] = T1.fy[iVec];
- t[iVec]->TFirst[2] = T1.ftx[iVec];
- t[iVec]->TFirst[3] = T1.fty[iVec];
- t[iVec]->TFirst[4] = T1.fqp[iVec];
- t[iVec]->TFirst[5] = T1.fz[iVec];
- t[iVec]->TFirst[6] = T1.ft[iVec];
-
- t[iVec]->CFirst[0] = T1.C00[iVec];
- t[iVec]->CFirst[1] = T1.C10[iVec];
- t[iVec]->CFirst[2] = T1.C11[iVec];
- t[iVec]->CFirst[3] = T1.C20[iVec];
- t[iVec]->CFirst[4] = T1.C21[iVec];
- t[iVec]->CFirst[5] = T1.C22[iVec];
- t[iVec]->CFirst[6] = T1.C30[iVec];
- t[iVec]->CFirst[7] = T1.C31[iVec];
- t[iVec]->CFirst[8] = T1.C32[iVec];
- t[iVec]->CFirst[9] = T1.C33[iVec];
- t[iVec]->CFirst[10] = T1.C40[iVec];
- t[iVec]->CFirst[11] = T1.C41[iVec];
- t[iVec]->CFirst[12] = T1.C42[iVec];
- t[iVec]->CFirst[13] = T1.C43[iVec];
- t[iVec]->CFirst[14] = T1.C44[iVec];
- t[iVec]->CFirst[15] = T1.C50[iVec];
- t[iVec]->CFirst[16] = T1.C51[iVec];
- t[iVec]->CFirst[17] = T1.C52[iVec];
- t[iVec]->CFirst[18] = T1.C53[iVec];
- t[iVec]->CFirst[19] = T1.C54[iVec];
- t[iVec]->CFirst[20] = T1.C55[iVec];
-
- t[iVec]->chi2 = T1.chi2[iVec];
- t[iVec]->NDF = (int) T1.NDF[iVec];
- }
-
// extrapolate to the PV region
- L1TrackParFit T1PV = T1;
- T1PV.Extrapolate(0.f, T1PV.fqp, fld);
- // T1PV.ExtrapolateLine(0.f);
+ L1TrackParFit T1pv = T1;
+ {
+ L1UMeasurementInfo vtxInfoX;
+ vtxInfoX.cos_phi = 1.;
+ vtxInfoX.sin_phi = 0.;
+ vtxInfoX.sigma2 = fvec(1.e-4);
+
+ L1UMeasurementInfo vtxInfoY;
+ vtxInfoY.cos_phi = 0;
+ vtxInfoY.sin_phi = 1.;
+ vtxInfoY.sigma2 = fvec(1.e-4);
+
+ if (kGlobal == fTrackingMode) {
+ fvec qp00 = T1.fqp;
+ for (int vtxIter = 0; vtxIter < 2; vtxIter++) {
+ T1pv = T1;
+ T1pv.fqp = qp00;
+ T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fqp, fld, fvec(1.));
+ T1pv.Filter(vtxInfoX, fParameters.GetTargetPositionX());
+ T1pv.Filter(vtxInfoY, fParameters.GetTargetPositionY());
+ qp00 = T1pv.fqp;
+ }
+ }
+ else {
+ T1pv = T1;
+ T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fqp, fld, fvec(1.));
+ }
+ }
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec]->Tpv[0] = T1PV.fx[iVec];
- t[iVec]->Tpv[1] = T1PV.fy[iVec];
- t[iVec]->Tpv[2] = T1PV.ftx[iVec];
- t[iVec]->Tpv[3] = T1PV.fty[iVec];
- t[iVec]->Tpv[4] = T1PV.fqp[iVec];
- t[iVec]->Tpv[5] = T1PV.fz[iVec];
- t[iVec]->Tpv[6] = T1PV.ft[iVec];
-
- t[iVec]->Cpv[0] = T1PV.C00[iVec];
- t[iVec]->Cpv[1] = T1PV.C10[iVec];
- t[iVec]->Cpv[2] = T1PV.C11[iVec];
- t[iVec]->Cpv[3] = T1PV.C20[iVec];
- t[iVec]->Cpv[4] = T1PV.C21[iVec];
- t[iVec]->Cpv[5] = T1PV.C22[iVec];
- t[iVec]->Cpv[6] = T1PV.C30[iVec];
- t[iVec]->Cpv[7] = T1PV.C31[iVec];
- t[iVec]->Cpv[8] = T1PV.C32[iVec];
- t[iVec]->Cpv[9] = T1PV.C33[iVec];
- t[iVec]->Cpv[10] = T1PV.C40[iVec];
- t[iVec]->Cpv[11] = T1PV.C41[iVec];
- t[iVec]->Cpv[12] = T1PV.C42[iVec];
- t[iVec]->Cpv[13] = T1PV.C43[iVec];
- t[iVec]->Cpv[14] = T1PV.C44[iVec];
- t[iVec]->Cpv[15] = T1PV.C50[iVec];
- t[iVec]->Cpv[16] = T1PV.C51[iVec];
- t[iVec]->Cpv[17] = T1PV.C52[iVec];
- t[iVec]->Cpv[18] = T1PV.C53[iVec];
- t[iVec]->Cpv[19] = T1PV.C54[iVec];
- t[iVec]->Cpv[20] = T1PV.C55[iVec];
+ L1TrackParFit Tf = T1;
+ if (kGlobal == fTrackingMode) { Tf = T1pv; }
+
+ for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
+ t[iVec]->TFirst[0] = Tf.fx[iVec];
+ t[iVec]->TFirst[1] = Tf.fy[iVec];
+ t[iVec]->TFirst[2] = Tf.ftx[iVec];
+ t[iVec]->TFirst[3] = Tf.fty[iVec];
+ t[iVec]->TFirst[4] = Tf.fqp[iVec];
+ if (kGlobal == fTrackingMode) { t[iVec]->TFirst[4] = T1pv.fqp[iVec]; }
+ t[iVec]->TFirst[5] = Tf.fz[iVec];
+ t[iVec]->TFirst[6] = Tf.ft[iVec];
+
+ t[iVec]->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];
+
+ 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] = T1pv.fx[iVec];
+ t[iVec]->Tpv[1] = T1pv.fy[iVec];
+ t[iVec]->Tpv[2] = T1pv.ftx[iVec];
+ t[iVec]->Tpv[3] = T1pv.fty[iVec];
+ t[iVec]->Tpv[4] = T1pv.fqp[iVec];
+ t[iVec]->Tpv[5] = T1pv.fz[iVec];
+ t[iVec]->Tpv[6] = T1pv.ft[iVec];
+
+ t[iVec]->Cpv[0] = T1pv.C00[iVec];
+ t[iVec]->Cpv[1] = T1pv.C10[iVec];
+ t[iVec]->Cpv[2] = T1pv.C11[iVec];
+ t[iVec]->Cpv[3] = T1pv.C20[iVec];
+ t[iVec]->Cpv[4] = T1pv.C21[iVec];
+ t[iVec]->Cpv[5] = T1pv.C22[iVec];
+ t[iVec]->Cpv[6] = T1pv.C30[iVec];
+ t[iVec]->Cpv[7] = T1pv.C31[iVec];
+ t[iVec]->Cpv[8] = T1pv.C32[iVec];
+ t[iVec]->Cpv[9] = T1pv.C33[iVec];
+ t[iVec]->Cpv[10] = T1pv.C40[iVec];
+ t[iVec]->Cpv[11] = T1pv.C41[iVec];
+ t[iVec]->Cpv[12] = T1pv.C42[iVec];
+ t[iVec]->Cpv[13] = T1pv.C43[iVec];
+ t[iVec]->Cpv[14] = T1pv.C44[iVec];
+ t[iVec]->Cpv[15] = T1pv.C50[iVec];
+ t[iVec]->Cpv[16] = T1pv.C51[iVec];
+ t[iVec]->Cpv[17] = T1pv.C52[iVec];
+ t[iVec]->Cpv[18] = T1pv.C53[iVec];
+ t[iVec]->Cpv[19] = T1pv.C54[iVec];
+ t[iVec]->Cpv[20] = T1pv.C55[iVec];
+ }
if (iter == 1) continue; // only 1.5 iterations
+
// fit forward
- i = 0;
+ ista = 0;
FilterFirst(T, x_first, y_first, staFirst);
// FilterFirst( T1, x_first, y_first, time_first, time_er_first, staFirst);
FilterFirst(T1, x_first, y_first, time_first, time_er_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
- T1.FilterTime(time[i], timeEr[i], w_time[i], sta[i].timeInfo);
+ T1.FilterTime(time[ista], timeEr[ista], w_time[ista], sta[ista].timeInfo);
// fit.L1AddMaterial( T, sta[i].materialInfo, qp0, 1 );
qp0 = T.qp;
qp01 = T1.fqp;
- fldZ1 = z[i];
- sta[i].fieldSlice.GetFieldValue(T.x, T.y, fldB1);
- fldB1.Combine(fB[i], w[i]);
+ fldZ1 = z[ista];
+ sta[ista].fieldSlice.GetFieldValue(T.x, T.y, fldB1);
+ fldB1.Combine(fB[ista], w[ista]);
- fldZ2 = z[i + 2];
+ fldZ2 = z[ista + 2];
dz = fldZ2 - fldZ1;
- sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fldB2);
- fldB2.Combine(fB[i + 2], w[i + 2]);
+ sta[ista].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fldB2);
+ fldB2.Combine(fB[ista + 2], w[ista + 2]);
fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
- for (++i; i < nHits; i++) {
- fldZ0 = z[i];
+ for (++ista; ista < nStations; ista++) {
+ fldZ0 = z[ista];
dz = (fldZ1 - fldZ0);
- sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fldB0);
- fldB0.Combine(fB[i], w[i]);
+ sta[ista].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fldB0);
+ fldB0.Combine(fB[ista], w[ista]);
fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
- fmask initialised = (z[i] <= z_end) & (z_start < z[i]);
- fvec w1 = iif(initialised, w[i], fvec::Zero());
- fvec w1_time = iif(initialised, w_time[i], fvec::Zero());
+ fmask initialised = (z[ista] <= z_end) & (z_start < z[ista]);
+ fvec w1 = iif(initialised, w[ista], fvec::Zero());
+ fvec w1_time = iif(initialised, w_time[ista], fvec::Zero());
fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
- L1Extrapolate(T, z[i], qp0, fld, &w1);
+ L1Extrapolate(T, z[ista], qp0, fld, &w1);
- // L1ExtrapolateLine( T, z[i]);
+ // L1ExtrapolateLine( T, z[ista]);
- T1.Extrapolate(z[i], qp0, fld, &w1);
+ T1.Extrapolate(z[ista], qp0, fld, w1);
- // T1.ExtrapolateLine( z[i]);
+ // T1.ExtrapolateLine( z[ista]);
- if (i == fNstationsBeforePipe) {
+ if (ista == fNstationsBeforePipe) {
fit.L1AddPipeMaterial(T, qp0, wIn);
fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
@@ -713,43 +754,44 @@ void L1Algo::L1KFTrackFitter()
T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(-1.f), wIn);
}
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, wIn);
- fit.EnergyLossCorrection(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, fvec(-1.f), wIn);
+ fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, wIn);
+ fit.EnergyLossCorrection(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec(-1.f), wIn);
- T1.L1AddMaterial(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, wIn);
- T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, T1.fx, T1.fy), qp01, fvec(-1.f), wIn);
+ T1.L1AddMaterial(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, wIn);
+ T1.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, T1.fx, T1.fy), qp01, fvec(-1.f), wIn);
}
else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
+ fit.L1AddMaterial(T, sta[ista].materialInfo, qp0, wIn);
}
- L1Filter(T, sta[i].frontInfo, u[i], w1);
- L1Filter(T, sta[i].backInfo, v[i], w1);
+ L1Filter(T, sta[ista].frontInfo, u[ista], w1);
+ L1Filter(T, sta[ista].backInfo, v[ista], w1);
- L1UMeasurementInfo info = sta[i].frontInfo;
- info.sigma2 = d_u[i] * d_u[i];
+ L1UMeasurementInfo info = sta[ista].frontInfo;
+ info.sigma2 = d_u[ista] * d_u[ista];
- T1.Filter(info, u[i], w1);
+ T1.Filter(info, u[ista], w1);
- info = sta[i].backInfo;
- info.sigma2 = d_v[i] * d_v[i];
+ info = sta[ista].backInfo;
+ info.sigma2 = d_v[ista] * d_v[ista];
- T1.Filter(info, v[i], w1);
+ T1.Filter(info, v[ista], w1);
- T1.FilterTime(time[i], timeEr[i], w1_time, sta[i].timeInfo);
+ T1.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo);
fldB2 = fldB1;
fldZ2 = fldZ1;
fldB1 = fldB0;
fldZ1 = fldZ0;
}
- // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
+ // fit.L1AddHalfMaterial( T, sta[ista].materialInfo, qp0 );
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
+ for (int iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->TLast[0] = T1.fx[iVec];
t[iVec]->TLast[1] = T1.fy[iVec];
t[iVec]->TLast[2] = T1.ftx[iVec];
t[iVec]->TLast[3] = T1.fty[iVec];
t[iVec]->TLast[4] = T1.fqp[iVec];
+ if (kGlobal == fTrackingMode) { t[iVec]->TLast[4] = T1pv.fqp[iVec]; }
t[iVec]->TLast[5] = T1.fz[iVec];
t[iVec]->TLast[6] = T1.ft[iVec];
@@ -1018,7 +1060,7 @@ void L1Algo::L1KFTrackFitterMuch()
L1Extrapolate(T, z[i], qp0, fld, &w1);
- T1.Extrapolate(z[i], qp01, fld, &w1);
+ T1.Extrapolate(z[i], qp01, fld, w1);
if (i == fNstationsBeforePipe) {
fit.L1AddPipeMaterial(T, qp0, wIn);
@@ -1308,7 +1350,7 @@ void L1Algo::L1KFTrackFitterMuch()
// fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01));
- T1.Extrapolate(z[i], qp01, fld, &w1);
+ T1.Extrapolate(z[i], qp01, fld, w1);
// L1Extrapolate( T, z[i], qp0, fld,&w1 );
@@ -1340,7 +1382,6 @@ void L1Algo::L1KFTrackFitterMuch()
}
// fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
-
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
t[iVec]->TFirst[0] = T1.fx[iVec];
@@ -1377,40 +1418,39 @@ void L1Algo::L1KFTrackFitterMuch()
t[iVec]->NDF = (int) T1.NDF[iVec];
}
-
// extrapolate to the PV region
- L1TrackParFit T1PV = T1;
- T1PV.Extrapolate(0.f, T1PV.fqp, fld);
+ L1TrackParFit T1pv = T1;
+ T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fqp, fld, fvec(1.));
for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec]->Tpv[0] = T1PV.fx[iVec];
- t[iVec]->Tpv[1] = T1PV.fy[iVec];
- t[iVec]->Tpv[2] = T1PV.ftx[iVec];
- t[iVec]->Tpv[3] = T1PV.fty[iVec];
- t[iVec]->Tpv[4] = T1PV.fqp[iVec];
- t[iVec]->Tpv[5] = T1PV.fz[iVec];
- t[iVec]->Tpv[6] = T1PV.ft[iVec];
-
- t[iVec]->Cpv[0] = T1PV.C00[iVec];
- t[iVec]->Cpv[1] = T1PV.C10[iVec];
- t[iVec]->Cpv[2] = T1PV.C11[iVec];
- t[iVec]->Cpv[3] = T1PV.C20[iVec];
- t[iVec]->Cpv[4] = T1PV.C21[iVec];
- t[iVec]->Cpv[5] = T1PV.C22[iVec];
- t[iVec]->Cpv[6] = T1PV.C30[iVec];
- t[iVec]->Cpv[7] = T1PV.C31[iVec];
- t[iVec]->Cpv[8] = T1PV.C32[iVec];
- t[iVec]->Cpv[9] = T1PV.C33[iVec];
- t[iVec]->Cpv[10] = T1PV.C40[iVec];
- t[iVec]->Cpv[11] = T1PV.C41[iVec];
- t[iVec]->Cpv[12] = T1PV.C42[iVec];
- t[iVec]->Cpv[13] = T1PV.C43[iVec];
- t[iVec]->Cpv[14] = T1PV.C44[iVec];
- t[iVec]->Cpv[15] = T1PV.C50[iVec];
- t[iVec]->Cpv[16] = T1PV.C51[iVec];
- t[iVec]->Cpv[17] = T1PV.C52[iVec];
- t[iVec]->Cpv[18] = T1PV.C53[iVec];
- t[iVec]->Cpv[19] = T1PV.C54[iVec];
- t[iVec]->Cpv[20] = T1PV.C55[iVec];
+ t[iVec]->Tpv[0] = T1pv.fx[iVec];
+ t[iVec]->Tpv[1] = T1pv.fy[iVec];
+ t[iVec]->Tpv[2] = T1pv.ftx[iVec];
+ t[iVec]->Tpv[3] = T1pv.fty[iVec];
+ t[iVec]->Tpv[4] = T1pv.fqp[iVec];
+ t[iVec]->Tpv[5] = T1pv.fz[iVec];
+ t[iVec]->Tpv[6] = T1pv.ft[iVec];
+
+ t[iVec]->Cpv[0] = T1pv.C00[iVec];
+ t[iVec]->Cpv[1] = T1pv.C10[iVec];
+ t[iVec]->Cpv[2] = T1pv.C11[iVec];
+ t[iVec]->Cpv[3] = T1pv.C20[iVec];
+ t[iVec]->Cpv[4] = T1pv.C21[iVec];
+ t[iVec]->Cpv[5] = T1pv.C22[iVec];
+ t[iVec]->Cpv[6] = T1pv.C30[iVec];
+ t[iVec]->Cpv[7] = T1pv.C31[iVec];
+ t[iVec]->Cpv[8] = T1pv.C32[iVec];
+ t[iVec]->Cpv[9] = T1pv.C33[iVec];
+ t[iVec]->Cpv[10] = T1pv.C40[iVec];
+ t[iVec]->Cpv[11] = T1pv.C41[iVec];
+ t[iVec]->Cpv[12] = T1pv.C42[iVec];
+ t[iVec]->Cpv[13] = T1pv.C43[iVec];
+ t[iVec]->Cpv[14] = T1pv.C44[iVec];
+ t[iVec]->Cpv[15] = T1pv.C50[iVec];
+ t[iVec]->Cpv[16] = T1pv.C51[iVec];
+ t[iVec]->Cpv[17] = T1pv.C52[iVec];
+ t[iVec]->Cpv[18] = T1pv.C53[iVec];
+ t[iVec]->Cpv[19] = T1pv.C54[iVec];
+ t[iVec]->Cpv[20] = T1pv.C55[iVec];
}
} // iter
}
@@ -1491,12 +1531,14 @@ void L1Algo::GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy,
fvec A0, A1 = ZERO, A2 = ZERO, A3 = ZERO, A4 = ZERO, A5 = ZERO, a0, a1 = ZERO, a2 = ZERO, b0, b1 = ZERO, b2 = ZERO;
fvec z0, x, y, z, S, w, wz, wS, time;
- time = 0;
+ time = fvec(0.);
int i = NHits - 1;
- if (zCur) z0 = *zCur;
- else
+ if (zCur) { z0 = *zCur; }
+ else {
z0 = zV[i];
+ }
+
w = wV[i];
A0 = w;
a0 = w * xV[i];
@@ -1524,9 +1566,13 @@ void L1Algo::GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy,
b0 += w * y;
b1 += wz * y;
b2 += wS * y;
- if (timeV) time += w_time[i] * (timeV[i] - sqrt(x * x + y * y + z * z) / 30.f);
+ if (timeV) {
+ fvec dx = xV[i] - fParameters.GetTargetPositionX();
+ fvec dy = yV[i] - fParameters.GetTargetPositionY();
+ fvec dz = zV[i] - fParameters.GetTargetPositionZ();
+ time += w_time[i] * (timeV[i] - sqrt(dx * dx + dy * dy + dz * dz) / fvec(30.));
+ }
}
- time = time / nhits;
fvec A3A3 = A3 * A3;
fvec A3A4 = A3 * A4;
@@ -1556,7 +1602,7 @@ void L1Algo::GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy,
t.fy = (A2 * b0 - A1 * b1) * det;
t.fty = (-A1 * b0 + A0 * b1) * det;
t.fqp = -L * c_light_i / sqrt(txtx1 + t.fty * t.fty);
- if (timeV) { t.ft = iif(nhits > fvec::Zero(), time, fvec::Zero()); }
+ if (timeV) { t.ft = iif(nhits > fvec::Zero(), time / nhits, fvec::Zero()); }
t.fz = z0;
}
@@ -1579,6 +1625,11 @@ void L1Algo::FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st)
track.C42 = ZERO;
track.C43 = ZERO;
track.C44 = ONE;
+ track.C50 = ZERO;
+ track.C51 = ZERO;
+ track.C52 = ZERO;
+ track.C53 = ZERO;
+ track.C55 = 2.6f * 2.6f;
track.x = x;
track.y = y;
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 75b1d4ec2b03a4b3bb361ce166e75fde43009142..b0c922bc9151a27dbca7fa780bc21a990a1ef6fa 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -149,7 +149,7 @@ inline void L1TrackPar::Print(int i) const
std::cout << " c55 " << C55[i] << std::endl;
}
PrintCorrelations(i);
- std::cout.flags( coutFlags );
+ std::cout.flags(coutFlags);
}
inline void L1TrackPar::PrintCorrelations(int i) const
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index 08a3c946d66e97a109701eb928ba43c7f637adf1..612265343590aa63b781754daafccdb9ebcf1bdc 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -37,7 +37,7 @@ void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w)
// correction to HCH is needed for the case when sigma2 is so small
// with respect to HCH that it disappears due to the roundoff error
//
- fvec wi = w / (info.sigma2 + 1.0000001f * HCH);
+ fvec wi = w / (info.sigma2 + fvec(1.0000001) * HCH);
fvec zetawi = w * zeta / (iif(maskDoFilter, info.sigma2, fvec::Zero()) + HCH);
// chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() );
@@ -395,7 +395,21 @@ void L1TrackParFit::ExtrapolateLine1(fvec z_out, fvec* w, fvec v)
void L1TrackParFit::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(fvec z_out, // extrapolate to this z position
fvec qp0, // use Q/p linearisation at this value
- const L1FieldRegion& F, fvec* w)
+ const L1FieldRegion& F, const fvec& w)
+{
+ fvec sgn = iif(fz < z_out, fvec(1.), fvec(-1.));
+ while (!(w * abs(z_out - fz) <= fvec(1.e-6)).isFull()) {
+ fvec zNew = fz + sgn * fvec(50.); // max. 50 cm step
+ zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out;
+ ExtrapolateStep(zNew, qp0, F, w);
+ }
+}
+
+void
+ L1TrackParFit::ExtrapolateStep // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
+ (fvec z_out, // extrapolate to this z position
+ fvec qp0, // use Q/p linearisation at this value
+ const L1FieldRegion& F, const fvec& w)
{
//
// Forth-order Runge-Kutta method for solution of the equation
@@ -424,20 +438,18 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja
cnst c_light = 0.000299792458;
- const fvec a[4] = {0.0f, 0.5f, 0.5f, 1.0f};
- const fvec c[4] = {1.0f / 6.0f, 1.0f / 3.0f, 1.0f / 3.0f, 1.0f / 6.0f};
- const fvec b[4] = {0.0f, 0.5f, 0.5f, 1.0f};
+ const fvec a[4] = {0., 0.5, 0.5, 1.};
+ const fvec c[4] = {1. / 6., 1. / 3., 1. / 3., 1. / 6.};
+ const fvec b[4] = {0., 0.5, 0.5, 1.};
- int step4;
- fvec k[20], x0[5], x[5], k1[20];
+ fvec k[20], x[5], k1[20];
fvec Ax[4], Ay[4], Ax_tx[4], Ay_tx[4], Ax_ty[4], Ay_ty[4], At[4], At_tx[4], At_ty[4];
//----------------------------------------------------------------
- if (w) { z_out = iif((fvec(0.f) < *w), z_out, fz); }
+ z_out = iif((fvec(0.f) < w), z_out, fz);
fvec qp_in = fqp;
- const fvec z_in = fz;
const fvec h = (z_out - fz);
// cout<<h<<" h"<<endl;
@@ -445,54 +457,54 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja
// cout<<fty<<" fty"<<endl;
fvec hC = h * c_light;
- x0[0] = fx;
- x0[1] = fy;
- x0[2] = ftx;
- x0[3] = fty;
- x0[4] = ft;
+ //std::cout << "fx " << fx << std::endl;
+ fvec x0[5];
+ x0[0] = fx;
+ x0[1] = fy;
+ x0[2] = ftx;
+ x0[3] = fty;
+ x0[4] = ft;
+
//
// Runge-Kutta step
//
- int step;
- int i;
+ for (int step = 0; step < 4; ++step) {
- fvec B[4][3];
- for (step = 0; step < 4; ++step) {
- F.Get(z_in + a[step] * h, B[step]);
- //std::cout << "extrapolation step " << step << " z " << z_in + a[step] * h << " field " << B[step][0] << " " << B[step][1] << " "
- // << B[step][2] << std::endl;
- }
-
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 5; ++i) {
+ for (int i = 0; i < 5; ++i) {
if (step == 0) { x[i] = x0[i]; }
else {
x[i] = x0[i] + b[step] * k[step * 5 - 5 + i];
}
}
+ fvec B[3];
+
+ F.Get(x[0], x[1], fz + a[step] * h, B);
+ //std::cout << "extrapolation step " << step << " z " << z_in + a[step] * h << " field " << B[0] << " " << B[1] << " "
+ // << B[2] << std::endl;
+
fvec tx = x[2];
fvec ty = x[3];
fvec tx2 = tx * tx;
fvec ty2 = ty * ty;
fvec txty = tx * ty;
- fvec tx2ty21 = 1.f + tx2 + ty2;
+ fvec tx2ty21 = fvec(1.) + tx2 + ty2;
// if( tx2ty21 > 1.e4 ) return 1;
- fvec I_tx2ty21 = 1.f / tx2ty21 * qp0;
+ fvec I_tx2ty21 = fvec(1.) / tx2ty21 * qp0;
fvec tx2ty2 = sqrt(tx2ty21);
// fvec I_tx2ty2 = qp0 * hC / tx2ty2 ; unsused ???
tx2ty2 *= hC;
fvec tx2ty2qp = tx2ty2 * qp0;
- // cout<<B[step][0]<<" B["<<step<<"][0] "<<B[step][2]<<" B["<<step<<"][2] "<<B[step][1]<<" B["<<step<<"][1]"<<endl;
- Ax[step] = (txty * B[step][0] + ty * B[step][2] - (1.f + tx2) * B[step][1]) * tx2ty2;
- Ay[step] = (-txty * B[step][1] - tx * B[step][2] + (1.f + ty2) * B[step][0]) * tx2ty2;
+ // cout<<B[0]<<" B["<<step<<"][0] "<<B[2]<<" B["<<step<<"][2] "<<B[1]<<" B["<<step<<"][1]"<<endl;
+ Ax[step] = (txty * B[0] + ty * B[2] - (fvec(1.) + tx2) * B[1]) * tx2ty2;
+ Ay[step] = (-txty * B[1] - tx * B[2] + (fvec(1.) + ty2) * B[0]) * tx2ty2;
- Ax_tx[step] = Ax[step] * tx * I_tx2ty21 + (ty * B[step][0] - 2.f * tx * B[step][1]) * tx2ty2qp;
- Ax_ty[step] = Ax[step] * ty * I_tx2ty21 + (tx * B[step][0] + B[step][2]) * tx2ty2qp;
- Ay_tx[step] = Ay[step] * tx * I_tx2ty21 + (-ty * B[step][1] - B[step][2]) * tx2ty2qp;
- Ay_ty[step] = Ay[step] * ty * I_tx2ty21 + (-tx * B[step][1] + 2.f * ty * B[step][0]) * tx2ty2qp;
+ Ax_tx[step] = Ax[step] * tx * I_tx2ty21 + (ty * B[0] - fvec(2.) * tx * B[1]) * tx2ty2qp;
+ Ax_ty[step] = Ax[step] * ty * I_tx2ty21 + (tx * B[0] + B[2]) * tx2ty2qp;
+ Ay_tx[step] = Ay[step] * tx * I_tx2ty21 + (-ty * B[1] - B[2]) * tx2ty2qp;
+ Ay_ty[step] = Ay[step] * ty * I_tx2ty21 + (-tx * B[1] + fvec(2.) * ty * B[0]) * tx2ty2qp;
fvec m2 = fMass2;
fvec vi = sqrt(fvec(1.) + m2 * qp0 * qp0) / fvec(29.9792458f);
@@ -503,12 +515,13 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja
// cout<<Ax[step]<<" Ax[step] "<<Ay[step]<<" ay "<<At[step]<<" At[step] "<<qp0<<" qp0 "<<h<<" h"<<endl;
- step4 = step * 5;
- k[step4] = tx * h;
- k[step4 + 1] = ty * h;
- k[step4 + 2] = Ax[step] * qp0;
- k[step4 + 3] = Ay[step] * qp0;
- k[step4 + 4] = At[step];
+ int step5 = step * 5;
+
+ k[step5 + 0] = tx * h;
+ k[step5 + 1] = ty * h;
+ k[step5 + 2] = Ax[step] * qp0;
+ k[step5 + 3] = Ay[step] * qp0;
+ k[step5 + 4] = At[step];
} // end of Runge-Kutta steps
{
@@ -520,6 +533,9 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja
// cout << "initialised = " << initialised << "; ";
// cout << "fx = " << fx;
+ //std::cout << " x : x0[0] " << x0[0] << " k[0] " << k[0] << " k[5] " << k[5] << " k[10] " << k[10] << " k[15] "
+ // << k[15] << std::endl;
+
fx = x0[0] + c[0] * k[0] + c[1] * k[5 + 0] + c[2] * k[10 + 0] + c[3] * k[15 + 0];
fy = x0[1] + c[0] * k[1] + c[1] * k[5 + 1] + c[2] * k[10 + 1] + c[3] * k[15 + 1];
ftx = x0[2] + c[0] * k[2] + c[1] * k[5 + 2] + c[2] * k[10 + 2] + c[3] * k[15 + 2];
@@ -531,133 +547,149 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja
}
// cout<<fx<<" fx"<<endl;
- //
- // Derivatives dx/dqp
- //
-
- x0[0] = 0.f;
- x0[1] = 0.f;
- x0[2] = 0.f;
- x0[3] = 0.f;
- x0[4] = 0.f;
+ fvec J[36]; // Jacobian of extrapolation
//
- // Runge-Kutta step for derivatives dx/dqp
-
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 5; ++i) {
- if (step == 0) { x[i] = x0[i]; }
- else {
- x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
- }
- }
- step4 = step * 5;
- k1[step4] = x[2] * h;
- k1[step4 + 1] = x[3] * h;
- k1[step4 + 2] = Ax[step] + Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
- k1[step4 + 3] = Ay[step] + Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
- k1[step4 + 4] = At[step] + At_tx[step] * x[2] + At_ty[step] * x[3];
+ // derivatives dx/dx and dx/dy
+ //
- } // end of Runge-Kutta steps for derivatives dx/dqp
+ for (int i = 0; i < 12; ++i) {
+ J[i] = fvec(0.);
+ }
- fvec J[36];
+ J[0] = fvec(1.);
+ J[7] = fvec(1.);
- for (i = 0; i < 4; ++i) {
- J[24 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i] + c[3] * k1[15 + i];
- }
- J[28] = 1.;
- J[29] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4] + c[3] * k1[15 + 4];
- //
- // end of derivatives dx/dqp
//
-
// Derivatives dx/tx
//
- x0[0] = 0.f;
- x0[1] = 0.f;
- x0[2] = 1.f;
- x0[3] = 0.f;
- x0[4] = 0.f;
+ x0[0] = fvec(0.);
+ x0[1] = fvec(0.);
+ x0[2] = fvec(1.);
+ x0[3] = fvec(0.);
+ x0[4] = fvec(0.);
//
// Runge-Kutta step for derivatives dx/dtx
//
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 5; ++i) {
+ for (int step = 0; step < 4; ++step) {
+ int step5 = step * 5;
+ for (int i = 0; i < 5; ++i) {
if (step == 0) { x[i] = x0[i]; }
else if (i != 2) {
- x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
+ x[i] = x0[i] + b[step] * k1[step5 - 5 + i];
}
}
- step4 = step * 5;
- k1[step4] = x[2] * h;
- k1[step4 + 1] = x[3] * h;
- // k1[step4+2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
- k1[step4 + 3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
- k1[step4 + 4] = At_tx[step] * x[2] + At_ty[step] * x[3];
+
+ k1[step5 + 0] = x[2] * h;
+ k1[step5 + 1] = x[3] * h;
+ // k1[step5+2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
+ k1[step5 + 3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
+ k1[step5 + 4] = At_tx[step] * x[2] + At_ty[step] * x[3];
} // end of Runge-Kutta steps for derivatives dx/dtx
- for (i = 0; i < 4; ++i) {
+ for (int i = 0; i < 4; ++i) {
if (i != 2) { J[12 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i] + c[3] * k1[15 + i]; }
}
- // end of derivatives dx/dtx
- J[14] = 1.f;
- J[16] = 0.f;
+
+ J[14] = fvec(1.);
+ J[16] = fvec(0.);
J[17] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4] + c[3] * k1[15 + 4];
+ // end of derivatives dx/dtx
+
+ //
// Derivatives dx/ty
//
- x0[0] = 0.f;
- x0[1] = 0.f;
- x0[2] = 0.f;
- x0[3] = 1.f;
- x0[4] = 0.f;
+ x0[0] = fvec(0.);
+ x0[1] = fvec(0.);
+ x0[2] = fvec(0.);
+ x0[3] = fvec(1.);
+ x0[4] = fvec(0.);
//
// Runge-Kutta step for derivatives dx/dty
//
- for (step = 0; step < 4; ++step) {
- for (i = 0; i < 5; ++i) {
+ for (int step = 0; step < 4; ++step) {
+ int step5 = step * 5;
+ for (int i = 0; i < 5; ++i) {
if (step == 0) {
x[i] = x0[i]; // ty fixed
}
else if (i != 3) {
- x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
+ x[i] = x0[i] + b[step] * k1[step5 - 5 + i];
}
}
- step4 = step * 5;
- k1[step4] = x[2] * h;
- k1[step4 + 1] = x[3] * h;
- k1[step4 + 2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
- // k1[step4+3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3]; // TODO: SG: check if the simplification below is ok
- k1[step4 + 4] = At_tx[step] * x[2] + At_ty[step] * x[3];
+ k1[step5 + 0] = x[2] * h;
+ k1[step5 + 1] = x[3] * h;
+ k1[step5 + 2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
+ // k1[step5+3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3]; // TODO: SG: check if the simplification below is ok
+ k1[step5 + 4] = At_tx[step] * x[2] + At_ty[step] * x[3];
} // end of Runge-Kutta steps for derivatives dx/dty
- for (i = 0; i < 3; ++i) {
+ for (int i = 0; i < 3; ++i) {
J[18 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i] + c[3] * k1[15 + i];
}
- // end of derivatives dx/dty
+
J[21] = fvec(1.);
J[22] = fvec(0.);
J[23] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4] + c[3] * k1[15 + 4];
+
+ // end of derivatives dx/dty
+
+ //
+ // Derivatives dx/dqp
//
- // derivatives dx/dx and dx/dy
- for (i = 0; i < 12; ++i) {
- J[i] = fvec(0.);
+ x0[0] = fvec(0.);
+ x0[1] = fvec(0.);
+ x0[2] = fvec(0.);
+ x0[3] = fvec(0.);
+ x0[4] = fvec(0.);
+
+ //
+ // Runge-Kutta step for derivatives dx/dqp
+ //
+
+ for (int step = 0; step < 4; ++step) {
+ for (int i = 0; i < 5; ++i) {
+ if (step == 0) { x[i] = x0[i]; }
+ else {
+ x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
+ }
+ }
+ int step5 = step * 5;
+ k1[step5 + 0] = x[2] * h;
+ k1[step5 + 1] = x[3] * h;
+ k1[step5 + 2] = Ax[step] + Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
+ k1[step5 + 3] = Ay[step] + Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
+ k1[step5 + 4] = At[step] + At_tx[step] * x[2] + At_ty[step] * x[3];
+
+ } // end of Runge-Kutta steps for derivatives dx/dqp
+
+ for (int i = 0; i < 4; ++i) {
+ J[24 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i] + c[3] * k1[15 + i];
}
+ J[28] = fvec(1.);
+ J[29] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4] + c[3] * k1[15 + 4];
- J[0] = fvec(1.);
- J[7] = fvec(1.);
- for (i = 30; i < 35; i++) {
+ // end of derivatives dx/dqp
+
+ //
+ // derivatives dx/dt
+ //
+
+ for (int i = 30; i < 35; i++) {
J[i] = fvec(0.);
}
J[35] = fvec(1.);
+ // end of derivatives dx/dt
+
fvec dqp = qp_in - qp0;
{ // update parameters
@@ -673,20 +705,26 @@ void L1TrackParFit::Extrapolate // extrapolates track parameters and returns ja
//cout<< (ft - ft_old)<<" ft dt "<<endl;
// // if ( C_in&&C_out ) CbmKFMath::multQtSQ( 5, J, C_in, C_out); // TODO
- // j(0,2) = J[5*2 + 0];
- // j(1,2) = J[5*2 + 1];
- // j(2,2) = J[5*2 + 2];
- // j(3,2) = J[5*2 + 3];
+ // j(0,2) = J[6*2 + 0];
+ // j(1,2) = J[6*2 + 1];
+ // j(2,2) = J[6*2 + 2];
+ // j(3,2) = J[6*2 + 3];
+ // j(4,2) = J[6*2 + 4];
+ // j(5,2) = J[6*2 + 5];
//
- // j(0,3) = J[5*3 + 0];
- // j(1,3) = J[5*3 + 1];
- // j(2,3) = J[5*3 + 2];
- // j(3,3) = J[5*3 + 3];
+ // j(0,3) = J[6*3 + 0];
+ // j(1,3) = J[6*3 + 1];
+ // j(2,3) = J[6*3 + 2];
+ // j(3,3) = J[6*3 + 3];
+ // j(4,3) = J[6*3 + 4];
+ // j(5,3) = J[6*3 + 5];
//
- // j(0,4) = J[5*4 + 0];
- // j(1,4) = J[5*4 + 1];
- // j(2,4) = J[5*4 + 2];
- // j(3,4) = J[5*4 + 3];
+ // j(0,4) = J[6*4 + 0];
+ // j(1,4) = J[6*4 + 1];
+ // j(2,4) = J[6*4 + 2];
+ // j(3,4) = J[6*4 + 3];
+ // j(4,4) = J[6*4 + 4];
+ // j(5,4) = J[6*4 + 5];
const fvec c42 = C42, c43 = C43;
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index 7350d34c4cd0fae9289fff7ab1d95e21d529d9a9..907ce9a9815030f64a246e044ec67cd877e489a0 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -109,9 +109,13 @@ public:
void FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y);
void FilterTime(fvec t0, fvec dt0, fvec w = 1., fvec timeInfo = 1.);
void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
- void Extrapolate(fvec z_out, fvec qp0, const L1FieldRegion& F, fvec* w = 0);
+
+ void Extrapolate(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);
+ void ExtrapolateStep(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);
+
void ExtrapolateLine(fvec z_out, fvec* w = 0);
void ExtrapolateLine1(fvec z_out, fvec* w = 0, fvec v = 0);
+
void Compare(L1TrackPar& T);
void EnergyLossCorrection(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
@@ -146,16 +150,25 @@ public:
inline void L1TrackParFit::Print(int i)
{
+ std::ios_base::fmtflags coutFlags(std::cout.flags());
std::cout.setf(std::ios::scientific, std::ios::floatfield);
+
if (i == -1) {
std::cout << "T = " << std::endl;
- std::cout << fx << std::endl;
- std::cout << fy << std::endl;
- std::cout << ftx << std::endl;
- std::cout << fty << std::endl;
- std::cout << fqp << std::endl;
- std::cout << fz << std::endl;
- std::cout << ft << std::endl;
+ std::cout << " x " << fx << std::endl;
+ std::cout << " y " << fy << std::endl;
+ std::cout << " tx " << ftx << std::endl;
+ std::cout << " ty " << fty << std::endl;
+ std::cout << " qp " << fqp << std::endl;
+ std::cout << " t " << ft << std::endl;
+ std::cout << " z " << fz << std::endl;
+ std::cout << "C = " << std::endl;
+ std::cout << " c00 " << C00 << std::endl;
+ std::cout << " c11 " << C11 << std::endl;
+ std::cout << " c22 " << C22 << std::endl;
+ std::cout << " c33 " << C33 << std::endl;
+ std::cout << " c44 " << C44 << std::endl;
+ std::cout << " c55 " << C55 << std::endl;
}
else {
std::cout << "T = ";
@@ -173,6 +186,7 @@ inline void L1TrackParFit::Print(int i)
std::cout << C44[i] << " ";
std::cout << C55[i] << std::endl;
}
+ std::cout.flags(coutFlags);
}
diff --git a/reco/qa/CbmTrackingTrdQa.cxx b/reco/qa/CbmTrackingTrdQa.cxx
index 7a261e425246ab8e1d16c1327d3d142b4f07e14d..078e179bdead36898d844edfe51339ef6fa1ec7b 100644
--- a/reco/qa/CbmTrackingTrdQa.cxx
+++ b/reco/qa/CbmTrackingTrdQa.cxx
@@ -262,6 +262,8 @@ void CbmTrackingTrdQa::Exec(Option_t* /*opt*/)
// Bool_t isRec = kFALSE;
if (globalTrackId < 0) continue;
if (trdTrackId < 0) continue;
+ CbmGlobalTrack* globalTrack = (CbmGlobalTrack*) fGlobalTracks->At(globalTrackId);
+
CbmTrdTrack* trdTrack = (CbmTrdTrack*) fTrdTracks->At(trdTrackId);
assert(trdTrack);
assert(quali >= fQuota);
@@ -273,6 +275,14 @@ void CbmTrackingTrdQa::Exec(Option_t* /*opt*/)
Int_t nFake = 0;
Int_t nAllHits = trdTrack->GetNofHits();
assert(nTrue + nWrong + nFake == nAllHits);
+
+ double qp = globalTrack->GetParamFirst()->GetQp();
+ //double q = (qp >= 1.) ? 1. : -1.;
+ double p = (fabs(qp) > 1. / 100.) ? 1. / fabs(qp) : 100.;
+ double tx = globalTrack->GetParamFirst()->GetTx();
+ double ty = globalTrack->GetParamFirst()->GetTy();
+ double pt = sqrt((tx * tx + ty * ty) / (1. + tx * tx + ty * ty)) * p;
+
// Verbose output
LOG(debug1) << GetName() << ": MCTrack " << iMcTrack << ", stations " << nStations << ", hits " << nAllHits
<< ", true hits " << nTrue;
@@ -287,6 +297,7 @@ void CbmTrackingTrdQa::Exec(Option_t* /*opt*/)
fhPtRecPrim->Fill(info.fPt);
fhPidPtY["prmE"][Pdg2Idx(info.fPdg)]->Fill(info.fY, info.fPt);
fhNpRecPrim->Fill(Double_t(nAllHits));
+ if (info.fPt > 0.001) { fhPtResPrim->Fill((pt / info.fPt - 1.)); }
}
else {
nRecSec++;
@@ -715,6 +726,9 @@ void CbmTrackingTrdQa::CreateHistos()
fHistoList->Add(fhNhClones);
fHistoList->Add(fhNhGhosts);
+ fhPtResPrim = new TH1F("hPtPrim", "Resolution Pt Primaries [100%]", 100, -1., 1.);
+ fHistoList->Add(fhPtResPrim);
+
TIter next(fHistoList);
while (TH1* histo = ((TH1*) next())) {
fOutFolder.Add(histo);
diff --git a/reco/qa/CbmTrackingTrdQa.h b/reco/qa/CbmTrackingTrdQa.h
index b6efeb075ca76258de1d03aa277bfa62b61a0dbf..f2d1414f58a53e21ae021b8c5d3a2f7678a1f2a9 100644
--- a/reco/qa/CbmTrackingTrdQa.h
+++ b/reco/qa/CbmTrackingTrdQa.h
@@ -214,6 +214,8 @@ private:
TH1F* fhNhClones = nullptr;
TH1F* fhNhGhosts = nullptr;
+ // Pt resolution
+ TH1F* fhPtResPrim = nullptr;
/** List of histograms **/
TList* fHistoList = nullptr;