From 7a393ba162621fb88ffd4531a7b85cf361657cb5 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Wed, 2 Nov 2022 15:17:30 +0000
Subject: [PATCH] L1: switch on the triplet refit
---
reco/L1/CbmL1.cxx | 134 ++++----------
reco/L1/L1Algo/L1Algo.h | 5 +-
reco/L1/L1Algo/L1CAIteration.cxx | 2 +
reco/L1/L1Algo/L1CAIteration.h | 35 ++--
reco/L1/L1Algo/L1CATrackFinder.cxx | 274 ++++++++++++++---------------
reco/L1/L1Algo/L1ConfigRW.cxx | 1 +
reco/L1/L1Algo/L1TrackFitter.cxx | 74 ++++----
reco/L1/L1Algo/L1TrackPar.h | 27 ++-
reco/L1/L1Algo/L1TrackParFit.cxx | 16 +-
reco/L1/L1Algo/L1TrackParFit.h | 10 +-
10 files changed, 259 insertions(+), 319 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 413f6a8cd5..2c62edb76b 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -219,12 +219,13 @@ InitStatus CbmL1::Init()
fUseTOF = 0;
fInitManager.DevSetUseOfOriginalField();
- //fInitManager.DevSetIgnoreHitSearchAreas(true);
- //fInitManager.DevSetIsMatchDoubletsViaMc(true);
- //fInitManager.DevSetIsMatchTripletsViaMc(true);
+ fInitManager.DevSetIgnoreHitSearchAreas(true);
fInitManager.SetMaxTripletPerDoublets(1000);
}
+ //fInitManager.DevSetIsMatchDoubletsViaMc(true);
+ //fInitManager.DevSetIsMatchTripletsViaMc(true);
+
CheckDetectorPresence();
fpStsPoints = 0;
@@ -634,154 +635,87 @@ InitStatus CbmL1::Init()
// TODO: Need to provide a selection: default iterations input (these hard-coded ones), input from file or input
// from running macro (S.Zharko)
- auto trackingIterFastPrim = L1CAIteration("FastPrimIter");
- trackingIterFastPrim.SetMinNhits(4);
+
+ auto trIterDefault = L1CAIteration("default");
+ trIterDefault.SetMinNhits(4);
+ trIterDefault.SetMinNhitsStation0(4);
+ trIterDefault.SetDoubletChi2Cut(7.56327f); // = 1.3449 * 2.f / 3.f; // prob = 0.1
+ trIterDefault.SetTripletChi2Cut(23.4450f); // = 7.815 * 3; // prob = 0.05
+ trIterDefault.SetTripletFinalChi2Cut(12.5);
+ trIterDefault.SetTripletLinkChi2(10.);
+ trIterDefault.SetTrackChi2Cut(10.);
+ trIterDefault.SetMaxSlope(2.748f);
+ trIterDefault.SetMaxSlopePV(1.1f);
+ trIterDefault.SetMaxDZ(0.1);
+ trIterDefault.SetPickGather(4.0f);
+ trIterDefault.SetExtendTracksFlag(true);
+
+ auto trackingIterFastPrim = L1CAIteration(trIterDefault, "FastPrimIter");
trackingIterFastPrim.SetMinNhitsStation0(3);
- trackingIterFastPrim.SetTrackChi2Cut(10.f);
- trackingIterFastPrim.SetTripletChi2Cut(23.4450f); // = 7.815 * 3; // prob = 0.05
- trackingIterFastPrim.SetDoubletChi2Cut(7.56327f); // = 1.3449 * 2.f / 3.f; // prob = 0.1
trackingIterFastPrim.SetPickGather(3.0f);
- trackingIterFastPrim.SetTripletLinkChi2(25.);
trackingIterFastPrim.SetMaxInvMom(1.0 / 0.5);
- trackingIterFastPrim.SetMaxSlopePV(1.1f);
- trackingIterFastPrim.SetMaxSlope(2.748f);
trackingIterFastPrim.SetMaxDZ(0);
trackingIterFastPrim.SetTargetPosSigmaXY(1, 1);
trackingIterFastPrim.SetPrimaryFlag(true);
- trackingIterFastPrim.SetExtendTracksFlag(true);
- auto trackingIterAllPrim = L1CAIteration("AllPrimIter");
- trackingIterAllPrim.SetMinNhits(4);
+ auto trackingIterAllPrim = L1CAIteration(trIterDefault, "AllPrimIter");
trackingIterAllPrim.SetMinNhitsStation0(3);
- trackingIterAllPrim.SetTrackChi2Cut(10.f);
- trackingIterAllPrim.SetTripletChi2Cut(23.4450f);
- trackingIterAllPrim.SetDoubletChi2Cut(7.56327f);
- trackingIterAllPrim.SetPickGather(4.0f);
- trackingIterAllPrim.SetTripletLinkChi2(25.);
trackingIterAllPrim.SetMaxInvMom(1.0 / 0.05);
- trackingIterAllPrim.SetMaxSlopePV(1.1f);
- trackingIterAllPrim.SetMaxSlope(2.748f);
- trackingIterAllPrim.SetMaxDZ(0.1);
trackingIterAllPrim.SetTargetPosSigmaXY(1, 1);
trackingIterAllPrim.SetPrimaryFlag(true);
- trackingIterAllPrim.SetExtendTracksFlag(true);
- auto trackingIterFastPrim2 = L1CAIteration("FastPrim2Iter");
- trackingIterFastPrim2.SetMinNhits(4);
- trackingIterFastPrim2.SetMinNhitsStation0(4);
- trackingIterFastPrim2.SetTrackChi2Cut(10.f);
+ auto trackingIterFastPrim2 = L1CAIteration(trIterDefault, "FastPrim2Iter");
trackingIterFastPrim2.SetTripletChi2Cut(21.1075f);
- trackingIterFastPrim2.SetDoubletChi2Cut(7.56327f);
trackingIterFastPrim2.SetPickGather(3.0f);
- trackingIterFastPrim2.SetTripletLinkChi2(25.);
trackingIterFastPrim2.SetMaxInvMom(1.0 / 0.5);
- trackingIterFastPrim2.SetMaxSlopePV(1.1f);
- trackingIterFastPrim2.SetMaxSlope(2.748f);
trackingIterFastPrim2.SetMaxDZ(0);
trackingIterFastPrim2.SetTargetPosSigmaXY(5, 5);
trackingIterFastPrim2.SetPrimaryFlag(true);
- trackingIterFastPrim2.SetExtendTracksFlag(true);
- auto trackingIterAllSec = L1CAIteration("AllSecIter");
- trackingIterAllSec.SetMinNhits(4);
- trackingIterAllSec.SetMinNhitsStation0(4);
- trackingIterAllSec.SetTrackChi2Cut(10.f);
+ auto trackingIterAllSec = L1CAIteration(trIterDefault, "AllSecIter");
trackingIterAllSec.SetTripletChi2Cut(18.7560f); // = 6.252 * 3; // prob = 0.1
- trackingIterAllSec.SetDoubletChi2Cut(7.56327f);
- trackingIterAllSec.SetPickGather(4.0f);
- trackingIterAllSec.SetTripletLinkChi2(25.);
trackingIterAllSec.SetMaxInvMom(1.0 / 0.1);
- trackingIterAllSec.SetMaxSlopePV(1.5f);
- trackingIterAllSec.SetMaxSlope(2.748f);
- trackingIterAllSec.SetMaxDZ(0.1);
trackingIterAllSec.SetTargetPosSigmaXY(10, 10);
trackingIterAllSec.SetPrimaryFlag(false);
- trackingIterAllSec.SetExtendTracksFlag(true);
- auto trackingIterFastPrimJump = L1CAIteration("FastPrimJumpIter");
- trackingIterFastPrimJump.SetMinNhits(4);
- trackingIterFastPrimJump.SetMinNhitsStation0(4);
- trackingIterFastPrimJump.SetTrackChi2Cut(10.f);
+ auto trackingIterFastPrimJump = L1CAIteration(trIterDefault, "FastPrimJumpIter");
trackingIterFastPrimJump.SetTripletChi2Cut(21.1075f); // prob = 0.01
- trackingIterFastPrimJump.SetDoubletChi2Cut(7.56327f);
trackingIterFastPrimJump.SetPickGather(3.0f);
- trackingIterFastPrimJump.SetTripletLinkChi2(25.);
trackingIterFastPrimJump.SetMaxInvMom(1.0 / 0.5);
- trackingIterFastPrimJump.SetMaxSlopePV(1.1f);
- trackingIterFastPrimJump.SetMaxSlope(2.748f);
trackingIterFastPrimJump.SetMaxDZ(0);
trackingIterFastPrimJump.SetTargetPosSigmaXY(5, 5);
trackingIterFastPrimJump.SetPrimaryFlag(true);
trackingIterFastPrimJump.SetJumpedFlag(true);
- trackingIterFastPrimJump.SetExtendTracksFlag(true);
- auto trackingIterAllPrimJump = L1CAIteration("AllPrimJumpIter");
- trackingIterAllPrimJump.SetMinNhits(4);
- trackingIterAllPrimJump.SetMinNhitsStation0(4);
- trackingIterAllPrimJump.SetTrackChi2Cut(10.f);
+ auto trackingIterAllPrimJump = L1CAIteration(trIterDefault, "AllPrimJumpIter");
trackingIterAllPrimJump.SetTripletChi2Cut(18.7560f);
- trackingIterAllPrimJump.SetDoubletChi2Cut(7.56327f);
- trackingIterAllPrimJump.SetPickGather(4.0f);
- trackingIterAllPrimJump.SetTripletLinkChi2(25.);
trackingIterAllPrimJump.SetMaxInvMom(1.0 / 0.1);
- trackingIterAllPrimJump.SetMaxSlopePV(1.1f);
- trackingIterAllPrimJump.SetMaxSlope(2.748f);
- trackingIterAllPrimJump.SetMaxDZ(0.1);
trackingIterAllPrimJump.SetTargetPosSigmaXY(5, 5);
trackingIterAllPrimJump.SetPrimaryFlag(true);
trackingIterAllPrimJump.SetJumpedFlag(true);
- trackingIterAllPrimJump.SetExtendTracksFlag(true);
- auto trackingIterAllSecJump = L1CAIteration("AllSecJumpIter");
- trackingIterAllSecJump.SetMinNhits(4);
- trackingIterAllSecJump.SetMinNhitsStation0(4);
- trackingIterAllSecJump.SetTrackChi2Cut(10.f);
+ auto trackingIterAllSecJump = L1CAIteration(trIterDefault, "AllSecJumpIter");
trackingIterAllSecJump.SetTripletChi2Cut(21.1075f);
- trackingIterAllSecJump.SetDoubletChi2Cut(7.56327f);
- trackingIterAllSecJump.SetPickGather(4.0f);
- trackingIterAllSecJump.SetTripletLinkChi2(25.);
trackingIterAllSecJump.SetMaxInvMom(1.0 / 0.1);
trackingIterAllSecJump.SetMaxSlopePV(1.5f);
- trackingIterAllSecJump.SetMaxSlope(2.748f);
- trackingIterAllSecJump.SetMaxDZ(0.1);
trackingIterAllSecJump.SetTargetPosSigmaXY(10, 10);
trackingIterAllSecJump.SetPrimaryFlag(false);
trackingIterAllSecJump.SetJumpedFlag(true);
- trackingIterAllSecJump.SetExtendTracksFlag(true);
- auto trackingIterAllPrimE = L1CAIteration("AllPrimEIter");
- trackingIterAllPrimE.SetMinNhits(4);
+ auto trackingIterAllPrimE = L1CAIteration(trIterDefault, "AllPrimEIter");
trackingIterAllPrimE.SetMinNhitsStation0(3);
- trackingIterAllPrimE.SetTrackChi2Cut(10.f);
- trackingIterAllPrimE.SetTripletChi2Cut(23.4450f);
- trackingIterAllPrimE.SetDoubletChi2Cut(7.56327f);
- trackingIterAllPrimE.SetPickGather(4.0f);
- trackingIterAllPrimE.SetTripletLinkChi2(25.);
trackingIterAllPrimE.SetMaxInvMom(1.0 / 0.05);
- trackingIterAllPrimE.SetMaxSlopePV(1.1f);
- trackingIterAllPrimE.SetMaxSlope(2.748f);
- trackingIterAllPrimE.SetMaxDZ(0.1);
trackingIterAllPrimE.SetTargetPosSigmaXY(1, 1);
trackingIterAllPrimE.SetPrimaryFlag(true);
trackingIterAllPrimE.SetElectronFlag(true);
- trackingIterAllPrimE.SetExtendTracksFlag(true);
- auto trackingIterAllSecE = L1CAIteration("AllSecEIter");
- trackingIterAllSecE.SetMinNhits(4);
- trackingIterAllSecE.SetMinNhitsStation0(4);
- trackingIterAllSecE.SetTrackChi2Cut(10.f);
+ auto trackingIterAllSecE = L1CAIteration(trIterDefault, "AllSecEIter");
trackingIterAllSecE.SetTripletChi2Cut(18.7560f);
- trackingIterAllSecE.SetDoubletChi2Cut(7.56327f);
- trackingIterAllSecE.SetPickGather(4.0f);
- trackingIterAllSecE.SetTripletLinkChi2(25.);
trackingIterAllSecE.SetMaxInvMom(1.0 / 0.05);
trackingIterAllSecE.SetMaxSlopePV(1.5f);
- trackingIterAllSecE.SetMaxSlope(2.748f);
- trackingIterAllSecE.SetMaxDZ(0.1);
trackingIterAllSecE.SetTargetPosSigmaXY(10, 10);
trackingIterAllSecE.SetPrimaryFlag(false);
trackingIterAllSecE.SetElectronFlag(true);
- trackingIterAllSecE.SetExtendTracksFlag(true);
// Select default track finder
if (L1Algo::TrackingMode::kMcbm == fTrackingMode) {
@@ -863,13 +797,13 @@ InitStatus CbmL1::Init()
auto trd2dIter3 = L1CAIteration("trd2dIter3");
trd2dIter3.SetMinNhits(4);
trd2dIter3.SetMinNhitsStation0(4);
- trd2dIter3.SetTrackChi2Cut(100 * 7.f); //10.f
- trd2dIter3.SetTripletChi2Cut(0.5 * 23.4450f); // = 7.815 * 3; // prob = 0.05
- trd2dIter3.SetDoubletChi2Cut(0.7 * 7.56327f); // = 1.3449 * 2.f / 3.f; // prob = 0.1
+ trd2dIter3.SetTrackChi2Cut(30); //10.f
+ trd2dIter3.SetTripletChi2Cut(15. * 2.); // = 7.815 * 3; // prob = 0.05
+ trd2dIter3.SetDoubletChi2Cut(15.); // = 1.3449 * 2.f / 3.f; // prob = 0.1
trd2dIter3.SetPickGather(3.0f);
- trd2dIter3.SetTripletLinkChi2(10 * 300.);
- trd2dIter3.SetMaxInvMom(1.0 / 0.05); //(1.0 / 0.5);
- trd2dIter3.SetMaxSlopePV(0 * .5f);
+ trd2dIter3.SetTripletLinkChi2(200.);
+ trd2dIter3.SetMaxInvMom(1.0 / 0.1); //(1.0 / 0.5);
+ trd2dIter3.SetMaxSlopePV(.01);
trd2dIter3.SetMaxSlope(.4); //.5f);
trd2dIter3.SetMaxDZ(0.05);
trd2dIter3.SetTargetPosSigmaXY(7 * 10, 6 * 10); //(1, 1);
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index c8fbb22ea0..f6a561ed1c 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -307,8 +307,8 @@ public:
/// Refit Triplets.
void findTripletsStep2( // input
- Tindex n3, int istal, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3, int nIterations = 0);
+ Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
+ L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3, int nIterations = 0);
/// Select triplets. Save them into vTriplets.
void findTripletsStep3( // input
@@ -528,6 +528,7 @@ private:
std::map<int, int> threadNumberToCpuMap {};
float fTrackChi2Cut {10.f};
+ float fTripletFinalChi2Cut {10.f};
float fTripletChi2Cut {5.f}; // cut for selecting triplets before collecting tracks.per one DoF
float fDoubletChi2Cut {5.f};
float fTimeCut1 {0.f}; // TODO: please, specify "1" and "2" (S.Zharko)
diff --git a/reco/L1/L1Algo/L1CAIteration.cxx b/reco/L1/L1Algo/L1CAIteration.cxx
index f643ba9ef4..5fe919b75d 100644
--- a/reco/L1/L1Algo/L1CAIteration.cxx
+++ b/reco/L1/L1Algo/L1CAIteration.cxx
@@ -32,6 +32,7 @@ bool L1CAIteration::Check() const
// debug purposes. In future, these values will be strengthened.
res = this->CheckValueLimits<float>("track_chi2_cut", fTrackChi2Cut, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("triplet_chi2_cut", fTripletChi2Cut, 0.f, kMaxFloat) && res;
+ res = this->CheckValueLimits<float>("triplet_final_chi2_cut", fTripletFinalChi2Cut, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("doublet_chi2_cut", fDoubletChi2Cut, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("pick_gather", fPickGather, 0.f, kMaxFloat) && res;
res = this->CheckValueLimits<float>("triplet_link_chi2", fTripletLinkChi2, 0.f, kMaxFloat) && res;
@@ -80,6 +81,7 @@ std::string L1CAIteration::ToString(int indentLevel) const
aStream << indent << indCh << "Min n hits for trscs at station 0: " << fMinNhitsStation0 << '\n';
aStream << indent << indCh << "Track chi2 cut: " << fTrackChi2Cut << '\n';
aStream << indent << indCh << "Triplet chi2 cut: " << fTripletChi2Cut << '\n';
+ aStream << indent << indCh << "Triplet final chi2 cut: " << fTripletFinalChi2Cut << '\n';
aStream << indent << indCh << "Doublet chi2 cut: " << fDoubletChi2Cut << '\n';
aStream << indent << indCh << "Pick gather: " << fPickGather << '\n';
aStream << indent << indCh << "Triplet link chi2: " << fTripletLinkChi2 << '\n';
diff --git a/reco/L1/L1Algo/L1CAIteration.h b/reco/L1/L1Algo/L1CAIteration.h
index 199d93a408..9892efe445 100644
--- a/reco/L1/L1Algo/L1CAIteration.h
+++ b/reco/L1/L1Algo/L1CAIteration.h
@@ -35,6 +35,9 @@ public:
/// Copy constructor
L1CAIteration(const L1CAIteration& other) = default;
+ /// Copy constructor
+ L1CAIteration(const L1CAIteration& other, const std::string& name) : L1CAIteration(other) { SetName(name); }
+
/// Move constructor
L1CAIteration(L1CAIteration&& other) noexcept = default;
@@ -110,6 +113,9 @@ public:
/// Gets triplet chi2 upper cut
float GetTripletChi2Cut() const { return fTripletChi2Cut; }
+ /// Gets triplet chi2 upper cut
+ float GetTripletFinalChi2Cut() const { return fTripletFinalChi2Cut; }
+
/// (DEBUG!) Sets flag:
/// true:
/// all the triplets found on this iteration will be converted to tracks,
@@ -196,6 +202,9 @@ public:
/// Sets triplet chi2 upper cut
void SetTripletChi2Cut(float input) { fTripletChi2Cut = input; }
+ /// Sets triplet chi2 upper cut
+ void SetTripletFinalChi2Cut(float input) { fTripletFinalChi2Cut = input; }
+
/// String representation of the class contents
/// \param indentLevel Level of indentation for the text (in terms of \t symbols)
std::string ToString(int indentLevel = 0) const;
@@ -209,18 +218,19 @@ private:
// NOTE: For each new cut one should not forget to create a setter and a getter, insert the value
// initialization in the copy constructor and the Swap operator as well as a string repre-
// sentation to the ToString method (S.Zharko)
- float fTrackChi2Cut = 10.f; ///< Track chi2 upper cut
- float fTripletChi2Cut = 21.1075f; ///< Triplet chi2 upper cut
- float fDoubletChi2Cut = 11.3449 * 2.f / 3.f; ///< Doublet chi2 upper cut
- float fPickGather = 3.0; ///< Size of region to attach new hits to the created track
- float fTripletLinkChi2 = 25.0; ///< Min value of dp^2/dp_error^2, for which two tiplets are neighbours
- float fMaxInvMom = 1.0 / 0.5; ///< Max considered q/p for tracks
- float fMaxSlopePV = 1.1; ///< Max slope (tx\ty) in primary vertex
- float fMaxSlope = 2.748; ///< Max slope (tx\ty) in 3D hit position of a triplet
- float fMaxDZ = 0.f; ///< Correction for accounting overlaping and iff z [cm]
- float fTargetPosSigmaX = 0; ///< Constraint on target position in X direction [cm]
- float fTargetPosSigmaY = 0; ///< Constraint on target position in Y direction [cm]
- int fFirstStationIndex = 0; ///< First station, used for tracking
+ float fTrackChi2Cut = 10.f; ///< Track chi2 upper cut
+ float fTripletChi2Cut = 21.1075f; ///< Triplet chi2 upper cut
+ float fTripletFinalChi2Cut = 21.1075f; ///< Triplet chi2 upper cut
+ float fDoubletChi2Cut = 11.3449 * 2.f / 3.f; ///< Doublet chi2 upper cut
+ float fPickGather = 3.0; ///< Size of region to attach new hits to the created track
+ float fTripletLinkChi2 = 25.0; ///< Min value of dp^2/dp_error^2, for which two tiplets are neighbours
+ float fMaxInvMom = 1.0 / 0.5; ///< Max considered q/p for tracks
+ float fMaxSlopePV = 1.1; ///< Max slope (tx\ty) in primary vertex
+ float fMaxSlope = 2.748; ///< Max slope (tx\ty) in 3D hit position of a triplet
+ float fMaxDZ = 0.f; ///< Correction for accounting overlaping and iff z [cm]
+ float fTargetPosSigmaX = 0; ///< Constraint on target position in X direction [cm]
+ float fTargetPosSigmaY = 0; ///< Constraint on target position in Y direction [cm]
+ int fFirstStationIndex = 0; ///< First station, used for tracking
bool fIsPrimary = false; ///< Flag: true - only primary tracks are searched for
bool fIsElectron = false; ///< Flag: true - only electrons are searched for
@@ -250,6 +260,7 @@ private:
ar& fName;
ar& fTrackChi2Cut;
ar& fTripletChi2Cut;
+ ar& fTripletFinalChi2Cut;
ar& fDoubletChi2Cut;
ar& fPickGather;
ar& fTripletLinkChi2;
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 20ec8e43c1..d885033dc8 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -21,17 +21,16 @@
#include "L1Algo.h"
#include "L1Assert.h"
+#include "L1Branch.h"
#include "L1Extrapolation.h"
#include "L1Filtration.h"
#include "L1Fit.h"
-#include "L1HitPoint.h"
-#include "L1Track.h"
-#include "L1TrackPar.h"
-//#include "TDHelper.h"
-#include "L1Branch.h"
#include "L1Grid.h"
#include "L1HitArea.h"
+#include "L1HitPoint.h"
#include "L1Portion.h"
+#include "L1Track.h"
+#include "L1TrackPar.h"
#ifdef _OPENMP
#include "omp.h"
@@ -80,19 +79,16 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
if (r.GetMSta() != l.GetRSta()) return false;
if (r.GetLSta() != l.GetMSta()) return false;
- fscal dqp = fabs(l.GetQp() - r.GetQp());
- fscal Cqp = l.GetCqp() + r.GetCqp();
-
- fscal dtx = fabs(l.GetTx() - r.GetTx());
- fscal Ctx = l.GetCtx() + r.GetCtx();
-
- fscal dty = fabs(l.GetTy() - r.GetTy());
- fscal Cty = l.GetCty() + r.GetCty();
if (r.IsMomentumFitted()) {
assert(l.IsMomentumFitted());
+
+ fscal dqp = l.GetQp() - r.GetQp();
+ fscal Cqp = l.GetCqp() + r.GetCqp();
+
if (!std::isfinite(dqp)) return false;
if (!std::isfinite(Cqp)) return false;
+
if (dqp * dqp > fTripletLinkChi2 * Cqp) {
return false; // bad neighbour // CHECKME why do we need recheck it?? (it really change result)
}
@@ -100,6 +96,12 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
}
else {
+ fscal dtx = l.GetTx() - r.GetTx();
+ fscal Ctx = l.GetCtx() + r.GetCtx();
+
+ fscal dty = l.GetTy() - r.GetTy();
+ fscal Cty = l.GetCty() + r.GetCty();
+
// it shouldn't happen, but happens sometimes
if (!std::isfinite(dtx)) return false;
@@ -953,26 +955,38 @@ inline void L1Algo::findTripletsStep1( // input
}
-/// Refit Triplets.
-inline void L1Algo::findTripletsStep2( // input // TODO not updated after gaps introduction
- Tindex n3, int istal, L1Vector<L1TrackPar>& T_3, L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
- L1Vector<L1HitIndex_t>& hitsr_3, int nIterations)
+inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar, L1Vector<L1TrackPar>& T_3,
+ L1Vector<L1HitIndex_t>& hitsl_3, L1Vector<L1HitIndex_t>& hitsm_3,
+ L1Vector<L1HitIndex_t>& hitsr_3, int nIterations)
{
+
+ nIterations = 2;
+
+ /// Refit Triplets
+
L1Fit fit;
- fit.SetParticleMass(fDefaultMass);
+ fit.SetParticleMass(GetDefaultParticleMass());
- const int NHits = 3; // triplets
+ L1TrackParFit fitNew;
+ fitNew.SetParticleMass(GetDefaultParticleMass());
- // TODO: SG: middle and right station numbers might be different for different triplets!!!
+ const int NHits = 3; // triplets
// prepare data
- int ista[NHits] = {istal, istal + 1, istal + 2};
+ int ista[NHits] = {istal, istam, istar};
L1Station sta[3];
for (int is = 0; is < NHits; ++is) {
sta[is] = fParameters.GetStation(ista[is]);
};
+ const L1Station& stal = fParameters.GetStation(istal);
+ const L1Station& stam = fParameters.GetStation(istam);
+ const L1Station& star = fParameters.GetStation(istar);
+
+ bool isMomentumFitted = ((stal.fieldStatus != 0) || (stam.fieldStatus != 0) || (star.fieldStatus != 0));
+ fvec ndf = isMomentumFitted ? -5 : -4;
+
for (int i3 = 0; i3 < n3; ++i3) {
int i3_V = i3 / fvec::size();
int i3_4 = i3 % fvec::size();
@@ -981,8 +995,6 @@ inline void L1Algo::findTripletsStep2( // input // TODO not updated after gaps
//if (!T3.IsEntryConsistent(false, i3_4)) continue;
- fscal qp0 = T3.qp[i3_4];
-
// prepare data
L1HitIndex_t ihit[NHits] = {(*RealIHitP)[hitsl_3[i3] + HitsUnusedStartIndex[ista[0]]],
(*RealIHitP)[hitsm_3[i3] + HitsUnusedStartIndex[ista[1]]],
@@ -995,41 +1007,23 @@ inline void L1Algo::findTripletsStep2( // input // TODO not updated after gaps
if ((mc1 != mc2) || (mc1 != mc3)) { continue; }
}
- fscal u[NHits], v[NHits], x[NHits], y[NHits], z[NHits];
+ fscal u[NHits], v[NHits], t[NHits], x[NHits], y[NHits], z[NHits];
+ fscal du2[NHits], dv2[NHits], dt2[NHits];
+
for (int ih = 0; ih < NHits; ++ih) {
const L1Hit& hit = fInputData.GetHit(ihit[ih]);
u[ih] = hit.u;
v[ih] = hit.v;
+ t[ih] = hit.t;
std::tie(x[ih], y[ih]) = sta[ih].ConvUVtoXY<fscal>(u[ih], v[ih]);
z[ih] = hit.z;
+ du2[ih] = hit.du * hit.du;
+ dv2[ih] = hit.dv * hit.dv;
+ dt2[ih] = hit.dt * hit.dt;
};
- // initialize parameters
- L1TrackPar T;
-
- const fvec vINF = .1;
- T.x = x[0];
- T.y = y[0];
-
- fvec dz01 = 1. / (z[1] - z[0]);
- T.tx = (x[1] - x[0]) * dz01;
- T.ty = (y[1] - y[0]) * dz01;
+ // find the field along the track
- T.qp = qp0;
- T.z = z[0];
- T.chi2 = 0.;
- T.NDF = 2.;
- T.C00 = sta[0].XYInfo.C00;
- T.C10 = sta[0].XYInfo.C10;
- T.C11 = sta[0].XYInfo.C11;
-
- T.C20 = T.C21 = 0;
- T.C30 = T.C31 = T.C32 = 0;
- T.C40 = T.C41 = T.C42 = T.C43 = 0;
- T.C22 = T.C33 = vINF;
- T.C44 = 1.;
-
- // find field along the track
L1FieldValue B[3] _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -1042,89 +1036,90 @@ inline void L1Algo::findTripletsStep2( // input // TODO not updated after gaps
fld.Set(B[0], sta[0].z, B[1], sta[1].z, B[2], sta[2].z);
- // fit
- for (int ih = 1; ih < NHits; ++ih) {
- L1Extrapolate(T, z[ih], T.qp, fld);
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista[ih], T.x, T.y), T.qp, fvec::One());
- }
- else {
- fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, fvec::One());
- }
- if (ista[ih] == fNstationsBeforePipe - 1) { fit.L1AddPipeMaterial(T, T.qp, fvec::One()); }
- L1Filter(T, sta[ih].frontInfo, u[ih]);
- L1Filter(T, sta[ih].backInfo, v[ih]);
+
+ L1TrackPar& T = fitNew.fTr;
+
+ // initial parameters
+ {
+ fvec dz01 = 1. / (z[1] - z[0]);
+ T.tx = (x[1] - x[0]) * dz01;
+ T.ty = (y[1] - y[0]) * dz01;
+ T.qp = 0.; //(fscal) T3.qp[i3_4];
}
- // repeat several times in order to improve precision
+ // repeat several times in order to improve the precision
for (int iiter = 0; iiter < nIterations; ++iiter) {
- // fit backward
- // keep tx,ty,q/p
- int ih = NHits - 1;
- T.x = x[ih];
- T.y = y[ih];
- T.z = z[ih];
- T.chi2 = 0.;
- T.NDF = 2.;
- T.C00 = sta[ih].XYInfo.C00;
- T.C10 = sta[ih].XYInfo.C10;
- T.C11 = sta[ih].XYInfo.C11;
-
- T.C20 = T.C21 = 0;
- T.C30 = T.C31 = T.C32 = 0;
- T.C40 = T.C41 = T.C42 = T.C43 = 0;
- T.C22 = T.C33 = vINF;
- T.C44 = 1.;
-
- // L1Filter( T, sta[ih].frontInfo, u[ih] );
- // L1Filter( T, sta[ih].backInfo, v[ih] );
- for (ih = NHits - 2; ih >= 0; ih--) {
- L1Extrapolate(T, z[ih], T.qp, fld);
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista[ih], T.x, T.y), T.qp, fvec::One());
- }
- else {
- fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, fvec::One());
- }
- if (ista[ih] == fNstationsBeforePipe) fit.L1AddPipeMaterial(T, T.qp, fvec::One());
-
- L1Filter(T, sta[ih].frontInfo, u[ih]);
- L1Filter(T, sta[ih].backInfo, v[ih]);
- }
// fit forward
- ih = 0;
- T.x = x[ih];
- T.y = y[ih];
- T.z = z[ih];
- T.chi2 = 0.;
- T.NDF = 2.;
- T.C00 = sta[ih].XYInfo.C00;
- T.C10 = sta[ih].XYInfo.C10;
- T.C11 = sta[ih].XYInfo.C11;
-
- T.C20 = T.C21 = 0;
- T.C30 = T.C31 = T.C32 = 0;
- T.C40 = T.C41 = T.C42 = T.C43 = 0;
- T.C22 = T.C33 = vINF;
- T.C44 = 1.;
-
- // L1Filter( T, sta[ih].frontInfo, u[ih] );
- // L1Filter( T, sta[ih].backInfo, v[ih] );
- for (ih = 1; ih < NHits; ++ih) {
- L1Extrapolate(T, z[ih], T.qp, fld);
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista[ih], T.x, T.y), T.qp, fvec::One());
- }
- else {
- fit.L1AddMaterial(T, sta[ih].materialInfo, T.qp, fvec::One());
+ {
+ fvec qp0 = T.qp;
+ if (qp0[0] > GetMaxInvMom()) { qp0 = GetMaxInvMom(); }
+ if (qp0[0] < -GetMaxInvMom()) { qp0 = -GetMaxInvMom(); }
+
+ T.ResetErrors(200., 200., 1., 1., 100., 1.e4);
+ T.NDF = ndf;
+ T.qp = 0.;
+ int ih0 = 0;
+ T.x = x[ih0];
+ T.y = y[ih0];
+ T.z = z[ih0];
+ T.t = t[ih0];
+ T.C55 = dt2[ih0];
+
+ fitNew.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One());
+ fitNew.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
+
+ for (int ih = 1; ih < NHits; ++ih) {
+ fitNew.Extrapolate(z[ih], qp0, fld, fvec::One());
+ if constexpr (L1Constants::control::kIfUseRadLengthTable) {
+ fitNew.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
+ }
+ else {
+ fitNew.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
+ }
+ if (ista[ih] == fNstationsBeforePipe) { fitNew.AddPipeMaterial(qp0, fvec::One()); }
+ fitNew.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
+ fitNew.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
}
- if (ista[ih] == fNstationsBeforePipe + 1) fit.L1AddPipeMaterial(T, T.qp, fvec::One());
+ }
- L1Filter(T, sta[ih].frontInfo, u[ih]);
- L1Filter(T, sta[ih].backInfo, v[ih]);
+ if (iiter == nIterations - 1) break;
+
+ // fit backward
+ {
+ fvec qp0 = T.qp;
+ if (qp0[0] > GetMaxInvMom()) { qp0 = GetMaxInvMom(); }
+ if (qp0[0] < -GetMaxInvMom()) { qp0 = -GetMaxInvMom(); }
+
+ T.ResetErrors(200., 200., 1., 1., 100., 1.e4);
+ T.NDF = ndf;
+ T.qp = 0.;
+ int ih0 = NHits - 1;
+ T.x = x[ih0];
+ T.y = y[ih0];
+ T.z = z[ih0];
+ T.t = t[ih0];
+ T.C55 = dt2[ih0];
+
+ fitNew.Filter(sta[ih0].frontInfo, u[ih0], du2[ih0], fvec::One());
+ fitNew.Filter(sta[ih0].backInfo, v[ih0], dv2[ih0], fvec::One());
+
+ for (int ih = NHits - 2; ih >= 0; --ih) {
+ fitNew.Extrapolate(z[ih], qp0, fld, fvec::One());
+ if constexpr (L1Constants::control::kIfUseRadLengthTable) {
+ fitNew.AddMaterial(fParameters.GetMaterialThickness(ista[ih], T.x, T.y), qp0, fvec::One());
+ }
+ else {
+ fitNew.AddMaterial(sta[ih].materialInfo, qp0, fvec::One());
+ }
+ if (ista[ih] == fNstationsBeforePipe - 1) { fitNew.AddPipeMaterial(qp0, fvec::One()); }
+ fitNew.Filter(sta[ih].frontInfo, u[ih], du2[ih], fvec::One());
+ fitNew.Filter(sta[ih].backInfo, v[ih], dv2[ih], fvec::One());
+ }
}
} // for iiter
+ //cout << " chi2 " << T3.chi2[i3_4] << " " << T.chi2[0] << endl;
+ //T.chi2 = (fscal) T3.chi2[i3_4];
T3.SetOneEntry(i3_4, T, i3_4);
} //i3
} // findTripletsStep2
@@ -1192,7 +1187,7 @@ inline void L1Algo::findTripletsStep3( // input
ihitl_prev = 1 + hitsl_3[i3];
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
- if (chi2 > fTripletChi2Cut) { continue; }
+ if (chi2 > fTripletFinalChi2Cut) { continue; }
}
// assert(std::isfinite(chi2) && chi2 > 0);
@@ -1475,7 +1470,7 @@ inline void L1Algo::CreatePortionOfTriplets(
L1HitPoint* vHits_r = 0;
vHits_r = &((*vHitPointsUnused)[0]) + HitsUnusedStartIndex[istar];
- Tindex n3 = 0, n3_V;
+ Tindex n3 = 0;
/// Add the middle hits to parameters estimation. Propagate to right station.
@@ -1523,8 +1518,6 @@ inline void L1Algo::CreatePortionOfTriplets(
n3, T_3, hitsl_3, hitsm_3, hitsr_3, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER);
- n3_V = (n3 + fvec::size() - 1) / fvec::size();
-
for (Tindex i = 0; i < static_cast<Tindex>(hitsl_3.size()); ++i)
L1_assert(hitsl_3[i] < HitsUnusedStopIndex[istal] - HitsUnusedStartIndex[istal]);
for (Tindex i = 0; i < static_cast<Tindex>(hitsm_3.size()); ++i)
@@ -1533,14 +1526,16 @@ inline void L1Algo::CreatePortionOfTriplets(
L1_assert(hitsr_3[i] < HitsUnusedStopIndex[istar] - HitsUnusedStartIndex[istar]);
/// Add the right hits to parameters estimation.
+ /*
+ Tindex n3_V = (n3 + fvec::size() - 1) / fvec::size();
findTripletsStep1( // input
n3_V, star, u_front3, u_back3, z_pos3, du3, dv3, timeR, timeER,
// output
T_3);
-
+*/
/// refit
- // findTripletsStep2( n3, istal, _RealIHit, T_3, hitsl_3, hitsm_3, hitsr_3, 0 );
+ findTripletsStep2(n3, istal, istam, istar, T_3, hitsl_3, hitsm_3, hitsr_3, 1);
#ifdef TRIP_PERFORMANCE
L1HitIndex_t* RealIHitL = &((*RealIHitP)[HitsUnusedStartIndex[istal]]);
@@ -1555,7 +1550,7 @@ inline void L1Algo::CreatePortionOfTriplets(
#else
fL1Eff_triplets->AddOne(iHits);
#endif
- if (T_3[i_V].chi2[i_4] < fTripletChi2Cut) { fL1Eff_triplets2->AddOne(iHits); }
+ if (T_3[i_V].chi2[i_4] < fTripletFinalChi2Cut) { fL1Eff_triplets2->AddOne(iHits); }
}
#endif // TRIP_PERFORMANCE
@@ -1840,15 +1835,16 @@ void L1Algo::CATrackFinder()
{
// --- SET PARAMETERS FOR THE ITERATION ---
- fFirstCAstation = caIteration.GetFirstStationIndex();
- fTrackChi2Cut = caIteration.GetTrackChi2Cut();
- fDoubletChi2Cut = caIteration.GetDoubletChi2Cut(); //11.3449 * 2.f / 3.f; // prob = 0.1
- fTripletChi2Cut = caIteration.GetTripletChi2Cut(); //21.1075; // prob = 0.01%
- fPickGather = caIteration.GetPickGather(); //3.0;
- fTripletLinkChi2 = caIteration.GetTripletLinkChi2(); //5.0;
- fMaxInvMom = caIteration.GetMaxInvMom(); //1.0 / 0.5; // max considered q/p
- fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
- fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
+ fFirstCAstation = caIteration.GetFirstStationIndex();
+ fTrackChi2Cut = caIteration.GetTrackChi2Cut();
+ fDoubletChi2Cut = caIteration.GetDoubletChi2Cut(); //11.3449 * 2.f / 3.f; // prob = 0.1
+ fTripletChi2Cut = caIteration.GetTripletChi2Cut(); //21.1075; // prob = 0.01%
+ fTripletFinalChi2Cut = caIteration.GetTripletFinalChi2Cut();
+ fPickGather = caIteration.GetPickGather(); //3.0;
+ fTripletLinkChi2 = caIteration.GetTripletLinkChi2(); //5.0;
+ fMaxInvMom = caIteration.GetMaxInvMom(); //1.0 / 0.5; // max considered q/p
+ fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
+ fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
// define the target
fTargX = fParameters.GetTargetPositionX();
diff --git a/reco/L1/L1Algo/L1ConfigRW.cxx b/reco/L1/L1Algo/L1ConfigRW.cxx
index 79c5ac68fe..562df5ef82 100644
--- a/reco/L1/L1Algo/L1ConfigRW.cxx
+++ b/reco/L1/L1Algo/L1ConfigRW.cxx
@@ -63,6 +63,7 @@ void L1ConfigRW::ReadCAIterations(const YAML::Node& node)
auto caIter = L1CAIteration(input["name"].as<std::string>());
caIter.SetTrackChi2Cut(input["track_chi2_cut"].as<float>(caIter.GetTrackChi2Cut()));
caIter.SetTripletChi2Cut(input["triplet_chi2_cut"].as<float>(caIter.GetTripletChi2Cut()));
+ caIter.SetTripletFinalChi2Cut(input["triplet_final_chi2_cut"].as<float>(caIter.GetTripletFinalChi2Cut()));
caIter.SetDoubletChi2Cut(input["doublet_chi2_cut"].as<float>(caIter.GetDoubletChi2Cut()));
caIter.SetPickGather(input["pick_gather"].as<float>(caIter.GetPickGather()));
caIter.SetTripletLinkChi2(input["triplet_link_chi2"].as<float>(caIter.GetTripletLinkChi2()));
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 5aa4f50863..bfcc7d40fa 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -540,26 +540,19 @@ void L1Algo::L1KFTrackFitter()
fit.Extrapolate(z[ista], qp01, fld1, wExtr);
if (ista == fNstationsBeforePipe - 1) {
- fit.L1AddPipeMaterial(qp01, wExtr);
+ fit.AddPipeMaterial(qp01, wExtr);
fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(1.f), wExtr);
}
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(1.f), wExtr);
}
else {
- fit.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr);
+ fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr);
}
- L1UMeasurementInfo info = sta[ista].frontInfo;
- info.sigma2 = d_u[ista] * d_u[ista];
-
- fit.Filter(info, u[ista], w1);
-
- info = sta[ista].backInfo;
- info.sigma2 = d_v[ista] * d_v[ista];
-
- fit.Filter(info, v[ista], w1);
+ fit.Filter(sta[ista].frontInfo, u[ista], d_u[ista] * d_u[ista], w1);
+ fit.Filter(sta[ista].backInfo, v[ista], d_v[ista] * d_v[ista], w1);
fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo);
fldB2 = fldB1;
@@ -588,8 +581,8 @@ void L1Algo::L1KFTrackFitter()
fitpv = fit;
fitpv.fTr.qp = qp00;
fitpv.Extrapolate(fParameters.GetTargetPositionZ(), fitpv.fTr.qp, fld, fvec(1.));
- fitpv.Filter(vtxInfoX, fParameters.GetTargetPositionX());
- fitpv.Filter(vtxInfoY, fParameters.GetTargetPositionY());
+ fitpv.Filter(vtxInfoX, fParameters.GetTargetPositionX(), fvec(1.e-8), fvec::One());
+ fitpv.Filter(vtxInfoY, fParameters.GetTargetPositionY(), fvec(1.e-8), fvec::One());
qp00 = fitpv.fTr.qp;
}
}
@@ -704,26 +697,19 @@ void L1Algo::L1KFTrackFitter()
fit.Extrapolate(z[ista], qp01, fld, w1);
if (ista == fNstationsBeforePipe) {
- fit.L1AddPipeMaterial(qp01, wExtr);
+ fit.AddPipeMaterial(qp01, wExtr);
fit.EnergyLossCorrection(fit.fPipeRadThick, qp01, fvec(-1.f), wExtr);
}
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- fit.L1AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
+ fit.AddMaterial(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, wExtr);
fit.EnergyLossCorrection(fParameters.GetMaterialThickness(ista, tr.x, tr.y), qp01, fvec(-1.f), wExtr);
}
else {
- fit.L1AddMaterial(sta[ista].materialInfo, qp01, wExtr);
+ fit.AddMaterial(sta[ista].materialInfo, qp01, wExtr);
}
- L1UMeasurementInfo info = sta[ista].frontInfo;
- info.sigma2 = d_u[ista] * d_u[ista];
-
- fit.Filter(info, u[ista], w1);
-
- info = sta[ista].backInfo;
- info.sigma2 = d_v[ista] * d_v[ista];
-
- fit.Filter(info, v[ista], w1);
+ fit.Filter(sta[ista].frontInfo, u[ista], d_u[ista] * d_u[ista], w1);
+ fit.Filter(sta[ista].backInfo, v[ista], d_v[ista] * d_v[ista], w1);
fit.FilterTime(time[ista], timeEr[ista], w1_time, sta[ista].timeInfo);
fldB2 = fldB1;
@@ -1014,7 +1000,7 @@ void L1Algo::L1KFTrackFitterMuch()
fit.L1AddPipeMaterial(T, qp0, wIn);
fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
- T1.L1AddPipeMaterial(qp01, wIn);
+ T1.AddPipeMaterial(qp01, wIn);
T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(-1.f), wIn);
}
@@ -1026,8 +1012,8 @@ void L1Algo::L1KFTrackFitterMuch()
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(-1.f), wIn);
- T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick,
- 1);
+ T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick,
+ 1);
}
else {
// L1AddMaterial( T, sta[i].materialInfo, qp0, wIn );
@@ -1036,13 +1022,13 @@ void L1Algo::L1KFTrackFitterMuch()
L1UMeasurementInfo info = sta[i].frontInfo;
info.sigma2 = d_xx[i];
L1Filter(T, info, u[i], w1);
- T1.Filter(info, u[i], w1);
+ T1.Filter(info, u[i], d_xx[i], w1);
info = sta[i].backInfo;
info.sigma2 = d_yy[i];
L1Filter(T, info, v[i], w1);
- T1.Filter(info, v[i], w1);
+ T1.Filter(info, v[i], d_yy[i], w1);
T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
}
@@ -1100,8 +1086,8 @@ void L1Algo::L1KFTrackFitterMuch()
T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
fvec(-1.f), wIn);
- T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, wIn,
- sta[i].materialInfo.thick / (nofSteps + fvec(1)), 1);
+ T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, wIn,
+ sta[i].materialInfo.thick / (nofSteps + fvec(1)), 1);
wIn = iif(!maskLastStep, wIn, fvec::Zero());
}
@@ -1115,12 +1101,12 @@ void L1Algo::L1KFTrackFitterMuch()
L1UMeasurementInfo info = sta[i].frontInfo;
info.sigma2 = d_xx[i];
L1Filter(T, info, u[i], w1);
- T1.Filter(info, u[i], w1);
+ T1.Filter(info, u[i], d_xx[i], w1);
info = sta[i].backInfo;
info.sigma2 = d_yy[i];
L1Filter(T, info, v[i], w1);
- T1.Filter(info, v[i], w1);
+ T1.Filter(info, v[i], d_yy[i], w1);
T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
}
}
@@ -1181,7 +1167,7 @@ void L1Algo::L1KFTrackFitterMuch()
// T1.EnergyLossCorrectionAl( fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(1.f), wIn);
//
- // T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, ONE, sta[i].materialInfo.thick, 0);
+ // T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, ONE, sta[i].materialInfo.thick, 0);
for (--i; i >= 0; i--) {
@@ -1236,8 +1222,8 @@ void L1Algo::L1KFTrackFitterMuch()
T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
fvec(1.f), w2);
- T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, w2,
- sta[i].materialInfo.thick / (nofSteps + fvec(1)), 0);
+ T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, w2,
+ sta[i].materialInfo.thick / (nofSteps + fvec(1)), 0);
w2 = iif(!maskLastStep, w2, fvec::Zero());
}
@@ -1257,13 +1243,13 @@ void L1Algo::L1KFTrackFitterMuch()
info.sigma2 = d_xx[i];
L1Filter(T, info, u[i], w1);
- T1.Filter(info, u[i], w1);
+ T1.Filter(info, u[i], d_xx[i], w1);
info = sta[i].backInfo;
info.sigma2 = d_yy[i];
L1Filter(T, info, v[i], w1);
- T1.Filter(info, v[i], w1);
+ T1.Filter(info, v[i], d_yy[i], w1);
T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
}
@@ -1300,8 +1286,8 @@ void L1Algo::L1KFTrackFitterMuch()
if constexpr (L1Constants::control::kIfUseRadLengthTable) {
T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(1.f), wIn);
- T1.L1AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick,
- 0);
+ T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick,
+ 0);
}
else {
fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
@@ -1309,11 +1295,11 @@ void L1Algo::L1KFTrackFitterMuch()
L1UMeasurementInfo info = sta[i].frontInfo;
// info.sigma2 = d_u[i] * d_u[i];
- T1.Filter(info, u[i], w1);
+ T1.Filter(info, u[i], info.sigma2, w1);
info = sta[i].backInfo;
// info.sigma2 = d_v[i] * d_v[i];
- T1.Filter(info, v[i], w1);
+ T1.Filter(info, v[i], info.sigma2, w1);
T1.FilterTime(time[i], timeEr[i], w1, sta[i].timeInfo);
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index b0c922bc91..7afaa11752 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -16,7 +16,6 @@ public:
fvec x, y, tx, ty, qp, z, t, C00, C10, C11, C20, C21, C22, C30, C31, C32, C33, C40, C41, C42, C43, C44, C50, C51, C52,
C53, C54, C55, chi2, NDF;
-
L1TrackPar()
: x(0)
, y(0)
@@ -48,6 +47,7 @@ public:
, C55(0)
, chi2(0)
, NDF(0) {};
+
L1TrackPar(double* T, double* C)
: x(T[0])
, y(T[1])
@@ -97,14 +97,23 @@ public:
bool IsConsistent(bool printWhenWrong, int nFilled) const;
- // void L1Extrapolate
- // (
- // // L1TrackPar &T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
- // fvec z_out , // extrapolate to this z position
- // fvec qp0 , // use Q/p linearisation at this value
- // L1FieldRegion &F
- // );
-
+ void ResetErrors(fvec c00, fvec c11, fvec c22, fvec c33, fvec c44, fvec c55)
+ {
+ C10 = 0.;
+ C20 = C21 = 0.;
+ C30 = C31 = C32 = 0.;
+ C40 = C41 = C42 = C43 = 0.;
+ C50 = C51 = C52 = C53 = C54 = 0.;
+
+ C00 = c00;
+ C11 = c11;
+ C22 = c22;
+ C33 = c33;
+ C44 = c44;
+ C55 = c55;
+ chi2 = 0.;
+ NDF = 0;
+ }
} _fvecalignment;
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index 7c40030a5b..fd310acdd4 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -9,7 +9,7 @@
#define cnst const fvec
-void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w)
+void L1TrackParFit::Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& sigma2, const fvec& w)
{
fvec zeta, HCH;
fvec F0, F1, F2, F3, F4, F5;
@@ -28,14 +28,14 @@ void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w)
F4 = info.cos_phi * fTr.C40 + info.sin_phi * fTr.C41;
F5 = info.cos_phi * fTr.C50 + info.sin_phi * fTr.C51;
- const fmask maskDoFilter = (HCH < info.sigma2 * 16.f);
+ const fmask maskDoFilter = (HCH < sigma2 * 16.f);
//const fvec maskDoFilter = _f32vec4_true;
// 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 + fvec(1.0000001) * HCH);
- fvec zetawi = w * zeta / (iif(maskDoFilter, info.sigma2, fvec::Zero()) + HCH);
+ fvec wi = w / (sigma2 + fvec(1.0000001) * HCH);
+ fvec zetawi = w * zeta / (iif(maskDoFilter, sigma2, fvec::Zero()) + HCH);
// fTr.chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() );
fTr.chi2 += zeta * zeta * wi;
@@ -792,7 +792,7 @@ void
fTr.C55 = cj55 + cj25 * J[17] + cj35 * J[23] + cj45 * J[29];
}
-void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w)
+void L1TrackParFit::AddPipeMaterial(fvec qp0, fvec w)
{
cnst ONE = 1.f;
@@ -822,7 +822,7 @@ void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w)
}
-void L1TrackParFit::L1AddMaterial(const fvec& radThick, fvec qp0, fvec w)
+void L1TrackParFit::AddMaterial(const fvec& radThick, fvec qp0, fvec w)
{
cnst ONE = 1.;
@@ -847,7 +847,7 @@ void L1TrackParFit::L1AddMaterial(const fvec& radThick, fvec qp0, fvec w)
fTr.C33 += w * (ONE + tyty) * a;
}
-void L1TrackParFit::L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream)
+void L1TrackParFit::AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream)
{
cnst ONE = 1.;
@@ -887,7 +887,7 @@ void L1TrackParFit::L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thi
}
-void L1TrackParFit::L1AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w)
+void L1TrackParFit::AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w)
{
cnst ONE = 1.f;
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index e66c7eb0b0..cb965a9b47 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -47,7 +47,7 @@ public:
/// get the particle mass squared
fvec GetParticleMass2() const { return fMass2; }
- void Filter(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
+ void Filter(const L1UMeasurementInfo& info, const fvec& u, const fvec& sigma2, const fvec& w);
void FilterXY(const L1XYMeasurementInfo& info, fvec x, fvec y);
void FilterTime(fvec t0, fvec dt0, fvec w = 1., fvec timeInfo = 1.);
void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
@@ -70,12 +70,12 @@ public:
void EnergyLossCorrectionCarbon(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
void EnergyLossCorrectionAl(const fvec& radThick, fvec& qp0, fvec direction, fvec w);
- void L1AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w);
+ void AddMaterial(const L1MaterialInfo& info, fvec qp0, fvec w);
- void L1AddMaterial(const fvec& radThick, fvec qp0, fvec w = 1);
+ void AddMaterial(const fvec& radThick, fvec qp0, fvec w = 1);
- void L1AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
- void L1AddPipeMaterial(fvec qp0, fvec w = 1);
+ void AddThickMaterial(fvec radThick, fvec qp0, fvec w, fvec thickness, bool fDownstream);
+ void AddPipeMaterial(fvec qp0, fvec w);
// void L1Extrapolate
// (
--
GitLab