diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 68e7d8df887c325c0582dd96387230b1af70397e..10863478a5c0bff9725438569bdeb70b7c2f0331 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -335,13 +335,11 @@ InitStatus CbmL1::Init()
if (nullptr == fpStsPoints) LOG(fatal) << GetName() << ": No StsPoint data!";
}
- if (!fUseTRD) {
- fpTrdPoints = 0;
- fpTrdHitMatches = 0;
- }
+ fpTrdPoints = mcManager->InitBranch("TrdPoint");
+
+ if (!fUseTRD) { fpTrdHitMatches = 0; }
else {
fpTrdHitMatches = (TClonesArray*) fairManager->GetObject("TrdHitMatch");
- fpTrdPoints = mcManager->InitBranch("TrdPoint");
}
if (!fUseMUCH) {
@@ -349,7 +347,6 @@ InitStatus CbmL1::Init()
fpMuchHitMatches = 0;
}
else {
-
fpMuchDigis = (TClonesArray*) fairManager->GetObject("MuchDigi");
fpMuchDigiMatches = (TClonesArray*) fairManager->GetObject("MuchDigiMatch");
fpMuchClusters = (TClonesArray*) fairManager->GetObject("MuchCluster");
@@ -357,12 +354,10 @@ InitStatus CbmL1::Init()
fpMuchHitMatches = L1_DYNAMIC_CAST<TClonesArray*>(fairManager->GetObject("MuchPixelHitMatch"));
}
- if (!fUseTOF) {
- fpTofPoints = 0;
- fpTofHitMatches = 0;
- }
+ fpTofPoints = mcManager->InitBranch("TofPoint");
+
+ if (!fUseTOF) { fpTofHitMatches = 0; }
else {
- fpTofPoints = mcManager->InitBranch("TofPoint");
fpTofHitMatches = static_cast<TClonesArray*>(fairManager->GetObject("TofHitMatch"));
}
}
@@ -1040,18 +1035,19 @@ void CbmL1::Reconstruct(CbmEvent* event)
for (L1Vector<L1Track>::iterator it = fpAlgo->fTracks.begin(); it != fpAlgo->fTracks.end(); it++) {
CbmL1Track t;
- for (int i = 0; i < 7; i++)
- t.T[i] = it->TFirst[i];
- for (int i = 0; i < 21; i++)
- t.C[i] = it->CFirst[i];
- for (int i = 0; i < 7; i++)
+
+ for (int i = 0; i < L1TrackPar::kNparTr; i++) {
+ t.T[i] = it->TFirst[i];
t.TLast[i] = it->TLast[i];
- for (int i = 0; i < 21; i++)
+ t.Tpv[i] = it->Tpv[i];
+ }
+
+ for (int i = 0; i < L1TrackPar::kNparCov; i++) {
+ t.C[i] = it->CFirst[i];
t.CLast[i] = it->CLast[i];
- for (int k = 0; k < 7; k++)
- t.Tpv[k] = it->Tpv[k];
- for (int k = 0; k < 21; k++)
- t.Cpv[k] = it->Cpv[k];
+ t.Cpv[i] = it->Cpv[i];
+ }
+
t.chi2 = it->chi2;
t.NDF = it->NDF;
//t.T[4] = it->Momentum;
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 7cf53c95b7ff63be02f29e42659677bc101da8b9..20995c11b8714546506c5ad27b26cfd5f6687ffe 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -420,6 +420,8 @@ private:
/// \note Should be called only after CbmL1::Performance()
void TrackFitPerformance();
+ void FillFitHistos(L1TrackPar& tr, const CbmL1MCPoint& mc, bool isTimeFitted, TH1F* h[]);
+
/// Fills performance histograms
void HistoPerformance();
@@ -485,7 +487,7 @@ public:
L1InitManager fInitManager; ///< Tracking parameters data manager
L1IODataManager fIODataManager; ///< Input-output data manager
- bool fMissingHits = false; ///< Turns on several ad-hock settings for "mcbm_beam_2021_07_surveyed.100ev" setup
+ bool fMissingHits = false; ///< Turns on several ad-hock settings for "mcbm_beam_2021_07_surveyed.100ev" setup
L1Algo::TrackingMode fTrackingMode = L1Algo::TrackingMode::kSts; ///< Tracking mode: kSts, kMcbm or kGlobal
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 4b9f07235f4d614daa978fc493d2a41161692b19..5a2479c3227d94e65a9e6478dd5d6021ceb5ccfc 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1077,8 +1077,8 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitFa
void CbmL1::TrackFitPerformance()
{
- const int Nh_fit = 14;
- static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit], *h_fitPV[Nh_fit],
+ const int Nh_fit = 17;
+ static TH1F *h_fit[Nh_fit], *h_fitL[Nh_fit], *h_fitSV[Nh_fit], *h_fitPV[Nh_fit], *h_fitTrd[Nh_fit], *h_fitTof[Nh_fit],
*h_fit_chi2; //, *h_smoothed[12][Nh_fit];
static TH2F *PRes2D, *PRes2DPrim, *PRes2DSec;
@@ -1104,6 +1104,8 @@ void CbmL1::TrackFitPerformance()
L1Fit fit;
fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
fit.SetMask(fmask::One());
+ //fit.SetMaxExtrapolationStep(10.);
+ fit.SetDoFitVelocity(true);
if (first_call) {
first_call = 0;
@@ -1139,12 +1141,8 @@ void CbmL1::TrackFitPerformance()
Double_t l, r;
} Table[Nh_fit] = {{"x", "Residual X [#mum]", 140, -40., 40.},
{"y", "Residual Y [#mum]", 100, -450., 450.},
- //{"y", "Residual Y [#mum]", 100, -45., 45.},
{"tx", "Residual Tx [mrad]", 100, -4., 4.},
- //{"tx", "Residual Tx [mrad]", 100, -7., 7.},
- //{"tx", "Residual Tx [mrad]", 100, -2.5, 2.5},
{"ty", "Residual Ty [mrad]", 100, -3.5, 3.5},
- //{"ty", "Residual Ty [mrad]", 100, -2.5, 2.5},
{"P", "Resolution P/Q [%]", 100, -15., 15.},
{"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
{"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
@@ -1154,7 +1152,10 @@ void CbmL1::TrackFitPerformance()
{"QPreco", "Reco Q/P ", 100, -10., 10.},
{"QPmc", "MC Q/P ", 100, -10., 10.},
{"time", "Residual Time [ns]", 100, -10., 10.},
- {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.}};
+ {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.},
+ {"vi", "Residual Vi [1/c]", 100, -4., 4.},
+ {"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
+ {"distrVi", "Vi distribution [1/c]", 100, 0., 4.}};
if (L1Algo::kGlobal == fpAlgo->fTrackingMode) {
Table[4].l = -1.;
@@ -1176,15 +1177,18 @@ void CbmL1::TrackFitPerformance()
//{"ty", "Residual Ty [mrad]", 100, -3., 3.},
{"ty", "Residual Ty [mrad]", 100, -2., 2.},
{"P", "Resolution P/Q [%]", 100, -15., 15.},
- {"px", "Pull X [residual/estimated_error]", 100, -6., 6.},
- {"py", "Pull Y [residual/estimated_error]", 100, -6., 6.},
- {"ptx", "Pull Tx [residual/estimated_error]", 100, -6., 6.},
- {"pty", "Pull Ty [residual/estimated_error]", 100, -6., 6.},
- {"pQP", "Pull Q/P [residual/estimated_error]", 100, -6., 6.},
+ {"px", "Pull X [residual/estimated_error]", 100, -10., 10.},
+ {"py", "Pull Y [residual/estimated_error]", 100, -10., 10.},
+ {"ptx", "Pull Tx [residual/estimated_error]", 100, -10., 10.},
+ {"pty", "Pull Ty [residual/estimated_error]", 100, -10., 10.},
+ {"pQP", "Pull Q/P [residual/estimated_error]", 100, -10., 10.},
{"QPreco", "Reco Q/P ", 100, -10., 10.},
{"QPmc", "MC Q/P ", 100, -10., 10.},
- {"t", "Residual T [ns]", 100, -10., 10.},
- {"pt", "Pull T [residual/estimated_error]", 100, -6., 6.}};
+ {"time", "Residual Time [ns]", 100, -10., 10.},
+ {"ptime", "Pull Time [residual/estimated_error]", 100, -10., 10.},
+ {"vi", "Residual Vi [1/c]", 100, -10., 10.},
+ {"pvi", "Pull Vi [residual/estimated_error]", 100, -10., 10.},
+ {"distrVi", "Vi distribution [1/c]", 100, -10., 10.}};
if (L1Algo::kGlobal == fpAlgo->fTrackingMode) {
TableVertex[4].l = -1.;
@@ -1206,6 +1210,26 @@ void CbmL1::TrackFitPerformance()
sprintf(t, "Primary vertex point %s", TableVertex[i].title);
h_fitPV[i] = new TH1F(n, t, TableVertex[i].n, TableVertex[i].l, TableVertex[i].r);
+ sprintf(n, "trd_%s", TableVertex[i].name);
+ sprintf(t, "Last Trd point %s", TableVertex[i].title);
+ if (i == 12) { h_fitTrd[i] = new TH1F(n, t, Table[i].n, -30., 30.); }
+ else if (i == 13) {
+ h_fitTrd[i] = new TH1F(n, t, Table[i].n, -5., 5.);
+ }
+ else {
+ h_fitTrd[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
+ }
+
+ sprintf(n, "tof_%s", TableVertex[i].name);
+ sprintf(t, "Last Tof point %s", TableVertex[i].title);
+ if (i == 12) { h_fitTof[i] = new TH1F(n, t, Table[i].n, -30., 30.); }
+ else if (i == 13) {
+ h_fitTof[i] = new TH1F(n, t, Table[i].n, -5., 5.);
+ }
+ else {
+ h_fitTof[i] = new TH1F(n, t, Table[i].n, Table[i].l, Table[i].r);
+ }
+
for (int j = 0; j < 12; j++) {
sprintf(n, "smth_%s_sta_%i", Table[i].name, j);
sprintf(t, "Station %i %s", j, Table[i].title);
@@ -1221,7 +1245,16 @@ void CbmL1::TrackFitPerformance()
if (it->IsGhost()) continue;
- const int isTimeFitted = (fvHitStore[it->Hits.back()].iStation > fNMvdStations);
+ bool isTimeFitted = 0;
+
+ {
+ int nTimeMeasurenments = 0;
+ for (unsigned int ih = 0; ih < it->Hits.size(); ih++) {
+ int ista = fvHitStore[it->Hits[ih]].iStation;
+ if (fpAlgo->GetParameters()->GetStation(ista).timeInfo) { nTimeMeasurenments++; }
+ }
+ isTimeFitted = (nTimeMeasurenments > 1);
+ }
do { // first hit
@@ -1239,12 +1272,9 @@ void CbmL1::TrackFitPerformance()
const CbmL1MCPoint& mcP = fvMCPoints[imcPoint];
- fit.SetTrack(it->T, it->C);
- const L1TrackPar& tr = fit.Tr();
-
- L1FieldRegion fld _fvecalignment;
- fld.SetUseOriginalField();
- fit.Extrapolate(mcP.zIn, fld);
+ L1TrackPar tr;
+ tr.copyFromArrays(it->T, it->C);
+ FillFitHistos(tr, mcP, isTimeFitted, h_fit);
double dx = tr.x[0] - mcP.xIn;
double dy = tr.y[0] - mcP.yIn;
@@ -1252,11 +1282,6 @@ void CbmL1::TrackFitPerformance()
// make dca distance negative for the half of the tracks to ease the gaussian fit and the rms calculation
if (mcTrack.ID % 2) dca = -dca;
double pt = sqrt(mcP.px * mcP.px + mcP.py * mcP.py);
- h_fit[0]->Fill(dx * 1.e4);
- h_fit[1]->Fill(dy * 1.e4);
- h_fit[2]->Fill((tr.tx[0] - mcP.pxIn / mcP.pzIn) * 1.e3);
- h_fit[3]->Fill((tr.ty[0] - mcP.pyIn / mcP.pzIn) * 1.e3);
-
double dP = (fabs(1. / tr.qp[0]) - mcP.p) / mcP.p;
PRes2D->Fill(mcP.p, 100. * dP);
@@ -1283,88 +1308,23 @@ void CbmL1::TrackFitPerformance()
else {
PRes2DSec->Fill(mcP.p, 100. * dP);
}
-
- if (std::isfinite(tr.C00[0]) && tr.C00[0] > 0) h_fit[5]->Fill((tr.x[0] - mcP.xIn) / sqrt(tr.C00[0]));
- if (std::isfinite(tr.C11[0]) && tr.C11[0] > 0) h_fit[6]->Fill((tr.y[0] - mcP.yIn) / sqrt(tr.C11[0]));
- if (std::isfinite(tr.C22[0]) && tr.C22[0] > 0) h_fit[7]->Fill((tr.tx[0] - mcP.pxIn / mcP.pzIn) / sqrt(tr.C22[0]));
- if (std::isfinite(tr.C33[0]) && tr.C33[0] > 0) h_fit[8]->Fill((tr.ty[0] - mcP.pyIn / mcP.pzIn) / sqrt(tr.C33[0]));
- if (std::isfinite(tr.C44[0]) && tr.C44[0] > 0) h_fit[9]->Fill((tr.qp[0] - mcP.q / mcP.p) / sqrt(tr.C44[0]));
- h_fit[10]->Fill(tr.qp[0]);
- h_fit[11]->Fill(mcP.q / mcP.p);
- if (isTimeFitted) {
- h_fit[12]->Fill(tr.t[0] - mcP.time);
- if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0.) { h_fit[13]->Fill((tr.t[0] - mcP.time) / sqrt(tr.C55[0])); }
- }
} while (0);
- { // last hit
+ { // last hit
+
int iMC = fvHitPointIndexes[it->Hits.back()]; // TODO2: adapt to linking
if (iMC < 0) continue;
-
-#define L1FSTPARAMEXTRAPOLATE
-#ifdef L1FSTPARAMEXTRAPOLATE
-
- const int last_station = fvHitStore[it->Hits.back()].iStation;
-
- CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
-
- L1FieldRegion fld _fvecalignment;
- fld.SetUseOriginalField();
-
- fit.SetTrack(it->TLast, it->CLast);
- const L1TrackPar& tr = fit.Tr();
-
CbmL1MCPoint& mcP = fvMCPoints[iMC];
- fit.Extrapolate(mcP.zOut, fld);
-
- h_fitL[0]->Fill((tr.x[0] - mcP.xOut) * 1.e4);
- h_fitL[1]->Fill((tr.y[0] - mcP.yOut) * 1.e4);
- h_fitL[2]->Fill((tr.tx[0] - mcP.pxOut / mcP.pzOut) * 1.e3);
- h_fitL[3]->Fill((tr.ty[0] - mcP.pyOut / mcP.pzOut) * 1.e3);
- h_fitL[4]->Fill(100. * (fabs(1. / tr.qp[0]) / mcP.p - 1.));
- if (last_station > fNMvdStations) h_fitL[12]->Fill(tr.t[0] - mcP.time);
-
-
- if (std::isfinite(tr.C00[0]) && tr.C00[0] > 0) h_fitL[5]->Fill((tr.x[0] - mcP.xOut) / sqrt(tr.C00[0]));
- if (std::isfinite(tr.C11[0]) && tr.C11[0] > 0) h_fitL[6]->Fill((tr.y[0] - mcP.yOut) / sqrt(tr.C11[0]));
- if (std::isfinite(tr.C22[0]) && tr.C22[0] > 0)
- h_fitL[7]->Fill((tr.tx[0] - mcP.pxOut / mcP.pzOut) / sqrt(tr.C22[0]));
- if (std::isfinite(tr.C33[0]) && tr.C33[0] > 0)
- h_fitL[8]->Fill((tr.ty[0] - mcP.pyOut / mcP.pzOut) / sqrt(tr.C33[0]));
- if (std::isfinite(tr.C44[0]) && tr.C44[0] > 0) h_fitL[9]->Fill((tr.qp[0] - mcP.q / mcP.p) / sqrt(tr.C44[0]));
- h_fitL[10]->Fill(tr.qp[0]);
- h_fitL[11]->Fill(mcP.q / mcP.p);
- if (last_station > fNMvdStations)
- if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0) h_fitL[13]->Fill((tr.t[0] - mcP.time) / sqrt(tr.C55[0]));
-#else
- CbmL1MCPoint& mc = fvMCPoints[iMC];
-
- h_fitL[0]->Fill((it->TLast[0] - mc.x) * 1.e4);
- h_fitL[1]->Fill((it->TLast[1] - mc.y) * 1.e4);
- h_fitL[2]->Fill((it->TLast[2] - mc.px / mc.pz) * 1.e3);
- h_fitL[3]->Fill((it->TLast[3] - mc.py / mc.pz) * 1.e3);
- h_fitL[4]->Fill(100. * (fabs(1. / it->TLast[4]) / mc.p - 1.));
- if (std::isfinite(it->CLast[0]) && it->CLast[0] > 0) h_fitL[5]->Fill((it->TLast[0] - mc.x) / sqrt(it->CLast[0]));
- if (std::isfinite(it->CLast[2]) && it->CLast[2] > 0) h_fitL[6]->Fill((it->TLast[1] - mc.y) / sqrt(it->CLast[2]));
- if (std::isfinite(it->CLast[5]) && it->CLast[5] > 0)
- h_fitL[7]->Fill((it->TLast[2] - mc.px / mc.pz) / sqrt(it->CLast[5]));
- if (std::isfinite(it->CLast[9]) && it->CLast[9] > 0)
- h_fitL[8]->Fill((it->TLast[3] - mc.py / mc.pz) / sqrt(it->CLast[9]));
- if (std::isfinite(it->CLast[14]) && it->CLast[14] > 0)
- h_fitL[9]->Fill((it->TLast[4] - mc.q / mc.p) / sqrt(it->CLast[14]));
- h_fitL[10]->Fill(it->TLast[4]);
- h_fitL[11]->Fill(mc.q / mc.p);
- h_fitL[12]->Fill(it->TLast[6] - mc.time);
- if (std::isfinite(it->CLast[20]) && it->CLast[20] > 0)
- h_fitL[13]->Fill((it->TLast[6] - mc.time) / sqrt(it->CLast[20]));
-
-#endif
+ L1TrackPar tr;
+ tr.copyFromArrays(it->TLast, it->CLast);
+ FillFitHistos(tr, mcP, isTimeFitted, h_fitL);
}
{ // vertex
CbmL1MCTrack mc = *(it->GetMCTracks()[0]);
fit.SetTrack(it->T, it->C);
+
const L1TrackPar& tr = fit.Tr();
// if (mc.mother_ID != -1) { // secondary
@@ -1414,13 +1374,17 @@ void CbmL1::TrackFitPerformance()
if (std::isfinite(tr.C44[0]) && tr.C44[0] > 0) h_fitSV[9]->Fill((tr.qp[0] - mc.q / mc.p) / sqrt(tr.C44[0]));
h_fitSV[10]->Fill(tr.qp[0]);
h_fitSV[11]->Fill(mc.q / mc.p);
- h_fitSV[12]->Fill(tr.t[0] - mc.time);
- if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0) h_fitSV[13]->Fill((tr.t[0] - mc.time) / sqrt(tr.C55[0]));
+ if (isTimeFitted) {
+ h_fitSV[12]->Fill(tr.t[0] - mc.time);
+ if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0.) { h_fitSV[13]->Fill((tr.t[0] - mc.time) / sqrt(tr.C55[0])); }
+ double dvi = tr.vi[0] - sqrt(1. + mc.mass * mc.mass / mc.p / mc.p) * L1TrackPar::kClightNsInv;
+ h_fitSV[14]->Fill(dvi * L1TrackPar::kClightNs);
+ if (std::isfinite(tr.C66[0]) && tr.C66[0] > 0.) { h_fitSV[15]->Fill(dvi / sqrt(tr.C66[0])); }
+ h_fitSV[16]->Fill(tr.vi[0] * L1TrackPar::kClightNs);
+ }
}
else { // primary
-#define L1EXTRAPOLATE
-#ifdef L1EXTRAPOLATE
{ // extrapolate to vertex
L1FieldRegion fld _fvecalignment;
fld.SetUseOriginalField();
@@ -1479,51 +1443,110 @@ void CbmL1::TrackFitPerformance()
if (std::isfinite(tr.C44[0]) && tr.C44[0] > 0) h_fitPV[9]->Fill((mc.q / mc.p - tr.qp[0]) / sqrt(tr.C44[0]));
h_fitPV[10]->Fill(tr.qp[0]);
h_fitPV[11]->Fill(mc.q / mc.p);
- h_fitPV[12]->Fill(mc.time - tr.t[0]);
- if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0) h_fitPV[13]->Fill((mc.time - tr.t[0]) / sqrt(tr.C55[0]));
-#else
- FairTrackParam fTP;
-
- CbmKFMath::CopyTC2TrackParam(&fTP, it->T, it->C);
+ if (isTimeFitted) {
+ h_fitPV[12]->Fill(tr.t[0] - mc.time);
+ if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0.) { h_fitPV[13]->Fill((tr.t[0] - mc.time) / sqrt(tr.C55[0])); }
+ double dvi = tr.vi[0] - sqrt(1. + mc.mass * mc.mass / mc.p / mc.p) * L1TrackPar::kClightNsInv;
+ h_fitPV[14]->Fill(dvi * L1TrackPar::kClightNs);
+ if (std::isfinite(tr.C66[0]) && tr.C66[0] > 0.) { h_fitPV[15]->Fill(dvi / sqrt(tr.C66[0])); }
+ h_fitPV[16]->Fill(tr.vi[0] * L1TrackPar::kClightNs);
+ }
+ }
+ }
- CbmKFTrack kfTr;
- kfTr.SetTrackParam(fTP);
+ h_fit_chi2->Fill(it->chi2 / it->NDF);
- kfTr.Extrapolate(mc.z);
- CbmL1Track it2;
- for (int ipar = 0; ipar < 6; ipar++)
- it2.T[ipar] = kfTr.GetTrack()[ipar];
- for (int ipar = 0; ipar < 15; ipar++)
- it2.C[ipar] = kfTr.GetCovMatrix()[ipar];
+ // last TRD point
+ do {
+ if (!fpTrdPoints) break;
+ const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]);
+ int nTrdPoints = fpTrdPoints->Size(mcTrack.iFile, mcTrack.iEvent);
+ int lastP = -1;
+ double lastZ = -1.e8;
+ for (int iPoint = 0; iPoint < nTrdPoints; iPoint++) {
+ const CbmTrdPoint* pt = dynamic_cast<CbmTrdPoint*>(fpTrdPoints->Get(mcTrack.iFile, mcTrack.iEvent, iPoint));
+ if (!pt) { continue; }
+ if (pt->GetTrackID() != mcTrack.ID) { continue; }
+ if (lastP < 0 || lastZ < pt->GetZ()) {
+ lastP = iPoint;
+ lastZ = pt->GetZ();
+ }
+ }
+ CbmL1MCPoint mcP;
+ if (CbmL1::ReadMCPoint(&mcP, lastP, mcTrack.iFile, mcTrack.iEvent, 3)) { break; }
+ L1TrackPar tr;
+ tr.copyFromArrays(it->TLast, it->CLast);
+ FillFitHistos(tr, mcP, isTimeFitted, h_fitTrd);
+ } while (0); // Trd point
- // calculate pulls
- // h_fitPV[0]->Fill( (mc.x-it2.T[0]) *1.e4);
- // h_fitPV[1]->Fill( (mc.y-it2.T[1]) *1.e4);
- h_fitPV[0]->Fill((mc.x - it2.T[0]));
- h_fitPV[1]->Fill((mc.y - it2.T[1]));
- h_fitPV[2]->Fill((mc.px / mc.pz - it2.T[2]) * 1.e3);
- h_fitPV[3]->Fill((mc.py / mc.pz - it2.T[3]) * 1.e3);
- h_fitPV[4]->Fill(100. * (it2.T[4] / mc.q * mc.p - 1.));
- if (std::isfinite(it2.C[0]) && it2.C[0] > 0) h_fitPV[5]->Fill((mc.x - it2.T[0]) / sqrt(it2.C[0]));
- if (std::isfinite(it2.C[2]) && it2.C[2] > 0) h_fitPV[6]->Fill((mc.y - it2.T[1]) / sqrt(it2.C[2]));
- if (std::isfinite(it2.C[5]) && it2.C[5] > 0) h_fitPV[7]->Fill((mc.px / mc.pz - it2.T[2]) / sqrt(it2.C[5]));
- if (std::isfinite(it2.C[9]) && it2.C[9] > 0) h_fitPV[8]->Fill((mc.py / mc.pz - it2.T[3]) / sqrt(it2.C[9]));
- if (std::isfinite(it2.C[14]) && it2.C[14] > 0) h_fitPV[9]->Fill((mc.q / mc.p - it2.T[4]) / sqrt(it2.C[14]));
- h_fitPV[10]->Fill(it2.T[4]);
- h_fitPV[11]->Fill(mc.q / mc.p);
- h_fitPV[12]->Fill(mc.time - it2.T[6]);
- if (std::isfinite(it2.C[20]) && it2.C[20] > 0) h_fitPV[13]->Fill((mc.time - it2.T[6]) / sqrt(it2.C[20]));
-#endif // L1EXTRAPOLATE
+ // last TOF point
+ do {
+ if (!fpTofPoints) break;
+ const CbmL1MCTrack& mcTrack = *(it->GetMCTracks()[0]);
+ int nTofPoints = fpTofPoints->Size(mcTrack.iFile, mcTrack.iEvent);
+ int lastP = -1;
+ double lastZ = -1.e8;
+ for (int iPoint = 0; iPoint < nTofPoints; iPoint++) {
+ const CbmTofPoint* pt = dynamic_cast<CbmTofPoint*>(fpTofPoints->Get(mcTrack.iFile, mcTrack.iEvent, iPoint));
+ if (!pt) { continue; }
+ if (pt->GetTrackID() != mcTrack.ID) { continue; }
+ if (lastP < 0 || lastZ < pt->GetZ()) {
+ lastP = iPoint;
+ lastZ = pt->GetZ();
+ }
}
- }
+ CbmL1MCPoint mcP;
+ if (CbmL1::ReadMCPoint(&mcP, lastP, mcTrack.iFile, mcTrack.iEvent, 4)) { break; }
+ L1TrackPar tr;
+ tr.copyFromArrays(it->TLast, it->CLast);
+ FillFitHistos(tr, mcP, isTimeFitted, h_fitTof);
+ } while (0); // tof point
- h_fit_chi2->Fill(it->chi2 / it->NDF);
- }
+ } // reco track
} // void CbmL1::TrackFitPerformance()
+void CbmL1::FillFitHistos(L1TrackPar& track, const CbmL1MCPoint& mcP, bool isTimeFitted, TH1F* h[])
+{
+ L1Fit fit;
+ fit.SetParticleMass(fpAlgo->GetDefaultParticleMass());
+ fit.SetMask(fmask::One());
+ //fit.SetMaxExtrapolationStep(10.);
+ fit.SetDoFitVelocity(true);
+ fit.SetTrack(track);
+ L1FieldRegion fld _fvecalignment;
+ fld.SetUseOriginalField();
+ fit.Extrapolate(mcP.zOut, fld);
+ track = fit.Tr();
+
+ const L1TrackPar& tr = track;
+
+ h[0]->Fill((tr.x[0] - mcP.xOut) * 1.e4);
+ h[1]->Fill((tr.y[0] - mcP.yOut) * 1.e4);
+ h[2]->Fill((tr.tx[0] - mcP.pxOut / mcP.pzOut) * 1.e3);
+ h[3]->Fill((tr.ty[0] - mcP.pyOut / mcP.pzOut) * 1.e3);
+ h[4]->Fill(100. * (fabs(1. / tr.qp[0]) / mcP.p - 1.));
+
+ if (std::isfinite(tr.C00[0]) && tr.C00[0] > 0) h[5]->Fill((tr.x[0] - mcP.xOut) / sqrt(tr.C00[0]));
+ if (std::isfinite(tr.C11[0]) && tr.C11[0] > 0) h[6]->Fill((tr.y[0] - mcP.yOut) / sqrt(tr.C11[0]));
+ if (std::isfinite(tr.C22[0]) && tr.C22[0] > 0) h[7]->Fill((tr.tx[0] - mcP.pxOut / mcP.pzOut) / sqrt(tr.C22[0]));
+ if (std::isfinite(tr.C33[0]) && tr.C33[0] > 0) h[8]->Fill((tr.ty[0] - mcP.pyOut / mcP.pzOut) / sqrt(tr.C33[0]));
+ if (std::isfinite(tr.C44[0]) && tr.C44[0] > 0) h[9]->Fill((tr.qp[0] - mcP.q / mcP.p) / sqrt(tr.C44[0]));
+ h[10]->Fill(tr.qp[0]);
+ h[11]->Fill(mcP.q / mcP.p);
+ if (isTimeFitted) {
+ h[12]->Fill(tr.t[0] - mcP.time);
+ if (std::isfinite(tr.C55[0]) && tr.C55[0] > 0.) { h[13]->Fill((tr.t[0] - mcP.time) / sqrt(tr.C55[0])); }
+ double dvi = tr.vi[0] - sqrt(1. + mcP.mass * mcP.mass / mcP.p / mcP.p) * L1TrackPar::kClightNsInv;
+ h[14]->Fill(dvi * L1TrackPar::kClightNs);
+ if (std::isfinite(tr.C66[0]) && tr.C66[0] > 0.) { h[15]->Fill(dvi / sqrt(tr.C66[0])); }
+ h[16]->Fill(tr.vi[0] * L1TrackPar::kClightNs);
+ }
+}
+
+
void CbmL1::FieldApproxCheck()
{
TDirectory* curr = gDirectory;
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 2dbc50b3089d2cd614eb939061890a38562ccc13..22a86e698267d1e03ae30a22cd016a2628bde890 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -88,7 +88,7 @@ struct TmpHit {
double ty; ///< MC slopes of the mc point
int iMC; ///< index of MCPoint in the fvMCPoints array
double time = 0.; ///< time of the hit [ns]
- double dt = 1.e10; ///< time error of the hit [ns]
+ double dt = 1.e4; ///< time error of the hit [ns]
int Det;
int id; ///< index of hit before hits sorting
int track;
@@ -1358,6 +1358,8 @@ void CbmL1::Fill_vMCTracks()
auto header = dynamic_cast<FairMCEventHeader*>(fpMcEventHeader->Get(iFile, iEvent));
assert(header);
+ double eventTime = fpEventList->GetEventTime(iEvent, iFile);
+
Int_t nMCTrack = fpMCTracks->Size(iFile, iEvent);
/* Loop over MC tracks */
for (Int_t iMCTrack = 0; iMCTrack < nMCTrack; iMCTrack++) {
@@ -1385,7 +1387,7 @@ void CbmL1::Fill_vMCTracks()
Int_t iTrack = fvMCTracks.size(); //or iMCTrack?
CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg, processID);
// CbmL1MCTrack T(mass, q, vr, vp, iMCTrack, mother_ID, pdg);
- T.time = MCTrack->GetStartT();
+ T.time = eventTime + MCTrack->GetStartT();
T.iFile = iFile;
T.iEvent = iEvent;
// signal: primary tracks, displaced from the primary vertex
@@ -1474,9 +1476,9 @@ bool CbmL1::ReadMCPoint(CbmL1MCPoint* MC, int iPoint, int file, int event, int M
if (!pt) return 1;
if (!fLegacyEventMode) {
- Double_t StartTime = fTimeSlice->GetStartTime(); // not used
- Double_t EndTime = fTimeSlice->GetEndTime(); // not used
- Double_t Time_MC_point = eventTime + pt->GetTime(); // not used
+ Double_t StartTime = fTimeSlice->GetStartTime(); // not used
+ Double_t EndTime = fTimeSlice->GetEndTime(); // not used
+ Double_t Time_MC_point = eventTime + pt->GetTime(); // not used
if (StartTime > 0)
if (Time_MC_point < StartTime) return 1;
diff --git a/reco/L1/CbmL1Track.h b/reco/L1/CbmL1Track.h
index 4286a905adb66e20960ef68cbc438764cf8a63d2..e5275831aa65903f7d64253ff4ca9231959db4ea 100644
--- a/reco/L1/CbmL1Track.h
+++ b/reco/L1/CbmL1Track.h
@@ -51,10 +51,10 @@ public:
static bool comparePChi2(const CbmL1Track* a, const CbmL1Track* b) { return (a->chi2 < b->chi2); }
- double Tpv[7], Cpv[21];
+ double Tpv[L1TrackPar::kNparTr], Cpv[L1TrackPar::kNparCov];
- double TLast[7], CLast[21];
+ double TLast[L1TrackPar::kNparTr], CLast[L1TrackPar::kNparCov];
vector<int> Hits;
int nStations;
int index;
diff --git a/reco/L1/CbmL1TrackPar.h b/reco/L1/CbmL1TrackPar.h
index d1d69306bc87d3fd9f9574c9cd59344718f34449..57fe79607bc74bd4bb512937114484cd34cbeecf 100644
--- a/reco/L1/CbmL1TrackPar.h
+++ b/reco/L1/CbmL1TrackPar.h
@@ -1,12 +1,14 @@
/* Copyright (C) 2006-2017 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
SPDX-License-Identifier: GPL-3.0-only
- Authors: Denis Bertini [committer] */
+ Authors: Sergey Gorbunov [committer] */
#ifndef CbmL1TrackPar_H
#define CbmL1TrackPar_H
#include "CbmKFTrackInterface.h"
+#include "L1TrackPar.h"
+
struct CbmL1TrackPar : public CbmKFTrackInterface {
public:
CbmL1TrackPar() : chi2(0), NDF(0), mass(0), is_electron(0) {}
@@ -17,7 +19,7 @@ public:
double GetMass() { return mass; }
bool IsElectron() { return is_electron; }
- double T[7], C[21], chi2;
+ double T[L1TrackPar::kNparTr], C[L1TrackPar::kNparCov], chi2;
int NDF;
double mass; // mass hypothesis
bool is_electron;
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 00d60be0de2bf0ec1f6a6dd3cd1e6de48d1348fd..9aad594f2904d6aecad6c552b8db88cf604550aa 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -70,8 +70,8 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool di
fvec dzi = fvec(1.) / (z1 - z0);
- T.x = x0;
- T.y = y0;
+ T.x = x0;
+ T.y = y0;
if (initParams) {
T.tx = (x1 - x0) * dzi;
T.ty = (y1 - y0) * dzi;
@@ -242,7 +242,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool dir,
const L1Hit& hit = (*vHitsUnused)[globalInd];
- if (sta.timeInfo) {
+ if (sta.timeInfo && tr.nTimeMeasurements[0] > 0) {
fscal dt = hit.t - tr.t[0];
if (dt * dt > (tr.C55[0] + hit.dt2) * 25) continue;
}
diff --git a/reco/L1/L1Algo/L1CloneMerger.cxx b/reco/L1/L1Algo/L1CloneMerger.cxx
index 16a469c6d2615c02c57f08c9bfe70a096d718689..cc5fd275e9745dd0465ffccb2809cc0d76230688 100644
--- a/reco/L1/L1Algo/L1CloneMerger.cxx
+++ b/reco/L1/L1Algo/L1CloneMerger.cxx
@@ -114,7 +114,10 @@ void L1CloneMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>& e
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;
+ if (Tf.nTimeMeasurements[0] > 0 && Tb.nTimeMeasurements[0] > 0) {
+ if (fabs(Tf.t[0] - Tb.t[0]) > 3 * sqrt(Tf.C55[0] + Tb.C55[0])) continue;
+ }
+
unsigned short stam;
frAlgo.GetParameters()->GetStation(staf).fieldSlice.GetFieldValue(Tf.x, Tf.y, fBf);
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index 3a99d903dbe6b6dd83d5cf418b5f1a23a9eaa174..8fe46a0534673e32b0b472161869f44ca4546da6 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -55,26 +55,32 @@ void L1Fit::Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2)
fTr.vi -= F6 * zetawi;
fTr.C00 -= F0 * F0 * wi;
+
fTr.C10 -= K1 * F0;
fTr.C11 -= K1 * F1;
+
fTr.C20 -= K2 * F0;
fTr.C21 -= K2 * F1;
fTr.C22 -= K2 * F2;
+
fTr.C30 -= K3 * F0;
fTr.C31 -= K3 * F1;
fTr.C32 -= K3 * F2;
fTr.C33 -= K3 * F3;
+
fTr.C40 -= K4 * F0;
fTr.C41 -= K4 * F1;
fTr.C42 -= K4 * F2;
fTr.C43 -= K4 * F3;
fTr.C44 -= K4 * F4;
+
fTr.C50 -= K5 * F0;
fTr.C51 -= K5 * F1;
fTr.C52 -= K5 * F2;
fTr.C53 -= K5 * F3;
fTr.C54 -= K5 * F4;
fTr.C55 -= K5 * F5;
+
fTr.C60 -= K6 * F0;
fTr.C61 -= K6 * F1;
fTr.C62 -= K6 * F2;
@@ -133,27 +139,35 @@ void L1Fit::FilterTime(cnst& t, cnst& dt2, cnst& timeInfo)
fTr.t -= F5 * zetawi;
fTr.vi -= F6 * zetawi;
+ fTr.nTimeMeasurements += w;
+
fTr.C00 -= F0 * F0 * wi;
+
fTr.C10 -= K1 * F0;
fTr.C11 -= K1 * F1;
+
fTr.C20 -= K2 * F0;
fTr.C21 -= K2 * F1;
fTr.C22 -= K2 * F2;
+
fTr.C30 -= K3 * F0;
fTr.C31 -= K3 * F1;
fTr.C32 -= K3 * F2;
fTr.C33 -= K3 * F3;
+
fTr.C40 -= K4 * F0;
fTr.C41 -= K4 * F1;
fTr.C42 -= K4 * F2;
fTr.C43 -= K4 * F3;
fTr.C44 -= K4 * F4;
+
fTr.C50 -= K5 * F0;
fTr.C51 -= K5 * F1;
fTr.C52 -= K5 * F2;
fTr.C53 -= K5 * F3;
fTr.C54 -= K5 * F4;
fTr.C55 -= K5 * F5;
+
fTr.C60 -= K6 * F0;
fTr.C61 -= K6 * F1;
fTr.C62 -= K6 * F2;
@@ -209,16 +223,22 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
K00 = F00 * S00 + F01 * S10;
K01 = F00 * S10 + F01 * S11;
+
K10 = F10 * S00 + F11 * S10;
K11 = F10 * S10 + F11 * S11;
+
K20 = F20 * S00 + F21 * S10;
K21 = F20 * S10 + F21 * S11;
+
K30 = F30 * S00 + F31 * S10;
K31 = F30 * S10 + F31 * S11;
+
K40 = F40 * S00 + F41 * S10;
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;
@@ -347,6 +367,188 @@ void L1Fit::FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& in
}
+void L1Fit::MeasureVelocityWithQp()
+{
+ // measure velocity using measured qp
+ // assuming particle mass == fMass;
+
+ cnst kClightNsInv = fvec(L1TrackPar::kClightNsInv);
+
+ fvec zeta, HCH;
+ fvec F0, F1, F2, F3, F4, F5, F6;
+ fvec K1, K2, K3, K4, K5, K6;
+
+ //FilterVi(sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv);
+ //return;
+
+ fvec vi0 = sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
+
+ fvec h = fMass2 * fQp0 / sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * kClightNsInv;
+
+ zeta = vi0 + h * (fTr.qp - fQp0) - fTr.vi;
+
+ fTr.vi = vi0;
+
+ // H = (0,0,0,0, h,0, -1)
+
+ // F = CH'
+
+ F0 = h * fTr.C40 - fTr.C60;
+ F1 = h * fTr.C41 - fTr.C61;
+ F2 = h * fTr.C42 - fTr.C62;
+ F3 = h * fTr.C43 - fTr.C63;
+ F4 = h * fTr.C44 - fTr.C64;
+ F5 = h * fTr.C54 - fTr.C65;
+ F6 = h * fTr.C64 - fTr.C66;
+
+ HCH = F4 * h - F6;
+
+ fvec wi = fMaskF / HCH;
+ fvec zetawi = fMaskF * zeta / HCH;
+ fTr.chi2 += zeta * zeta * wi;
+ fTr.NDF += fMaskF;
+
+ K1 = F1 * wi;
+ K2 = F2 * wi;
+ K3 = F3 * wi;
+ K4 = F4 * wi;
+ K5 = F5 * wi;
+ K6 = F6 * wi;
+
+ fTr.x -= F0 * zetawi;
+ fTr.y -= F1 * zetawi;
+ fTr.tx -= F2 * zetawi;
+ fTr.ty -= F3 * zetawi;
+ fTr.qp -= F4 * zetawi;
+ fTr.t -= F5 * zetawi;
+ fTr.vi -= F6 * zetawi;
+
+ fTr.C00 -= F0 * F0 * wi;
+
+ fTr.C10 -= K1 * F0;
+ fTr.C11 -= K1 * F1;
+
+ fTr.C20 -= K2 * F0;
+ fTr.C21 -= K2 * F1;
+ fTr.C22 -= K2 * F2;
+
+ fTr.C30 -= K3 * F0;
+ fTr.C31 -= K3 * F1;
+ fTr.C32 -= K3 * F2;
+ fTr.C33 -= K3 * F3;
+
+ fTr.C40 -= K4 * F0;
+ fTr.C41 -= K4 * F1;
+ fTr.C42 -= K4 * F2;
+ fTr.C43 -= K4 * F3;
+ fTr.C44 -= K4 * F4;
+
+ fTr.C50 -= K5 * F0;
+ fTr.C51 -= K5 * F1;
+ fTr.C52 -= K5 * F2;
+ fTr.C53 -= K5 * F3;
+ fTr.C54 -= K5 * F4;
+ fTr.C55 -= K5 * F5;
+
+ 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;
+
+ // fTr.vi( fTr.vi < fvec(L1TrackPar::kClightNsInv) ) = fvec(L1TrackPar::kClightNsInv);
+}
+
+void L1Fit::FilterVi(fvec vi)
+{
+ // set inverse velocity to vi
+
+ fvec zeta, HCH;
+ fvec F0, F1, F2, F3, F4, F5, F6;
+ fvec K1, K2, K3, K4, K5, K6;
+
+ zeta = fTr.vi - vi;
+
+ // H = (0,0,0,0, 0, 0, 1)
+
+ // F = CH'
+
+ F0 = fTr.C60;
+ F1 = fTr.C61;
+ F2 = fTr.C62;
+ F3 = fTr.C63;
+ F4 = fTr.C64;
+ F5 = fTr.C65;
+ F6 = fTr.C66;
+
+ HCH = F6;
+
+ fvec wi = fMaskF / HCH;
+ fvec zetawi = fMaskF * zeta / HCH;
+ fTr.chi2 += zeta * zeta * wi;
+ fTr.NDF += fMaskF;
+
+ K1 = F1 * wi;
+ K2 = F2 * wi;
+ K3 = F3 * wi;
+ K4 = F4 * wi;
+ K5 = F5 * wi;
+ K6 = F6 * wi;
+
+ fTr.x -= F0 * zetawi;
+ fTr.y -= F1 * zetawi;
+ fTr.tx -= F2 * zetawi;
+ fTr.ty -= F3 * zetawi;
+ fTr.qp -= F4 * zetawi;
+ fTr.t -= F5 * zetawi;
+ // fTr.vi -= F6 * zetawi;
+ fTr.vi = vi;
+
+ fTr.C00 -= F0 * F0 * wi;
+
+ fTr.C10 -= K1 * F0;
+ fTr.C11 -= K1 * F1;
+
+ fTr.C20 -= K2 * F0;
+ fTr.C21 -= K2 * F1;
+ fTr.C22 -= K2 * F2;
+
+ fTr.C30 -= K3 * F0;
+ fTr.C31 -= K3 * F1;
+ fTr.C32 -= K3 * F2;
+ fTr.C33 -= K3 * F3;
+
+ fTr.C40 -= K4 * F0;
+ fTr.C41 -= K4 * F1;
+ fTr.C42 -= K4 * F2;
+ fTr.C43 -= K4 * F3;
+ fTr.C44 -= K4 * F4;
+
+ fTr.C50 -= K5 * F0;
+ fTr.C51 -= K5 * F1;
+ fTr.C52 -= K5 * F2;
+ fTr.C53 -= K5 * F3;
+ fTr.C54 -= K5 * F4;
+ fTr.C55 -= K5 * F5;
+
+ //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;
+ fTr.C60 = fvec(0.);
+ fTr.C61 = fvec(0.);
+ fTr.C62 = fvec(0.);
+ fTr.C63 = fvec(0.);
+ fTr.C64 = fvec(0.);
+ fTr.C65 = fvec(0.);
+ fTr.C66 = fvec(1.e-8); // just for a case..
+}
+
void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& z_out, // extrapolate to this z position
const L1FieldRegion& F)
@@ -355,7 +557,7 @@ void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian f
fvec sgn = iif(fTr.z < z_out, fvec(1.), fvec(-1.));
while (!(fMaskF * abs(z_out - fTr.z) <= fvec(1.e-6)).isFull()) {
- fvec zNew = fTr.z + sgn * fvec(50.); // max. 50 cm step
+ fvec zNew = fTr.z + sgn * fMaxExtraplationStep; // max. 50 cm step
zNew(sgn * (z_out - zNew) <= fvec(0.)) = z_out;
ExtrapolateStepFull(zNew, F);
}
@@ -363,7 +565,7 @@ void L1Fit::Extrapolate // extrapolates track parameters and returns jacobian f
void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
(cnst& zOut, // extrapolate to this z position
- const L1FieldRegion& F)
+ const L1FieldRegion& Field)
{
// use Q/p linearisation at fQp0
// implementation of the Runge-Kutta method without optimization
@@ -406,235 +608,130 @@ void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns ja
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.)};
+ cnst h = (zMasked - fTr.z);
- 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]
+ cnst stepDz[5] = {0., 0., h * fvec(0.5), h * fvec(0.5), h};
- //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 f[5][7] = {{0.}}; // ( d*/dz ) [step]
+ fvec F[5][7][7] = {{0.}}; // ( d *new [step] / d *old )
- fvec k[5][7] = {{0.}}; // [ step-1 ] [ {x,y,tx,ty,qp,t,vi} ]
-
- // Runge-Kutta steps for track parameters
+ // Runge-Kutta steps
//
- {
- 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;
+ fvec r0[7] = {fTr.x, fTr.y, fTr.tx, fTr.ty, fQp0, fTr.t, fTr.vi};
+ fvec R0[7][7] = {{0.}};
+ for (int i = 0; i < 7; ++i) {
+ R0[i][i] = 1.;
}
+ for (int step = 1; step <= 4; ++step) {
- // 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.;
- }
+ fvec rstep[7] = {{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];
+ rstep[i] = r0[i] + stepDz[step] * f[step - 1][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.;
+ fvec z = fTr.z + stepDz[step];
+
+ fvec B[3];
+ cnst& Bx = B[0];
+ cnst& By = B[1];
+ cnst& Bz = B[2];
+ Field.Get(rstep[0], rstep[1], z, B);
+
+ fvec tx = rstep[2];
+ fvec ty = rstep[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;
+
+ f[step][0] = tx;
+ F[step][0][2] = 1.;
+
+ f[step][1] = ty;
+ F[step][1][3] = 1.;
+
+ fvec f2tmp = txty * Bx - (fvec(1.) + tx2) * By + ty * Bz;
+ f[step][2] = cLqp0 * f2tmp;
+
+ F[step][2][2] = cLqp0 * (tx * f2tmp * L2i + ty * Bx - fvec(2.) * tx * By);
+ F[step][2][3] = cLqp0 * (ty * f2tmp * L2i + tx * Bx + Bz);
+ F[step][2][4] = cL * f2tmp;
+
+ fvec f3tmp = -txty * By - tx * Bz + (fvec(1.) + ty2) * Bx;
+ f[step][3] = cLqp0 * f3tmp;
+ F[step][3][2] = cLqp0 * (tx * f3tmp * L2i - ty * By - Bz);
+ F[step][3][3] = cLqp0 * (ty * f3tmp * L2i + fvec(2.) * ty * Bx - tx * By);
+ F[step][3][4] = cL * f3tmp;
+
+ f[step][4] = 0.;
+
+ if (fDoFitVelocity) {
+ fvec vi = rstep[6];
+ f[step][5] = vi * L;
+ F[step][5][2] = vi * tx / L;
+ F[step][5][3] = vi * ty / L;
+ F[step][5][4] = 0.;
+ F[step][5][5] = 0.;
+ F[step][5][6] = L;
}
- 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];
+ else {
+ fvec vi = sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * fvec(L1TrackPar::kClightNsInv);
+ f[step][5] = vi * L;
+ F[step][5][2] = vi * tx / L;
+ F[step][5][3] = vi * ty / L;
+ F[step][5][4] = fMass2 * fQp0 * L / sqrt(fvec(1.) + fMass2 * fQp0 * fQp0) * fvec(L1TrackPar::kClightNsInv);
+ F[step][5][5] = 0.;
+ F[step][5][6] = 0.;
}
- }
- //
- // 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];
- }
- }
+ f[step][6] = 0.;
+ } // end of Runge-Kutta step
- //
- // 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
+ fvec r[7] = {{0.}}; // extrapolated parameters
+ fvec R[7][7] = {{0.}}; // Jacobian of the extrapolation
- //
- // 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];
+ cnst stepW[5] = {0., h / fvec(6.), h / fvec(3.), h / fvec(3.), h / fvec(6.)};
+
+ for (int i = 0; i < 7; i++) {
+ for (int j = 0; j < 7; j++) {
+ R[i][j] = R0[i][j];
+ for (int step = 1; step <= 4; step++) {
+ R[i][j] += stepW[step] * F[step][i][j];
+ }
}
}
+ fvec dqp = fTr.qp - fQp0;
- { // 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;
+ for (int i = 0; i < 7; i++) {
+ r[i] = r0[i];
+ for (int step = 1; step <= 4; step++) {
+ r[i] += stepW[step] * f[step][i];
+ }
+ // take into account linearisation at fQp0
+ r[i] += R[i][4] * dqp;
}
- // covariance matrix transport
+ // update parameters
- 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
- }
+ fTr.x = r[0];
+ fTr.y = r[1];
+ fTr.tx = r[2];
+ fTr.ty = r[3];
+ fTr.qp = r[4];
+ fTr.t = r[5];
+ fTr.vi = r[6];
+
+ //fTr.vi( fTr.vi < fvec(L1TrackPar::kClightNsInv) ) = fvec(L1TrackPar::kClightNsInv);
+ fTr.z = zMasked;
+
+ // covariance matrix transport
fvec C[7][7];
for (int i = 0; i < 7; i++) {
@@ -643,12 +740,12 @@ void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns ja
}
}
- fvec JC[7][7];
+ fvec RC[7][7];
for (int i = 0; i < 7; i++) {
for (int j = 0; j < 7; j++) {
- JC[i][j] = 0.;
+ RC[i][j] = 0.;
for (int m = 0; m < 7; m++) {
- JC[i][j] += J[i][m] * C[m][j];
+ RC[i][j] += R[i][m] * C[m][j];
}
}
}
@@ -656,7 +753,7 @@ void L1Fit::ExtrapolateStepFull // extrapolates track parameters and returns ja
for (int j = 0; j < 7; j++) {
fvec Cij = 0.;
for (int m = 0; m < 7; m++) {
- Cij += JC[i][m] * J[j][m];
+ Cij += RC[i][m] * R[j][m];
}
fTr.C(i, j) = Cij;
}
@@ -1849,7 +1946,7 @@ void L1Fit::GuessTrack(const fvec& trackZ, const fvec hitX[], const fvec hitY[],
fTr.qp = -L * c_light_i / sqrt(txtx1 + fTr.ty * fTr.ty);
fTr.t = time;
fTr.z = trackZ;
- fTr.vi = 0.;
+ fTr.vi = L1TrackPar::kClightNsInv;
fQp0 = fTr.qp;
}
diff --git a/reco/L1/L1Algo/L1Fit.h b/reco/L1/L1Algo/L1Fit.h
index 6ca7a9b1f93f87f1ee1773e0072107db01b3297f..458d9d576ac1989efc9c06b5987f69bc4a0cd182 100644
--- a/reco/L1/L1Algo/L1Fit.h
+++ b/reco/L1/L1Algo/L1Fit.h
@@ -83,11 +83,21 @@ public:
/// get the particle mass squared
fvec GetParticleMass2() const { return fMass2; }
+ /// set max extrapolation step [cm]
+ void SetMaxExtrapolationStep(fscal step) { fMaxExtraplationStep = fvec(step); }
+
+ /// get the particle mass
+ fvec GetMaxExtrapolationStep() const { return fMaxExtraplationStep; }
+
+
void Filter(const L1UMeasurementInfo& info, cnst& u, cnst& sigma2);
void FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y);
void FilterTime(cnst& t, cnst& dt2, cnst& timeInfo);
void FilterExtrapolatedXY(cnst& x, cnst& y, const L1XYMeasurementInfo& info, cnst& extrX, cnst& extrY, cnst Jx[6],
cnst Jy[6]);
+ void FilterVi(fvec vi);
+
+ void MeasureVelocityWithQp();
void Extrapolate(cnst& z_out, const L1FieldRegion& F);
@@ -143,6 +153,8 @@ private:
fvec fMass {0.10565800}; // muon mass
fvec fMass2 {fMass * fMass}; // mass squared
+ fvec fMaxExtraplationStep {50.}; // max extrapolation step [cm]
+
fvec fPipeRadThick {7.87e-3f}; // 0.7 mm Aluminium // TODO:
fvec fTargetRadThick {3.73e-2f * 2}; // 250 mum Gold // TODO:
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 4676463254b70d2e098154f8162a68efa39bf1db..6f8330d0ea5b79901056e0a373bfd8d6dee703a2 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -38,11 +38,10 @@ void L1Algo::L1KFTrackFitter()
const int nStations = fParameters.GetNstationsActive();
int nTracks_SIMD = fvec::size();
- L1Fit fit; // fitting parametr coresponding to current track
+ L1Fit fit; // fit parameters coresponding to the current track
L1TrackPar& tr = fit.Tr();
-
fit.SetParticleMass(GetDefaultParticleMass());
- fit.SetDoFitVelocity(false);
+ fit.SetDoFitVelocity(true);
L1Track* t[fvec::size()] {nullptr};
@@ -73,6 +72,7 @@ void L1Algo::L1KFTrackFitter()
fvec d_xy_fst;
fvec time_first;
+ fvec wtime_first;
fvec dt2_first;
fvec x_last;
@@ -82,6 +82,7 @@ void L1Algo::L1KFTrackFitter()
fvec d_xy_lst;
fvec time_last;
+ fvec wtime_last;
fvec dt2_last;
// fvec dt2_lst; /// TODO: Why are there two different variables for the time error on the last station?
@@ -154,7 +155,8 @@ void L1Algo::L1KFTrackFitter()
y[ista][iVec] = y_temp[iVec];
time[ista][iVec] = hit.t;
dt2[ista][iVec] = hit.dt2;
- z[ista][iVec] = hit.z;
+ if (!sta[ista].timeInfo) { dt2[ista][iVec] = 1.e4; }
+ z[ista][iVec] = hit.z;
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]);
@@ -162,25 +164,27 @@ void L1Algo::L1KFTrackFitter()
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
if (ih == 0) {
- z_start[iVec] = z[ista][iVec];
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
- time_first[iVec] = time[ista][iVec];
- dt2_first[iVec] = dt2[ista][iVec];
- d_xx_fst[iVec] = d_xx[ista][iVec];
- d_yy_fst[iVec] = d_yy[ista][iVec];
- d_xy_fst[iVec] = d_xy[ista][iVec];
+ z_start[iVec] = z[ista][iVec];
+ x_first[iVec] = x[ista][iVec];
+ y_first[iVec] = y[ista][iVec];
+ time_first[iVec] = time[ista][iVec];
+ wtime_first[iVec] = sta[ista].timeInfo ? 1. : 0.;
+ dt2_first[iVec] = dt2[ista][iVec];
+ d_xx_fst[iVec] = d_xx[ista][iVec];
+ d_yy_fst[iVec] = d_yy[ista][iVec];
+ d_xy_fst[iVec] = d_xy[ista][iVec];
}
else if (ih == nHitsTrack - 1) {
- z_end[iVec] = z[ista][iVec];
- x_last[iVec] = x[ista][iVec];
- y_last[iVec] = y[ista][iVec];
- d_xx_lst[iVec] = d_xx[ista][iVec];
- d_yy_lst[iVec] = d_yy[ista][iVec];
- d_xy_lst[iVec] = d_xy[ista][iVec];
- time_last[iVec] = time[ista][iVec];
- dt2_last[iVec] = dt2[ista][iVec];
+ z_end[iVec] = z[ista][iVec];
+ x_last[iVec] = x[ista][iVec];
+ y_last[iVec] = y[ista][iVec];
+ d_xx_lst[iVec] = d_xx[ista][iVec];
+ d_yy_lst[iVec] = d_yy[ista][iVec];
+ d_xy_lst[iVec] = d_xy[ista][iVec];
+ time_last[iVec] = time[ista][iVec];
+ dt2_last[iVec] = dt2[ista][iVec];
+ wtime_last[iVec] = sta[ista].timeInfo ? 1. : 0.;
}
}
@@ -208,13 +212,17 @@ void L1Algo::L1KFTrackFitter()
time_last = iif(w_time[ista], time_last, fvec::Zero());
dt2_last = iif(w_time[ista], dt2_last, fvec(1.e6));
- tr.ResetErrors(d_xx_lst, d_yy_lst, 0.1, 0.1, 1.0, dt2_last, 1.e2);
+ tr.ResetErrors(d_xx_lst, d_yy_lst, 0.1, 0.1, 1.0, dt2_last, 1.e-2);
tr.C10 = d_xy_lst;
tr.x = x_last;
tr.y = y_last;
tr.t = time_last;
+ tr.vi = L1TrackPar::kClightNsInv;
+ tr.InitVelocityRange(0.5);
tr.NDF = -3.0;
+ tr.nTimeMeasurements = wtime_last;
+
fldZ1 = z[ista];
sta[ista].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1);
@@ -286,6 +294,11 @@ void L1Algo::L1KFTrackFitter()
}
}
+ //fit.SetQp0(tr.qp);
+ //fit.SetMask(fmask::One());
+ //fit.MeasureVelocityWithQp();
+ //fit.FilterVi(L1TrackPar::kClightNsInv);
+
L1TrackPar Tf = fit.Tr();
if (kGlobal == fTrackingMode) { Tf.qp = fitpv.Tr().qp; }
@@ -299,20 +312,25 @@ void L1Algo::L1KFTrackFitter()
fitpv.Tr().copyToArrays(iVec, t[iVec]->Tpv, t[iVec]->Cpv);
}
- if (iter == 1) continue; // only 1.5 iterations
+ if (iter == 1) { break; } // only 1.5 iterations
// fit forward
ista = 0;
- tr.ResetErrors(d_xx_fst, d_yy_fst, 0.1, 0.1, 1., dt2_first, 1.e2);
+ tr.ResetErrors(d_xx_fst, d_yy_fst, 0.1, 0.1, 1., dt2_first, 1.e-2);
tr.C10 = d_xy_fst;
- tr.x = x_first;
- tr.y = y_first;
- tr.t = time_first;
+ tr.x = x_first;
+ tr.y = y_first;
+ tr.t = time_first;
+ tr.vi = L1TrackPar::kClightNsInv;
+ tr.InitVelocityRange(0.5);
+
tr.NDF = -3.0;
+ tr.nTimeMeasurements = wtime_first;
+
fit.SetQp0(tr.qp);
fldZ1 = z[ista];
@@ -349,6 +367,11 @@ void L1Algo::L1KFTrackFitter()
fldZ1 = fldZ0;
}
+ //fit.SetQp0(tr.qp);
+ //fit.SetMask(fmask::One());
+ //fit.MeasureVelocityWithQp();
+ //fit.FilterVi(L1TrackPar::kClightNsInv);
+
L1TrackPar Tl = fit.Tr();
if (kGlobal == fTrackingMode) { Tl.qp = fitpv.Tr().qp; }
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 25e4fb43a7d60be9a88519c5683282afde5b3986..10eb77615bebea405ac417f1adfe8b7b47f2515b 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -13,6 +13,9 @@ using namespace cbm::algo::ca;
class L1TrackPar {
public:
+ static constexpr int kNparTr {8};
+ static constexpr int kNparCov {28};
+
fvec x {0.}, y {0.}, tx {0.}, ty {0.}, qp {0.}, z {0.}, t {0.}, vi {0.};
fvec C00 {0.},
@@ -31,6 +34,8 @@ public:
fvec chi2 {0.}, NDF {0.};
+ fvec nTimeMeasurements {0.};
+
L1TrackPar() = default;
template<typename T>
@@ -58,7 +63,6 @@ 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;
@@ -70,20 +74,29 @@ public:
bool IsConsistent(bool printWhenWrong, int nFilled) const;
template<typename T>
- void copyToArrays(const int iVec, T p[7], T c[28]) const;
+ void copyToArrays(const int iVec, T p[kNparTr], T c[kNparCov]) const;
template<typename T>
- void copyFromArrays(const int iVec, const T p[7], const T c[28]);
+ void copyFromArrays(const int iVec, const T p[kNparTr], const T c[kNparCov]);
template<typename T>
- void copyFromArrays(const T p[7], const T c[28]);
+ void copyFromArrays(const T p[kNparTr], const T c[kNparCov]);
+
+ void InitVelocityRange(fscal minP);
+
+ static constexpr float kClightNs {29.9792458}; // the speed of light cm/ns
+ static constexpr float kClightNsInv {1. / 29.9792458}; // inverse speed of light
+ static constexpr float kProtonMass = 0.93827208816;
+ static constexpr float kPionMass = 0.1395703918;
+ static constexpr float kMuonMass = 0.105658375523;
+ static constexpr float kElectronMass = 0.0005109989500015;
} _fvecalignment;
// =============================================================================================
template<typename T>
-inline void L1TrackPar::copyToArrays(const int iVec, T p[7], T c[28]) const
+inline void L1TrackPar::copyToArrays(const int iVec, T p[kNparTr], T c[kNparCov]) const
{
p[0] = x[iVec];
p[1] = y[iVec];
@@ -95,13 +108,13 @@ inline void L1TrackPar::copyToArrays(const int iVec, T p[7], T c[28]) const
p[7] = vi[iVec];
const fvec* tC = &(C00);
- for (int i = 0; i < 28; i++) {
+ for (int i = 0; i < kNparCov; i++) {
c[i] = tC[i][iVec];
}
}
template<typename T>
-inline void L1TrackPar::copyFromArrays(const int iVec, const T p[7], const T c[28])
+inline void L1TrackPar::copyFromArrays(const int iVec, const T p[kNparTr], const T c[kNparCov])
{
x[iVec] = p[0];
y[iVec] = p[1];
@@ -113,7 +126,7 @@ inline void L1TrackPar::copyFromArrays(const int iVec, const T p[7], const T c[2
vi[iVec] = p[7];
fvec* tC = &(C00);
- for (int i = 0; i < 28; i++) {
+ for (int i = 0; i < kNparCov; i++) {
tC[i][iVec] = c[i];
}
@@ -122,7 +135,7 @@ inline void L1TrackPar::copyFromArrays(const int iVec, const T p[7], const T c[2
}
template<typename T>
-inline void L1TrackPar::copyFromArrays(const T p[7], const T c[28])
+inline void L1TrackPar::copyFromArrays(const T p[kNparTr], const T c[kNparCov])
{
x = p[0];
y = p[1];
@@ -134,7 +147,7 @@ inline void L1TrackPar::copyFromArrays(const T p[7], const T c[28])
vi = p[7];
fvec* tC = &(C00);
- for (int i = 0; i < 28; i++) {
+ for (int i = 0; i < kNparCov; i++) {
tC[i] = c[i];
}
@@ -196,15 +209,29 @@ inline void L1TrackPar::ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec
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;
+ C00 = c00;
+ C11 = c11;
+ C22 = c22;
+ C33 = c33;
+ C44 = c44;
+ C55 = c55;
+ C66 = c66;
+
chi2 = 0.;
NDF = -6.;
+ nTimeMeasurements = 0.;
+}
+
+inline void L1TrackPar::InitVelocityRange(fscal minP)
+{
+ // initialise the velocity range with respect to the minimal momentum minP {GeV/c}
+
+ fscal maxVi = sqrt(1. + (kProtonMass / minP) * (kProtonMass / minP)) * kClightNsInv;
+ fscal minVi = kClightNsInv;
+ fscal vmean = minVi + 0.4 * (maxVi - minVi);
+ fscal dvi = (maxVi - vmean) / 3.;
+ vi = vmean;
+ C66 = dvi * dvi;
}
#endif