diff --git a/core/base/CbmTrackingDetectorInterfaceBase.cxx b/core/base/CbmTrackingDetectorInterfaceBase.cxx
index 997182b13d3b98efa3f69ea2456257e4214251f4..6fd5b1bd5d872aa959d08d44c7dfc30c61dbd1cc 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.cxx
+++ b/core/base/CbmTrackingDetectorInterfaceBase.cxx
@@ -82,20 +82,6 @@ bool CbmTrackingDetectorInterfaceBase::Check() const
res = false && res;
}
- // Front strips spatial resolution
- auto rmsF = this->GetStripsSpatialRmsFront(iSt);
- if (rmsF < std::numeric_limits<double>::epsilon() || std::isnan(rmsF)) {
- msg << prefix << " zero, negative or NaN front strips spatial resolution (" << rmsF << " cm)\n";
- res = false && res;
- }
-
- // Back strips spatial resolution
- auto rmsB = this->GetStripsSpatialRmsBack(iSt);
- if (rmsB < std::numeric_limits<double>::epsilon() || std::isnan(rmsB)) {
- msg << prefix << " zero, negative or NaN back strips spatial resolution (" << rmsB << " cm)\n";
- res = false && res;
- }
-
// Time resolution
auto timeRes = this->GetTimeResolution(iSt);
if (timeRes < std::numeric_limits<double>::epsilon() || std::isnan(timeRes)) {
@@ -145,8 +131,6 @@ std::string CbmTrackingDetectorInterfaceBase::ToString() const
table << setw(12) << setfill(' ') << "res.time[ns]" << ' ';
table << setw(11) << setfill(' ') << "angleF[rad]" << ' ';
table << setw(11) << setfill(' ') << "angleB[rad]" << ' ';
- table << setw(10) << setfill(' ') << "res.F [cm]" << ' ';
- table << setw(10) << setfill(' ') << "res.B [cm]" << ' ';
table << setw(10) << setfill(' ') << "dz [cm]" << ' ';
table << setw(10) << setfill(' ') << "RL [cm]" << '\n';
for (int iSt = 0; iSt < GetNtrackingStations(); ++iSt) {
@@ -159,8 +143,6 @@ std::string CbmTrackingDetectorInterfaceBase::ToString() const
table << setw(12) << setfill(' ') << GetTimeResolution(iSt) << ' ';
table << setw(11) << setfill(' ') << GetStripsStereoAngleFront(iSt) << ' ';
table << setw(11) << setfill(' ') << GetStripsStereoAngleBack(iSt) << ' ';
- table << setw(10) << setfill(' ') << GetStripsSpatialRmsFront(iSt) << ' ';
- table << setw(10) << setfill(' ') << GetStripsSpatialRmsBack(iSt) << ' ';
table << setw(10) << setfill(' ') << GetThickness(iSt) << ' ';
table << setw(10) << setfill(' ') << GetRadLength(iSt) << '\n';
}
diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h
index 2c402a314f258235d3bea8d760cd640c3657df22..ad1c9fb5a7b5827405804e20887f2d7c1fa16a3b 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.h
+++ b/core/base/CbmTrackingDetectorInterfaceBase.h
@@ -58,16 +58,6 @@ public:
/// \return Absolute stereo angle for back strips [rad]
virtual double GetStripsStereoAngleBack(int stationId) const = 0;
- /// Gets spatial resolution (RMS) for front strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- virtual double GetStripsSpatialRmsFront(int stationId) const = 0;
-
- /// Gets spatial resolution (RMS) for back strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- virtual double GetStripsSpatialRmsBack(int stationId) const = 0;
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
diff --git a/core/detectors/much/CbmMuchTrackingInterface.h b/core/detectors/much/CbmMuchTrackingInterface.h
index 260c79d1d11555f37173914d41624ffdf32e9894..348d40e0723836fe6022a7279b996d1073d1e161 100644
--- a/core/detectors/much/CbmMuchTrackingInterface.h
+++ b/core/detectors/much/CbmMuchTrackingInterface.h
@@ -72,16 +72,6 @@ public:
/// \return Size of station inner radius [cm]
double GetRmin(int /*stationId*/) const { return 10.; }
- /// Gets spatial resolution (RMS) for back strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsBack(int /*stationId*/) const { return 0.35; }
-
- /// Gets spatial resolution (RMS) for front strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsFront(int /*stationId*/) const { return 0.35; }
-
/// Gets back strips stereo angle
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Absolute stereo angle for back strips [rad]
diff --git a/core/detectors/sts/CbmStsTrackingInterface.h b/core/detectors/sts/CbmStsTrackingInterface.h
index edc503e37e1994ffc0e1c1f1849177bb835c37d5..72af70ac787212c07ba382bd1b5bc3ab87841049 100644
--- a/core/detectors/sts/CbmStsTrackingInterface.h
+++ b/core/detectors/sts/CbmStsTrackingInterface.h
@@ -71,22 +71,6 @@ public:
/// \return Size of station inner radius [cm]
double GetRmin(int /*stationId*/) const { return 0.; }
- /// Gets spatial resolution (RMS) for back strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsBack(int stationId) const
- {
- return GetStsStation(stationId)->GetSensorPitch(0) / TMath::Sqrt(12.);
- }
-
- /// Gets spatial resolution (RMS) for front strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsFront(int stationId) const
- {
- return GetStsStation(stationId)->GetSensorPitch(0) / TMath::Sqrt(12.);
- }
-
/// Gets back strips stereo angle
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Absolute stereo angle for back strips [rad]
diff --git a/core/detectors/tof/CbmTofTrackingInterface.h b/core/detectors/tof/CbmTofTrackingInterface.h
index 2f57e16a4402850a0cdd6f99041d50deb08bbdea..c94fa74a84db6a9338d31ddb9186c544992cbd13 100644
--- a/core/detectors/tof/CbmTofTrackingInterface.h
+++ b/core/detectors/tof/CbmTofTrackingInterface.h
@@ -78,16 +78,6 @@ public:
/// \return Absolute stereo angle for front strips [rad]
double GetStripsStereoAngleFront(int /*stationId*/) const { return 0.; }
- /// Gets spatial resolution (RMS) for back strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsBack(int /*stationId*/) const { return 0.23; }
-
- /// Gets spatial resolution (RMS) for front strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsFront(int /*stationId*/) const { return 0.42; }
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
diff --git a/core/detectors/trd/CbmTrdTrackingInterface.h b/core/detectors/trd/CbmTrdTrackingInterface.h
index 197d4dcde85e597365c0dc6c0edbe87f7f71c22c..f6e94dfff38844dd793d6b70d13347749b7930df 100644
--- a/core/detectors/trd/CbmTrdTrackingInterface.h
+++ b/core/detectors/trd/CbmTrdTrackingInterface.h
@@ -78,16 +78,6 @@ public:
/// \return Absolute stereo angle for front strips [rad]
double GetStripsStereoAngleFront(int /*stationId*/) const { return 0.; }
- /// Gets spatial resolution (RMS) for back strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsBack(int /*stationId*/) const { return 0.15; }
-
- /// Gets spatial resolution (RMS) for front strips
- /// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
- /// \return Spatial resolution (RMS) for front strips [cm]
- double GetStripsSpatialRmsFront(int /*stationId*/) const { return 0.15; }
-
/// Gets station thickness along the Z-axis
/// \param stationId Tracking station ID in the setup (NOTE: must be in range [0..GetNstations()-1])
/// \return Station thickness [cm]
diff --git a/macro/L1/run_reco_L1global.C b/macro/L1/run_reco_L1global.C
index 979e42aaf602ab93b86e5074fe40d60bb38c63e0..05368ebdffab82641f0e6b8caf1025520b1ad539 100644
--- a/macro/L1/run_reco_L1global.C
+++ b/macro/L1/run_reco_L1global.C
@@ -359,7 +359,6 @@ void run_reco_L1global(TString input = "", Int_t nTimeSlices = -1, Int_t firstTi
l1 = new CbmL1("L1", 0);
}
l1->SetGlobalMode();
- //l1->SetUseHitErrors(false);
//l1->SetUseMcHit(0, 0, 0, 1, 0);
// --- Material budget file names
diff --git a/macro/beamtime/mcbm2021/mcbm_event_reco.C b/macro/beamtime/mcbm2021/mcbm_event_reco.C
index 469b49544af48335d8a348d9d03b0637c25574cf..f40a29949e6090422836d647c913b292dd729b5c 100644
--- a/macro/beamtime/mcbm2021/mcbm_event_reco.C
+++ b/macro/beamtime/mcbm2021/mcbm_event_reco.C
@@ -810,7 +810,6 @@ Bool_t mcbm_event_reco(UInt_t uRunId = 1588,
// CbmL1* l1 = new CbmL1();
// l1->SetLegacyEventMode(1);
// l1->SetMcbmMode();
- // l1->SetUseHitErrors(1);
// if (strcmp(geoSetupTag.data(), "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1);
//
// // --- Material budget file names
diff --git a/macro/beamtime/mcbm2022/mcbm_digievent_reco.C b/macro/beamtime/mcbm2022/mcbm_digievent_reco.C
index 6ff2580c8d409d5e3fe23f5ca2079a1b38566a20..75f695b2917ddc8b72c68bc600247e9db356bf07 100644
--- a/macro/beamtime/mcbm2022/mcbm_digievent_reco.C
+++ b/macro/beamtime/mcbm2022/mcbm_digievent_reco.C
@@ -459,7 +459,6 @@ Bool_t mcbm_digievent_reco(UInt_t uRunId = 2365,
// CbmL1* l1 = new CbmL1();
// l1->SetLegacyEventMode(1);
// l1->SetMcbmMode();
- // l1->SetUseHitErrors(1);
// if (strcmp(geoSetupTag.data(), "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1);
//
// // --- Material budget file names
diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco.C b/macro/beamtime/mcbm2022/mcbm_event_reco.C
index 46c5eeeba201408b69e00d65d6a9b43738eb2981..e4640c7645f99a4d851b08d96287c7fe052ec979 100644
--- a/macro/beamtime/mcbm2022/mcbm_event_reco.C
+++ b/macro/beamtime/mcbm2022/mcbm_event_reco.C
@@ -861,7 +861,6 @@ Bool_t mcbm_event_reco(UInt_t uRunId = 2391,
CbmL1* l1 = new CbmL1();
l1->SetLegacyEventMode(1);
l1->SetMcbmMode();
- l1->SetUseHitErrors(1);
// --- Material budget file names
TString mvdGeoTag;
diff --git a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
index 301766e970938704a8422dc615aff1758f834d7f..bb957967f51d7c1c8a39f4c57422b5bbad570912 100644
--- a/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
+++ b/macro/beamtime/mcbm2022/mcbm_event_reco_L1.C
@@ -569,7 +569,7 @@ Bool_t mcbm_event_reco_L1(UInt_t uRunId = 2391,
CbmL1* l1 = new CbmL1();
l1->SetLegacyEventMode(1);
l1->SetMcbmMode();
- l1->SetUseHitErrors(1);
+
// if (strcmp(geoSetupTag.data(), "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1);
// --- Material budget file names
diff --git a/macro/mcbm/mcbm_reco.C b/macro/mcbm/mcbm_reco.C
index f67bc7f87f6a782b6f902a4cfb0f83c11d5ec60a..a5c09f0c8c76a1c6b0be27db1def758efedc813b 100644
--- a/macro/mcbm/mcbm_reco.C
+++ b/macro/mcbm/mcbm_reco.C
@@ -431,7 +431,6 @@ void mcbm_reco(Int_t nEvents = 10, TString dataset = "data/test", TString sEvBui
CbmL1* l1 = new CbmL1();
l1->SetLegacyEventMode(1);
l1->SetMcbmMode();
- l1->SetUseHitErrors(1);
// --- Material budget file names
TString mvdGeoTag;
diff --git a/macro/mcbm/mcbm_reco_event.C b/macro/mcbm/mcbm_reco_event.C
index 447f3f25c671446cbaed362b9111f41215424c1d..fe761acf33a32074f181024f92eda7c0c1bcb41d 100644
--- a/macro/mcbm/mcbm_reco_event.C
+++ b/macro/mcbm/mcbm_reco_event.C
@@ -378,7 +378,7 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test",
auto l1 = (debugWithMC) ? new CbmL1("L1", 1, 3, 0) : new CbmL1();
l1->SetLegacyEventMode(1);
l1->SetMcbmMode();
- l1->SetUseHitErrors(1);
+
if (strcmp(setupName, "mcbm_beam_2021_07_surveyed") == 0) l1->SetMissingHits(1);
// --- Material budget file names
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 359a43e5ad52c3f717885e8fb39ab09f854d3bfa..7373ce8608e264a68d0c890382bcd78fa2c04292 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -34,6 +34,7 @@
#include <boost/filesystem.hpp>
// TODO: include of CbmSetup.h creates problems on Mac
// #include "CbmSetup.h"
+#include "CbmDigiManager.h"
#include "CbmMCDataObject.h"
#include "CbmStsFindTracks.h"
#include "CbmStsHit.h"
@@ -198,9 +199,10 @@ InitStatus CbmL1::Init()
{
CbmStsFindTracks* findTask = L1_DYNAMIC_CAST<CbmStsFindTracks*>(FairRunAna::Instance()->GetTask("STSFindTracks"));
if (findTask) fUseMVD = findTask->MvdUsage();
+ CbmDigiManager::Instance()->Init();
+ if (!CbmDigiManager::IsPresent(ECbmModuleId::kMvd)) { fUseMVD = false; }
}
-
if (L1Algo::TrackingMode::kMcbm == fTrackingMode) {
fUseMUCH = 1;
fUseTRD = 1;
@@ -474,9 +476,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZthickness(mvdInterface->GetThickness(iSt));
stationInfo.SetMaterialMap(std::move(materialTableMvd[iSt]));
// TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
- stationInfo.SetFrontBackStripsGeometry(
- (fscal) mvdInterface->GetStripsStereoAngleFront(iSt), (fscal) mvdInterface->GetStripsSpatialRmsFront(iSt),
- (fscal) mvdInterface->GetStripsStereoAngleBack(iSt), (fscal) mvdInterface->GetStripsSpatialRmsBack(iSt));
+ stationInfo.SetFrontBackStripsGeometry((fscal) mvdInterface->GetStripsStereoAngleFront(iSt),
+ (fscal) mvdInterface->GetStripsStereoAngleBack(iSt));
stationInfo.SetTrackingStatus(target.z < stationInfo.GetZdouble() ? true : false);
fInitManager.AddStation(stationInfo);
LOG(info) << "- MVD station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
@@ -502,9 +503,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZthickness(stsInterface->GetThickness(iSt));
stationInfo.SetMaterialMap(std::move(materialTableSts[iSt]));
// TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
- stationInfo.SetFrontBackStripsGeometry(
- (fscal) stsInterface->GetStripsStereoAngleFront(iSt), (fscal) stsInterface->GetStripsSpatialRmsFront(iSt),
- (fscal) stsInterface->GetStripsStereoAngleBack(iSt), (fscal) stsInterface->GetStripsSpatialRmsBack(iSt));
+ stationInfo.SetFrontBackStripsGeometry((fscal) stsInterface->GetStripsStereoAngleFront(iSt),
+ (fscal) stsInterface->GetStripsStereoAngleBack(iSt));
stationInfo.SetTrackingStatus(target.z < stationInfo.GetZdouble() ? true : false);
fInitManager.AddStation(stationInfo);
LOG(info) << "- STS station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
@@ -530,9 +530,8 @@ InitStatus CbmL1::Init()
stationInfo.SetZthickness(muchInterface->GetThickness(iSt));
stationInfo.SetMaterialMap(std::move(materialTableMuch[iSt]));
// TODO: The CA TF result is dependent from type of geometry settings. Should be understood (S.Zharko)
- stationInfo.SetFrontBackStripsGeometry(
- (fscal) muchInterface->GetStripsStereoAngleFront(iSt), (fscal) muchInterface->GetStripsSpatialRmsFront(iSt),
- (fscal) muchInterface->GetStripsStereoAngleBack(iSt), (fscal) muchInterface->GetStripsSpatialRmsBack(iSt));
+ stationInfo.SetFrontBackStripsGeometry((fscal) muchInterface->GetStripsStereoAngleFront(iSt),
+ (fscal) muchInterface->GetStripsStereoAngleBack(iSt));
stationInfo.SetTrackingStatus(target.z < stationInfo.GetZdouble() ? true : false);
fInitManager.AddStation(stationInfo);
LOG(info) << "- MuCh station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
@@ -557,17 +556,13 @@ InitStatus CbmL1::Init()
stationInfo.SetRmax(trdInterface->GetRmax(iSt));
stationInfo.SetZthickness(trdInterface->GetThickness(iSt));
stationInfo.SetMaterialMap(std::move(materialTableTrd[iSt]));
- fscal trdFrontPhi = trdInterface->GetStripsStereoAngleFront(iSt);
- fscal trdBackPhi = trdInterface->GetStripsStereoAngleBack(iSt);
- fscal trdFrontSigma = trdInterface->GetStripsSpatialRmsFront(iSt);
- fscal trdBackSigma = trdInterface->GetStripsSpatialRmsBack(iSt);
+ fscal trdFrontPhi = trdInterface->GetStripsStereoAngleFront(iSt);
+ fscal trdBackPhi = trdInterface->GetStripsStereoAngleBack(iSt);
if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
- trdFrontSigma = 1.1;
- trdBackSigma = 1.1;
// stationInfo.SetTimeResolution(1.e10);
stationInfo.SetTimeInfo(false);
}
- stationInfo.SetFrontBackStripsGeometry(trdFrontPhi, trdFrontSigma, trdBackPhi, trdBackSigma);
+ stationInfo.SetFrontBackStripsGeometry(trdFrontPhi, trdBackPhi);
stationInfo.SetTrackingStatus(target.z < stationInfo.GetZdouble() ? true : false);
if (iSt == 1 && L1Algo::TrackingMode::kMcbm == fTrackingMode && fMissingHits) {
stationInfo.SetTrackingStatus(false);
@@ -596,11 +591,9 @@ InitStatus CbmL1::Init()
stationInfo.SetYmax(tofInterface->GetYmax(iSt));
stationInfo.SetRmin(tofInterface->GetRmin(iSt));
stationInfo.SetRmax(tofInterface->GetRmax(iSt));
- fscal tofFrontPhi = tofInterface->GetStripsStereoAngleFront(iSt);
- fscal tofBackPhi = tofInterface->GetStripsStereoAngleBack(iSt);
- fscal tofFrontSigma = tofInterface->GetStripsSpatialRmsFront(iSt);
- fscal tofBackSigma = tofInterface->GetStripsSpatialRmsBack(iSt);
- stationInfo.SetFrontBackStripsGeometry(tofFrontPhi, tofFrontSigma, tofBackPhi, tofBackSigma);
+ fscal tofFrontPhi = tofInterface->GetStripsStereoAngleFront(iSt);
+ fscal tofBackPhi = tofInterface->GetStripsStereoAngleBack(iSt);
+ stationInfo.SetFrontBackStripsGeometry(tofFrontPhi, tofBackPhi);
stationInfo.SetTrackingStatus(target.z < stationInfo.GetZdouble() ? true : false);
fInitManager.AddStation(stationInfo);
LOG(info) << "- TOF station " << iSt << " at z = " << stationInfo.GetZdouble() << " cm";
@@ -845,7 +838,7 @@ InitStatus CbmL1::Init()
// Init L1 algo core
fpAlgo = &gAlgo;
- fpAlgo->Init(fUseHitErrors, fTrackingMode, fMissingHits);
+ fpAlgo->Init(fTrackingMode, fMissingHits);
//
// ** Send formed parameters object to L1Algo instance **
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 849f04022c681da2c97cd268effca4b9eb643b16..676a65878c3fae74629fc4d801c0e36103fea4ac 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -294,7 +294,6 @@ public:
void SetExtrapolateToTheEndOfSTS(bool b) { fExtrapolateToTheEndOfSTS = b; }
void SetLegacyEventMode(bool b) { fLegacyEventMode = b; }
- void SetUseHitErrors(bool value) { fUseHitErrors = value; }
void SetMissingHits(bool value) { fMissingHits = value; }
void SetStsOnlyMode() { fTrackingMode = L1Algo::TrackingMode::kSts; }
void SetMcbmMode() { fTrackingMode = L1Algo::TrackingMode::kMcbm; }
@@ -462,7 +461,6 @@ public:
L1InitManager fInitManager; ///< Tracking parameters data manager
L1IODataManager fIODataManager; ///< Input-output data manager
- bool fUseHitErrors = true; ///<
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 da1c6784a8fb6ed002d20b6a5c90c4f2533115ca..7fe6e9bf0d715d37d35875e3b5da0de329976d26 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -1964,21 +1964,11 @@ void CbmL1::InputPerformance()
mcPos.SetY(pt->GetYIn() + t * (pt->GetYOut() - pt->GetYIn()));
mcPos.SetZ(hitPos.Z());
- if (0) { // standard errors
- if (hitErr.X() > 1.e-8) pullXsts->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
- if (hitErr.Y() > 1.e-8) pullYsts->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
- }
- else if (1) { // qa errors
+ { // qa errors
if (sh->GetDx() > 1.e-8) pullXsts->Fill((hitPos.X() - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYsts->Fill((hitPos.Y() - mcPos.Y()) / sh->GetDy());
pullTsts->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
}
- else { // errors used in TF
- pullXsts->Fill((hitPos.X() - mcPos.X())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C00[0]));
- pullYsts->Fill((hitPos.Y() - mcPos.Y())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C11[0]));
- }
resXsts->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
resYsts->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
@@ -2018,15 +2008,8 @@ void CbmL1::InputPerformance()
mcPos.SetY((pt->GetY() + pt->GetYOut()) / 2.);
mcPos.SetZ(hitPos.Z());
- // if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / hitErr.X() ); // standard errors
- // if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / hitErr.Y() );
- // if (hitErr.X() != 0) pullX->Fill( (hitPos.X() - mcPos.X()) / sh->GetDx() ); // qa errors
- // if (hitErr.Y() != 0) pullY->Fill( (hitPos.Y() - mcPos.Y()) / sh->GetDy() );
- if (hitErr.X() != 0)
- pullXmvd->Fill((hitPos.X() - mcPos.X())
- / sqrt(fpAlgo->GetParameters()->GetStation(0).XYInfo.C00[0])); // errors used in TF
- if (hitErr.Y() != 0)
- pullYmvd->Fill((hitPos.Y() - mcPos.Y()) / sqrt(fpAlgo->GetParameters()->GetStation(0).XYInfo.C11[0]));
+ if (hitErr.X() != 0) pullXmvd->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
+ if (hitErr.Y() != 0) pullYmvd->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
resXmvd->Fill((hitPos.X() - mcPos.X()) * 10 * 1000);
resYmvd->Fill((hitPos.Y() - mcPos.Y()) * 10 * 1000);
@@ -2083,21 +2066,11 @@ void CbmL1::InputPerformance()
mcPos.SetY(0.5 * (pt->GetYIn() + pt->GetYOut()));
mcPos.SetZ(hitPos.Z());
- if (0) { // standard errors
- if (hitErr.X() > 1.e-8) pullXmuch->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
- if (hitErr.Y() > 1.e-8) pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
- }
- else if (1) { // qa errors
+ { // qa errors
if (sh->GetDx() > 1.e-8) pullXmuch->Fill((h.x - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYmuch->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTmuch->Fill((h.t - mcTime) / sh->GetTimeError());
}
- else { // errors used in TF
- pullXmuch->Fill((hitPos.X() - mcPos.X())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C00[0]));
- pullYmuch->Fill((hitPos.Y() - mcPos.Y())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C11[0]));
- }
resXmuch->Fill((h.x - mcPos.X()) * 10 * 1000);
resYmuch->Fill((h.y - mcPos.Y()) * 10 * 1000);
@@ -2156,21 +2129,11 @@ void CbmL1::InputPerformance()
mcPos.SetY((pt->GetYIn() + pt->GetYOut()) / 2.);
mcPos.SetZ(hitPos.Z());
- if (0) { // standard errors
- if (hitErr.X() > 1.e-8) pullXtrd->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
- if (hitErr.Y() > 1.e-8) pullYtrd->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
- }
- else if (1) { // qa errors
+ { // qa errors
if (sh->GetDx() > 1.e-8) pullXtrd->Fill((h.x - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYtrd->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTtrd->Fill((h.t - mcTime) / sh->GetTimeError());
}
- else { // errors used in TF
- pullXtrd->Fill((hitPos.X() - mcPos.X())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C00[0]));
- pullYtrd->Fill((hitPos.Y() - mcPos.Y())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C11[0]));
- }
resXtrd->Fill((h.x - mcPos.X()) * 10 * 1000);
resYtrd->Fill((h.y - mcPos.Y()) * 10 * 1000);
@@ -2230,21 +2193,11 @@ void CbmL1::InputPerformance()
mcPos.SetY((pt->GetY()));
mcPos.SetZ(hitPos.Z());
- if (0) { // standard errors
- if (hitErr.X() > 1.e-8) pullXmuch->Fill((hitPos.X() - mcPos.X()) / hitErr.X());
- if (hitErr.Y() > 1.e-8) pullYmuch->Fill((hitPos.Y() - mcPos.Y()) / hitErr.Y());
- }
- else if (1) { // qa errors
+ { // qa errors
if (sh->GetDx() > 1.e-8) pullXtof->Fill((h.x - mcPos.X()) / sh->GetDx());
if (sh->GetDy() > 1.e-8) pullYtof->Fill((h.y - mcPos.Y()) / sh->GetDy());
pullTtof->Fill((sh->GetTime() - mcTime) / sh->GetTimeError());
}
- else { // errors used in TF
- pullXtof->Fill((hitPos.X() - mcPos.X())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C00[0]));
- pullYtof->Fill((hitPos.Y() - mcPos.Y())
- / sqrt(fpAlgo->GetParameters()->GetStation(fNMvdStations).XYInfo.C11[0]));
- }
resXtof->Fill((h.x - mcPos.X()) * 10 * 1000);
resYtof->Fill((h.y - mcPos.Y()) * 10 * 1000);
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 79493eeafc19b561a88853bae0b5c4a01419f4ad..d4d7132fbcc44ba659bfc3960de073614fbdd19d 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -114,7 +114,7 @@ struct TmpHit {
/// \param NStrips
/// \param st reference to the station info object
void CreateHitFromPoint(const CbmL1MCPoint& point, int ip, int det, int nTmpHits, int& NStrips, const L1Station& st,
- bool doSmear = 1)
+ double du_, double dv_, bool doSmear)
{
ExtIndex = 0;
Det = det;
@@ -128,12 +128,15 @@ struct TmpHit {
iStripB = iStripF;
NStrips++;
- dx = sqrt(st.XYInfo.C00[0]);
- dy = sqrt(st.XYInfo.C11[0]);
- dxy = st.XYInfo.C10[0];
+ du = du_;
+ dv = dv_;
- du = sqrt(st.frontInfo.sigma2[0]);
- dv = sqrt(st.backInfo.sigma2[0]);
+ double c00, c10, c11;
+ std::tie(c00, c10, c11) = st.FormXYCovarianceMatrix(du_ * du_, dv_ * dv_);
+
+ dx = sqrt(c00);
+ dy = sqrt(c11);
+ dxy = c10;
SetHitFromPoint(point, ip, st, doSmear);
}
@@ -739,7 +742,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
const CbmL1MCPoint& p = fvMCPoints[ip];
TmpHit th;
int DetId = 0;
- th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation));
+ double du = 5.e-4;
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du,
+ du, true);
tmpHits.push_back(th);
nMvdHits++;
}
@@ -850,7 +855,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
// << p.time << " mc " << p.ID << " p " << p.p << endl;
TmpHit th;
int DetId = 1;
- th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation));
+ double du = 10.e-4;
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du,
+ du, true);
tmpHits.push_back(th);
nStsHits++;
}
@@ -963,7 +970,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
TmpHit th;
int DetId = 2;
- th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation));
+ double du = 100.e-4;
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du,
+ du, true);
tmpHits.push_back(th);
nMuchHits++;
@@ -1094,7 +1103,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
// const CbmL1MCTrack& t = fvMCTracks[mcTrack];
TmpHit th;
int DetId = 3;
- th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation));
+ double du = 0.1;
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du,
+ du, true);
tmpHits.push_back(th);
nTrdHits++;
}
@@ -1188,7 +1199,9 @@ void CbmL1::ReadEvent(float& TsStart, float& TsLength, float& /*TsOverlap*/, int
TmpHit th;
int DetId = 4;
- th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation));
+ double du = 0.1;
+ th.CreateHitFromPoint(p, ip, DetId, tmpHits.size(), NStrips, fpAlgo->GetParameters()->GetStation(p.iStation), du,
+ du, true);
tmpHits.push_back(th);
nTofHits++;
}
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 4699109dc3a9d6f96b4f7915fb1a42135f693c74..d6522a659f7e8a36a72253ab3dcc192107309e35 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -60,7 +60,7 @@ void L1Algo::SetNThreads(unsigned int n)
}
-void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool MissingHits)
+void L1Algo::Init(const TrackingMode mode, const bool MissingHits)
{
for (int iProc = 0; iProc < 4; iProc++) {
for (int i = 0; i < 8; i++) {
@@ -73,7 +73,6 @@ void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool M
}
}
- fUseHitErrors = UseHitErrors;
fTrackingMode = mode;
fMissingHits = MissingHits;
}
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 4079694853097b53f5fc9679ff947d17e4d70c2f..b7a7b54b07b30a38b301d39119f09056874ad3a1 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -344,14 +344,7 @@ public:
void GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur = 0,
fvec* timeV = 0, fvec* w_time = 0);
- void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st);
- void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st);
- void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st);
-
- void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& dx2, fvec& dy2,
- fvec& dxy);
- void FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& dx2, fvec& dy2,
- fvec& dxy);
+ void FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fvec& dx2, fvec& dy2, fvec& dxy);
#ifdef DRAW
@@ -366,7 +359,7 @@ public:
kMcbm
};
- void Init(const bool UseHitErrors, const TrackingMode mode, const bool MissingHits);
+ void Init(const TrackingMode mode, const bool MissingHits);
void Finish();
@@ -423,6 +416,8 @@ public:
L1Grid vGrid[L1Constants::size::kMaxNstations]; ///<
L1Grid vGridTime[L1Constants::size::kMaxNstations]; ///<
+ fscal fMaxDx[L1Constants::size::kMaxNstations];
+ fscal fMaxDy[L1Constants::size::kMaxNstations];
double fCATime {0.}; // time of track finding
@@ -465,7 +460,6 @@ public:
L1Vector<int> fStripToTrackB {"L1Algo::fStripToTrackB"}; // back strip to track pointers
int fNThreads {0};
- bool fUseHitErrors {true};
bool fMissingHits {0}; ///< TODO ???
TrackingMode fTrackingMode {kSts};
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.cxx b/reco/L1/L1Algo/L1BaseStationInfo.cxx
index ba2c03518e0a1ac4bbedae9960b641073370a97e..01a93f090b03b2fe1ffbe53b25b13647c2ea2b21 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.cxx
+++ b/reco/L1/L1Algo/L1BaseStationInfo.cxx
@@ -292,7 +292,7 @@ void L1BaseStationInfo::SetFieldStatus(int fieldStatus)
//------------------------------------------------------------------------------------------------------------------------
//
-void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double frontSigma, double backPhi, double backSigma)
+void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double backPhi)
{
//----- Original code from L1Algo -----------------------------------------------------------------------
double cFront = cos(frontPhi);
@@ -303,34 +303,23 @@ void L1BaseStationInfo::SetFrontBackStripsGeometry(double frontPhi, double front
// NOTE: Here additional double variables are used to save the precission
fL1Station.frontInfo.cos_phi = cFront;
fL1Station.frontInfo.sin_phi = sFront;
- fL1Station.frontInfo.sigma2 = frontSigma * frontSigma;
fL1Station.backInfo.cos_phi = cBack;
fL1Station.backInfo.sin_phi = sBack;
- fL1Station.backInfo.sigma2 = backSigma * backSigma;
double det = cFront * sBack - sFront * cBack;
assert(fabs(det) > 1.e-2);
- double det2 = det * det;
-
- fL1Station.XYInfo.C00 = (sBack * sBack * frontSigma * frontSigma + sFront * sFront * backSigma * backSigma) / det2;
- fL1Station.XYInfo.C10 = -(sBack * cBack * frontSigma * frontSigma + sFront * cFront * backSigma * backSigma) / det2;
- fL1Station.XYInfo.C11 = (cBack * cBack * frontSigma * frontSigma + cFront * cFront * backSigma * backSigma) / det2;
fL1Station.xInfo.cos_phi = sBack / det;
fL1Station.xInfo.sin_phi = -sFront / det;
- fL1Station.xInfo.sigma2 = fL1Station.XYInfo.C00;
fL1Station.yInfo.cos_phi = -cBack / det;
fL1Station.yInfo.sin_phi = cFront / det;
- fL1Station.yInfo.sigma2 = fL1Station.XYInfo.C11;
//-------------------------------------------------------------------------------------------------------
fInitController.SetFlag(EInitKey::kStripsFrontPhi);
- fInitController.SetFlag(EInitKey::kStripsFrontSigma);
fInitController.SetFlag(EInitKey::kStripsBackPhi);
- fInitController.SetFlag(EInitKey::kStripsBackSigma);
}
//------------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1BaseStationInfo.h b/reco/L1/L1Algo/L1BaseStationInfo.h
index 0f3a51fad3e0d95736d960e2301e8487e472f5cc..b896f27f96603290e34f64e70cd07136edd97278 100644
--- a/reco/L1/L1Algo/L1BaseStationInfo.h
+++ b/reco/L1/L1Algo/L1BaseStationInfo.h
@@ -46,20 +46,18 @@ public:
kXmax, ///< max size in X direction
kYmax, ///< max size in Y direction
// L1Station initialization
- kType, ///< station type
- kTimeInfo, ///< if time info is used (flag)
- kFieldStatus, ///< if station is placed in field (flag)
- kZ, ///< z coordinate of the station position
- kRmin, ///< internal radius of station (gap size)
- kRmax, ///< exteranl radius of station
- kZthickness, ///< Z thickness of the station
- kThicknessMap, ///< thickness map of the station (optional?)
- kFieldSlice, ///< L1Station.L1FieldSlice object initialization
- kStripsFrontPhi, ///< strips geometry initialization
- kStripsFrontSigma, ///<
- kStripsBackPhi, ///<
- kStripsBackSigma, ///<
- kTimeResolution, ///< time resolution
+ kType, ///< station type
+ kTimeInfo, ///< if time info is used (flag)
+ kFieldStatus, ///< if station is placed in field (flag)
+ kZ, ///< z coordinate of the station position
+ kRmin, ///< internal radius of station (gap size)
+ kRmax, ///< exteranl radius of station
+ kZthickness, ///< Z thickness of the station
+ kThicknessMap, ///< thickness map of the station (optional?)
+ kFieldSlice, ///< L1Station.L1FieldSlice object initialization
+ kStripsFrontPhi, ///< strips geometry initialization
+ kStripsBackPhi, ///<
+ kTimeResolution, ///< time resolution
// The last item is equal to the number of bits in fInitFlags
kEnd
};
@@ -183,12 +181,10 @@ public:
/// of magnetic field components in position
void SetFieldFunction(const std::function<void(const double (&xyz)[3], double (&B)[3])>& getFieldValue);
- /// Sets stereo angles and sigmas for front and back strips
+ /// Sets stereo angles for front and back strips
/// \param f_phi Stereoangle of front strips [rad]
- /// \param f_sigma Sigma of front strips [rad]
/// \param b_phi Stereoangle of back strips [rad]
- /// \param b_sigma Sigma of back strips [rad]
- void SetFrontBackStripsGeometry(double fPhi, double fSigma, double bPhi, double bSigma);
+ void SetFrontBackStripsGeometry(double fPhi, double bPhi);
/// Sets station thickness and radiation length
/// \param thickness Thickness of station [arb. units]
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index 2c3d46bf79d8909b45de9e12363d55080f2ce4fa..bb1cb0d445531a01a4e31ec41f57861ed7234b5b 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -84,11 +84,8 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
// T.t[0]=(hit0.t+hit1.t+hit2.t)/3;
T.chi2 = fvec(0.);
T.NDF = fvec(2.);
- T.C00 = sta0.XYInfo.C00;
- T.C10 = sta0.XYInfo.C10;
- T.C11 = sta0.XYInfo.C11;
- if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2); }
+ std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2);
T.C20 = T.C21 = 0;
T.C30 = T.C31 = T.C32 = 0;
@@ -138,16 +135,8 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
fvec u = hit.u;
fvec v = hit.v;
- L1UMeasurementInfo info = sta.frontInfo;
-
- if (fUseHitErrors) { info.sigma2 = hit.du2; }
- L1Filter(T, info, u);
-
- info = sta.backInfo;
-
- if (fUseHitErrors) { info.sigma2 = hit.dv2; }
- L1Filter(T, info, v);
-
+ L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
+ L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
FilterTime(T, hit.t, hit.dt2);
fldB0 = fldB1;
@@ -251,9 +240,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const fscal iz = 1.f / (T.z[0] - fParameters.GetTargetPositionZ()[0]);
L1HitAreaTime area(vGridTime[ista], T.x[0] * iz, T.y[0] * iz,
- (sqrt(fPickGather * (T.C00 + sta.XYInfo.C00)) + fMaxDZ * abs(T.tx))[0] * iz,
- (sqrt(fPickGather * (T.C11 + sta.XYInfo.C11)) + fMaxDZ * abs(T.ty))[0] * iz, T.t[0],
- sqrt(T.C55[0]));
+ (sqrt(fPickGather * T.C00) + fMaxDx[ista] + fMaxDZ * abs(T.tx))[0] * iz,
+ (sqrt(fPickGather * T.C11) + fMaxDy[ista] + fMaxDZ * abs(T.ty))[0] * iz, T.t[0], sqrt(T.C55[0]));
for (L1HitIndex_t ih = -1; true;) { // loop over the hits in the area
@@ -294,7 +282,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
fscal d2 = d_x * d_x + d_y * d_y;
if (d2 > r2_best) continue;
- fscal dxm_est2 = (pickGather2 * (abs(C00 + sta.XYInfo.C00)))[0];
+ fscal dxm_est2 = (pickGather2 * (abs(C00 + fMaxDx[ista] * fMaxDx[ista])))[0];
if (d_x * d_x > dxm_est2) continue;
r2_best = d2;
@@ -317,16 +305,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
L1ExtrapolateLine(T, z);
fit.L1AddMaterial(T, fParameters.GetMaterialThickness(ista, T.x, T.y), qp0, fvec::One());
- L1UMeasurementInfo info = sta.frontInfo;
-
- if (fUseHitErrors) { info.sigma2 = hit.du2; }
- L1Filter(T, info, u);
-
- info = sta.backInfo;
-
- if (fUseHitErrors) { info.sigma2 = hit.dv2; }
- L1Filter(T, info, v);
-
+ L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
+ L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
FilterTime(T, hit.t, hit.dt2);
fldB0 = fldB1;
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 2a8c9e8b684ff50fc4a42449b745e753ed4467d4..cd8d74c14415a183f20595bb2851044ee4356da1 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -162,10 +162,8 @@ inline void L1Algo::findSingletsStep0( // input
u_front_l[i1_V][i1_4] = hitl.U();
u_back_l[i1_V][i1_4] = hitl.V();
- if (fUseHitErrors) {
- du2_l[i1_V][i1_4] = hitl.dU2();
- dv2_l[i1_V][i1_4] = hitl.dV2();
- }
+ du2_l[i1_V][i1_4] = hitl.dU2();
+ dv2_l[i1_V][i1_4] = hitl.dV2();
zPos_l[i1_V][i1_4] = hitl.Z();
}
@@ -285,14 +283,11 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
T.C55 = timeEr2;
{ // add the target constraint
- T.x = xl;
- T.y = yl;
- T.z = zl;
- T.C00 = stal.XYInfo.C00;
- T.C10 = stal.XYInfo.C10;
- T.C11 = stal.XYInfo.C11;
+ T.x = xl;
+ T.y = yl;
+ T.z = zl;
- if (fUseHitErrors) { std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]); }
+ std::tie(T.C00, T.C10, T.C11) = stal.FormXYCovarianceMatrix(du2_l[i1_V], dv2_l[i1_V]);
//assert(T.IsConsistent(true, -1));
@@ -403,8 +398,8 @@ inline void L1Algo::findDoubletsStep0(
const fscal time = T1.t[i1_4];
L1HitAreaTime areaTime(vGridTime[iStaM], T1.x[i1_4] * iz, T1.y[i1_4] * iz,
- (sqrt(Pick_m22 * (T1.C00 + stam.XYInfo.C00)) + fMaxDZ * abs(T1.tx))[i1_4] * iz,
- (sqrt(Pick_m22 * (T1.C11 + stam.XYInfo.C11)) + fMaxDZ * abs(T1.ty))[i1_4] * iz, time,
+ (sqrt(Pick_m22 * T1.C00) + fMaxDx[iStaM] + fMaxDZ * abs(T1.tx))[i1_4] * iz,
+ (sqrt(Pick_m22 * T1.C11) + fMaxDy[iStaM] + fMaxDZ * abs(T1.ty))[i1_4] * iz, time,
sqrt(timeError2) * 5);
for (L1HitIndex_t imh = -1; true;) { // loop over the hits in the area
@@ -447,12 +442,10 @@ inline void L1Algo::findDoubletsStep0(
fvec y, C11;
L1ExtrapolateYC11Line(T1, zm, y, C11);
- fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + stam.XYInfo.C11[i1_4]);
-
/// Covariation matrix of the hit
auto [dxxScalMhit, dxyScalMhit, dyyScalMhit] = stam.FormXYCovarianceMatrix(hitm.dU2(), hitm.dV2());
- if (fUseHitErrors) { dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + dyyScalMhit); }
+ fscal dy_est2 = Pick_m22[i1_4] * fabs(C11[i1_4] + dyyScalMhit);
auto [xm, ym] = stam.ConvUVtoXY<fscal>(hitm.U(), hitm.V());
@@ -464,9 +457,7 @@ inline void L1Algo::findDoubletsStep0(
fvec x, C00;
L1ExtrapolateXC00Line(T1, zm, x, C00);
- fscal dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + stam.XYInfo.C00[i1_4]);
-
- if (fUseHitErrors) { dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + dxxScalMhit); }
+ fscal dx_est2 = Pick_m22[i1_4] * fabs(C00[i1_4] + dxxScalMhit);
const fscal dX = xm - x[i1_4];
@@ -477,22 +468,13 @@ inline void L1Algo::findDoubletsStep0(
L1ExtrapolateC10Line(T1, zm, C10);
fvec chi2 = T1.chi2;
- L1UMeasurementInfo info = stam.frontInfo;
-
- if (fUseHitErrors) info.sigma2 = hitm.dU2();
-
- L1FilterChi2XYC00C10C11(info, x, y, C00, C10, C11, chi2, hitm.U());
+ L1FilterChi2XYC00C10C11(stam.frontInfo, x, y, C00, C10, C11, chi2, hitm.U(), hitm.dU2());
if (!fpCurrentIteration->GetTrackFromTripletsFlag()) {
if (chi2[i1_4] > fDoubletChi2Cut) continue;
}
-
- info = stam.backInfo;
-
- if (fUseHitErrors) info.sigma2 = hitm.dV2();
-
- L1FilterChi2(info, x, y, C00, C10, C11, chi2, hitm.V());
+ L1FilterChi2(stam.backInfo, x, y, C00, C10, C11, chi2, hitm.V(), hitm.dV2());
// FilterTime(T1, hitm.T(), hitm.dT2());
@@ -631,57 +613,22 @@ inline void L1Algo::findTripletsStep0( // input
// L1TrackPar tStore1 = T2;
- L1UMeasurementInfo info = stam.frontInfo;
-
- if (fUseHitErrors) info.sigma2 = du2_2;
-
// TODO: SG: L1FilterNoField is wrong.
// TODO: If the field was present before,
// TODO: the momentum is correlated with the position and corresponding
// TODO: matrix elements must be up[dated
- if (istam < fNfieldStations) { L1Filter(T2, info, u_front_2); }
- else {
- L1FilterNoField(T2, info, u_front_2);
- }
-
- // L1TrackPar tStore2 = T2;
-
- //assert(T2.IsConsistent(true, n2_4));
-
- // a good place to debug the fit
-
- /*
- if (tStore1.IsConsistent(false, n2_4) && !tStore2.IsConsistent(false, n2_4)) {
- cout << " i2 " << i2 << " dx2 " << dx2 << endl;
- cout << "\n\n before filtration: \n\n" << endl;
- tStore1.Print();
- tStore1.IsConsistent(true, n2_4);
- cout << "\n\n after filtration: \n\n" << endl;
- tStore2.Print();
- tStore2.IsConsistent(true, n2_4);
- cout << "\n\n " << n2_4 << " vector elements are filled." << endl;
- cout << "measurement: cos " << info.cos_phi[0] << " sin " << info.sin_phi[0] << " u " << u_front_2 << " s2 "
- << sqrt(info.sigma2) << endl;
- cout << "\n\n" << endl;
- exit(0);
+ if (istam < fNfieldStations) {
+ L1Filter(T2, stam.frontInfo, u_front_2, du2_2, fvec::One());
+ L1Filter(T2, stam.backInfo, u_back_2, dv2_2, fvec::One());
}
- */
-
- info = stam.backInfo;
- if (fUseHitErrors) info.sigma2 = dv2_2;
-
- if (istam < fNfieldStations) { L1Filter(T2, info, u_back_2); }
else {
- L1FilterNoField(T2, info, u_back_2);
+ L1FilterNoField(T2, stam.frontInfo, u_front_2, du2_2, fvec::One());
+ L1FilterNoField(T2, stam.backInfo, u_back_2, dv2_2, fvec::One());
}
- // assert(T2.IsConsistent(true, n2_4));
-
FilterTime(T2, t_2, dt2_2, stam.timeInfo);
- // assert(T2.IsConsistent(true, n2_4));
-
if (kMcbm == fTrackingMode) {
fit.L1AddThickMaterial(T2, fParameters.GetMaterialThickness(istam, T2.x, T2.y), fMaxInvMom, fvec::One(),
stam.fZthick, 1);
@@ -695,13 +642,9 @@ inline void L1Algo::findTripletsStep0( // input
//if ((istar >= fNstationsBeforePipe) && (istam <= fNstationsBeforePipe - 1)) { fit.L1AddPipeMaterial(T2, T2.qp, 1); }
- // assert(T2.IsConsistent(true, n2_4));
-
fvec dz2 = star.fZ - T2.z;
L1ExtrapolateTime(T2, dz2, stam.timeInfo);
- // assert(T2.IsConsistent(true, n2_4));
-
// extrapolate to the right hit station
if (istar <= fNfieldStations) {
@@ -731,8 +674,8 @@ inline void L1Algo::findTripletsStep0( // input
const fscal iz = 1.f / (T2.z[i2_4] - fParameters.GetTargetPositionZ()[0]);
L1HitAreaTime area(vGridTime[&star - fParameters.GetStations().begin()], T2.x[i2_4] * iz, T2.y[i2_4] * iz,
- (sqrt(Pick_r22 * (T2.C00 + stam.XYInfo.C00)) + fMaxDZ * abs(T2.tx))[i2_4] * iz,
- (sqrt(Pick_r22 * (T2.C11 + stam.XYInfo.C11)) + fMaxDZ * abs(T2.ty))[i2_4] * iz, time,
+ (sqrt(Pick_r22 * T2.C00) + fMaxDx[iStaM] + fMaxDZ * abs(T2.tx))[i2_4] * iz,
+ (sqrt(Pick_r22 * T2.C11) + fMaxDy[iStaM] + fMaxDZ * abs(T2.ty))[i2_4] * iz, time,
sqrt(timeError2) * 5);
L1HitIndex_t irh = 0;
@@ -781,17 +724,11 @@ inline void L1Algo::findTripletsStep0( // input
// check lower boundary
fvec y, C11;
L1ExtrapolateYC11Line(T2, zr, y, C11);
- // cout << "sta " << iStaR << " dy = " << sqrt(C11) << endl;
- fscal dy_est2 =
- (Pick_r22[i2_4]
- * (fabs(
- C11[i2_4]
- + star.XYInfo.C11[i2_4]))); // TODO for FastPrim dx < dy - other sort is optimal. But not for doublets
/// Covariation matrix of the hit
auto [dxxScalRhit, dxyScalRhit, dyyScalRhit] = star.FormXYCovarianceMatrix(hitr.dU2(), hitr.dV2());
- if (fUseHitErrors) { dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + dyyScalRhit))); }
+ fscal dy_est2 = (Pick_r22[i2_4] * (fabs(C11[i2_4] + dyyScalRhit)));
const fscal dY = yr - y[i2_4];
const fscal dY2 = dY * dY;
@@ -802,11 +739,7 @@ inline void L1Algo::findTripletsStep0( // input
L1ExtrapolateXC00Line(T2, zr, x, C00);
- // cout << "sta " << iStaR << " dx = " << sqrt(C00) << endl;
-
- fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + star.XYInfo.C00[i2_4])));
-
- if (fUseHitErrors) { dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + dxxScalRhit))); }
+ fscal dx_est2 = (Pick_r22[i2_4] * (fabs(C00[i2_4] + dxxScalRhit)));
const fscal dX = xr - x[i2_4];
if (dX * dX > dx_est2) continue;
@@ -815,16 +748,9 @@ inline void L1Algo::findTripletsStep0( // input
L1ExtrapolateC10Line(T2, zr, C10);
fvec chi2 = T2.chi2;
- info = star.frontInfo;
+ L1FilterChi2XYC00C10C11(star.frontInfo, x, y, C00, C10, C11, chi2, hitr.U(), hitr.dU2());
- if (fUseHitErrors) info.sigma2 = hitr.dU2();
-
- L1FilterChi2XYC00C10C11(info, x, y, C00, C10, C11, chi2, hitr.U());
- info = star.backInfo;
-
- if (fUseHitErrors) info.sigma2 = hitr.dV2();
-
- L1FilterChi2(info, x, y, C00, C10, C11, chi2, hitr.V());
+ L1FilterChi2(star.backInfo, x, y, C00, C10, C11, chi2, hitr.V(), hitr.dV2());
FilterTime(T_cur, hitr.T(), hitr.dT2(), star.timeInfo);
@@ -906,42 +832,22 @@ inline void L1Algo::findTripletsStep1( // input
L1TrackPar& T3 = T_3[i3_V];
- // assert(T3.IsConsistent(true, -1));
-
L1ExtrapolateTime(T3, dz, star.timeInfo);
L1ExtrapolateLine(T3, z_Pos[i3_V]);
- // assert(T3.IsConsistent(true, -1));
-
- L1UMeasurementInfo info = star.frontInfo;
-
- if (fUseHitErrors) info.sigma2 = du2_3[i3_V];
-
bool noField = (&star - fParameters.GetStations().begin() >= fNfieldStations);
- if (noField) { L1FilterNoField(T3, info, u_front_[i3_V]); }
- else {
- L1Filter(T3, info, u_front_[i3_V]);
+ if (noField) {
+ L1FilterNoField(T3, star.frontInfo, u_front_[i3_V], du2_3[i3_V], fvec::One());
+ L1FilterNoField(T3, star.backInfo, u_back_[i3_V], dv2_3[i3_V], fvec::One());
}
-
- // assert(T3.IsConsistent(true, -1));
-
- info = star.backInfo;
-
- if (fUseHitErrors) info.sigma2 = dv2_3[i3_V];
-
- if (noField) { L1FilterNoField(T3, info, u_back_[i3_V]); }
else {
- L1Filter(T3, info, u_back_[i3_V]);
+ L1Filter(T3, star.frontInfo, u_front_[i3_V], du2_3[i3_V], fvec::One());
+ L1Filter(T3, star.backInfo, u_back_[i3_V], dv2_3[i3_V], fvec::One());
}
- // assert(T3.IsConsistent(true, -1));
-
- if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) { FilterTime(T3, t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
-
- // assert(T3.IsConsistent(true, -1));
- // FilterTime(T_3[i3_V], timeR[i3_V], timeER[i3_V]);
+ if (kMcbm != fTrackingMode) { FilterTime(T3, t_3[i3_V], dt2_3[i3_V], star.timeInfo); }
}
}
@@ -1746,11 +1652,21 @@ void L1Algo::CATrackFinder()
vHitsUnused = &vNotUsedHits_Buf;
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
+ const L1Station& st = fParameters.GetStation(iS);
+ fMaxDx[iS] = 0.;
+ fMaxDy[iS] = 0.;
bool timeUninitialised = 1;
fscal lasttime = 0;
fscal starttime = 0;
for (L1HitIndex_t ih = fInputData.GetStartHitIndex(iS); ih < fInputData.GetStopHitIndex(iS); ++ih) {
- const fscal time = fInputData.GetHit(ih).t;
+ const L1Hit& h = fInputData.GetHit(ih);
+ auto [dxx, dxy, dyy] = st.FormXYCovarianceMatrix(h.du2, h.dv2);
+ fscal dx = sqrt(dxx);
+ fscal dy = sqrt(dyy);
+ if (fMaxDx[iS] < dx) { fMaxDx[iS] = dx; }
+ if (fMaxDy[iS] < dy) { fMaxDy[iS] = dy; }
+
+ const fscal time = h.t;
assert(std::isfinite(time));
if (timeUninitialised || lasttime < time) { lasttime = time; }
if (timeUninitialised || starttime > time) { starttime = time; }
diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h
index 9797c173d8514307412d965d3b34ddb29f22b2f2..cf74fcbff6944868cf7c5c25e8c1827a011b87fa 100644
--- a/reco/L1/L1Algo/L1Filtration.h
+++ b/reco/L1/L1Algo/L1Filtration.h
@@ -83,7 +83,7 @@ inline void FilterTime(L1TrackPar& T, fvec t, fvec dt2, fvec timeInfo = fvec::On
}
-inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec w = 1.)
+inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec du2, fvec w)
{
// filter the track T with an 1-D measurement u
@@ -99,19 +99,19 @@ inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec
fvec HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
- // when the measurement error info.sigma2 is much smaller than the current track error HCH,
+ // when the measurement error du2 is much smaller than the current track error HCH,
// move track exactly to the measurement without filtering
// it helps to keep the initial track errors reasonably small
// the calculations in the covariance matrix are not affected
- const fmask maskDoFilter = (HCH < info.sigma2 * 16.f);
+ const fmask maskDoFilter = (HCH < du2 * 16.f);
//const fvec maskDoFilter = _f32vec4_true;
- // correction to HCH is needed for the case when sigma2 is so small
+ // correction to HCH is needed for the case when du2 is so small
// with respect to HCH that it disappears due to the roundoff error
//
- fvec wi = w / (info.sigma2 + 1.0000001f * HCH);
- fvec zetawi = w * zeta / (iif(maskDoFilter, info.sigma2, fvec::Zero()) + HCH);
+ fvec wi = w / (du2 + 1.0000001f * HCH);
+ fvec zetawi = w * zeta / (iif(maskDoFilter, du2, fvec::Zero()) + HCH);
// T.chi2 += iif( maskDoFilter, zeta * zetawi, fvec::Zero() );
T.chi2 += zeta * zeta * wi;
@@ -154,7 +154,7 @@ inline void L1Filter(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec
}
-inline void L1FilterNoField(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec w = 1.)
+inline void L1FilterNoField(L1TrackPar& T, const L1UMeasurementInfo& info, fvec u, fvec du2, fvec w)
{
fvec zeta, HCH;
fvec F0, F1, F2, F3, F4, F5;
@@ -173,14 +173,14 @@ inline void L1FilterNoField(L1TrackPar& T, const L1UMeasurementInfo& info, fvec
F4 = info.cos_phi * T.C40 + info.sin_phi * T.C41;
F5 = info.cos_phi * T.C50 + info.sin_phi * T.C51;
- //const fmask maskDoFilter = (HCH < info.sigma2 * 16.f);
+ //const fmask maskDoFilter = (HCH < du2 * 16.f);
const fmask maskDoFilter(fmask::One());
//TODO: SG: try this
- //fvec wi = w / (info.sigma2 + 1.0000001f * HCH);
+ //fvec wi = w / (du2 + 1.0000001f * HCH);
- fvec wi = w / (info.sigma2 + HCH);
- fvec zetawi = w * zeta / (iif(maskDoFilter, info.sigma2, fvec::Zero()) + HCH);
+ fvec wi = w / (du2 + HCH);
+ fvec zetawi = w * zeta / (iif(maskDoFilter, du2, fvec::Zero()) + HCH);
T.chi2 += zeta * zeta * wi;
T.NDF += w;
@@ -223,7 +223,7 @@ inline void L1FilterNoField(L1TrackPar& T, const L1UMeasurementInfo& info, fvec
}
inline void L1FilterChi2(const L1UMeasurementInfo& info, const fvec& x, const fvec& y, const fvec& C00, const fvec& C10,
- const fvec& C11, fvec& chi2, const fvec& u)
+ const fvec& C11, fvec& chi2, const fvec& u, const fvec& du2)
{
fvec zeta, HCH;
fvec F0, F1;
@@ -236,11 +236,11 @@ inline void L1FilterChi2(const L1UMeasurementInfo& info, const fvec& x, const fv
HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
- chi2 += zeta * zeta / (info.sigma2 + HCH);
+ chi2 += zeta * zeta / (du2 + HCH);
}
inline void L1FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fvec& y, fvec& C00, fvec& C10, fvec& C11,
- fvec& chi2, const fvec& u)
+ fvec& chi2, const fvec& u, const fvec& du2)
{
fvec wi, zeta, zetawi, HCH;
fvec F0, F1;
@@ -254,7 +254,7 @@ inline void L1FilterChi2XYC00C10C11(const L1UMeasurementInfo& info, fvec& x, fve
HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
- wi = fvec(1.) / (info.sigma2 + HCH);
+ wi = fvec(1.) / (du2 + HCH);
zetawi = zeta * wi;
chi2 += zeta * zetawi;
@@ -417,146 +417,6 @@ inline void L1FilterXY(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo&
T.C44 -= K40 * F40 + K41 * F41;
}
-/*
-inline void L1Filter1D( L1TrackPar &T, fvec &u, fvec &sigma2, fvec H[] )
-{
-
- cnst ZERO = 0.0, ONE = 1.;
-
- fvec wi, zeta, zetawi, HCH;
- fvec F0, F1, F2, F3, F4;
- fvec K1, K2, K3, K4;
-
- zeta = H[0]*T.x +H[1]*T.y +H[2]*T.tx +H[3]*T.ty +H[4]*T.qp - u;
-
- // F = CH'
- F0 = H[0]*T.C00 + H[1]*T.C10 + H[2]*T.C20+ H[3]*T.C30+ H[4]*T.C40;
- F1 = H[0]*T.C10 + H[1]*T.C11 + H[2]*T.C21+ H[3]*T.C31+ H[4]*T.C41;
- F2 = H[0]*T.C20 + H[1]*T.C21 + H[2]*T.C22+ H[3]*T.C32+ H[4]*T.C42;
- F3 = H[0]*T.C30 + H[1]*T.C31 + H[2]*T.C32+ H[3]*T.C33+ H[4]*T.C43;
- F4 = H[0]*T.C40 + H[1]*T.C41 + H[2]*T.C42+ H[3]*T.C43+ H[4]*T.C44;
-
- HCH = H[0]*F0 +H[1]*F1 +H[2]*F2 +H[3]*F3 +H[4]*F4 ;
-
- wi = rcp(fvec(sigma2 +HCH));
- zetawi = zeta * wi;
- T.chi2 += zeta * zetawi;
- T.NDF += ONE;
-
- K1 = F1*wi;
- K2 = F2*wi;
- K3 = F3*wi;
- K4 = F4*wi;
-
- T.x -= F0*zetawi;
- T.y -= F1*zetawi;
- T.tx -= F2*zetawi;
- T.ty -= F3*zetawi;
- T.qp -= F4*zetawi;
-
- T.C00-= F0*F0*wi;
- T.C10-= K1*F0;
- T.C11-= K1*F1;
- T.C20-= K2*F0;
- T.C21-= K2*F1;
- T.C22-= K2*F2;
- T.C30-= K3*F0;
- T.C31-= K3*F1;
- T.C32-= K3*F2;
- T.C33-= K3*F3;
- T.C40-= K4*F0;
- T.C41-= K4*F1;
- T.C42-= K4*F2;
- T.C43-= K4*F3;
- T.C44-= K4*F4;
-}
-
-
-inline void L1Filter2D( L1TrackPar &T, fvec &x, fvec &y, L1XYMeasurementInfo &info, fvec H[] )
-{
- cnst TWO = 2.;
-
- fvec zeta0, zeta1, S00, S10, S11, si;
- fvec F00, F10, F20, F30, F40, F01, F11, F21, F31, F41 ;
- fvec K00, K10, K20, K30, K40, K01, K11, K21, K31, K41;
-
- zeta0 = H[0]*T.x +H[1]*T.y +H[2]*T.tx +H[3]*T.ty +H[4]*T.qp - x;
- zeta1 = H[5]*T.x +H[6]*T.y +H[7]*T.tx +H[8]*T.ty +H[9]*T.qp - y;
-
- // F = CH'
- F00 = H[0]*T.C00 + H[1]*T.C10 + H[2]*T.C20+ H[3]*T.C30+ H[4]*T.C40;
- F10 = H[0]*T.C10 + H[1]*T.C11 + H[2]*T.C21+ H[3]*T.C31+ H[4]*T.C41;
- F20 = H[0]*T.C20 + H[1]*T.C21 + H[2]*T.C22+ H[3]*T.C32+ H[4]*T.C42;
- F30 = H[0]*T.C30 + H[1]*T.C31 + H[2]*T.C32+ H[3]*T.C33+ H[4]*T.C43;
- F40 = H[0]*T.C40 + H[1]*T.C41 + H[2]*T.C42+ H[3]*T.C43+ H[4]*T.C44;
- F01 = H[5]*T.C00 + H[6]*T.C10 + H[7]*T.C20+ H[8]*T.C30+ H[9]*T.C40;
- F11 = H[5]*T.C10 + H[6]*T.C11 + H[7]*T.C21+ H[8]*T.C31+ H[9]*T.C41;
- F21 = H[5]*T.C20 + H[6]*T.C21 + H[7]*T.C22+ H[8]*T.C32+ H[9]*T.C42;
- F31 = H[5]*T.C30 + H[6]*T.C31 + H[7]*T.C32+ H[8]*T.C33+ H[9]*T.C43;
- F41 = H[5]*T.C40 + H[6]*T.C41 + H[7]*T.C42+ H[8]*T.C43+ H[9]*T.C44;
-
- S00 = H[0]*F00 +H[1]*F10 +H[2]*F20 +H[3]*F30 +H[4]*F40 + info.C00;
- S10 = H[5]*F00 +H[6]*F10 +H[7]*F20 +H[8]*F30 +H[9]*F40 + info.C10;
- S11 = H[5]*F01 +H[6]*F11 +H[7]*F21 +H[8]*F31 +H[9]*F41 + info.C11;
-
- si = rcp(fvec(S00*S11 - S10*S10));
- fvec S00tmp = S00;
- S00 = si*S11;
- S10 = -si*S10;
- S11 = si*S00tmp;
-
- T.chi2+= zeta0*zeta0*S00 + 2.*zeta0*zeta1*S10 + zeta1*zeta1*S11;
- T.NDF += TWO;
-
- 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;
-
- T.x -= K00*zeta0 + K01*zeta1;
- T.y -= K10*zeta0 + K11*zeta1;
- T.tx -= K20*zeta0 + K21*zeta1;
- T.ty -= K30*zeta0 + K31*zeta1;
- T.qp -= K40*zeta0 + K41*zeta1;
-
- T.C00-= K00*F00 + K01*F01;
- T.C10-= K10*F00 + K11*F01;
- T.C11-= K10*F10 + K11*F11;
- T.C20-= K20*F00 + K21*F01;
- T.C21-= K20*F10 + K21*F11;
- T.C22-= K20*F20 + K21*F21;
- T.C30-= K30*F00 + K31*F01;
- T.C31-= K30*F10 + K31*F11;
- T.C32-= K30*F20 + K31*F21;
- T.C33-= K30*F30 + K31*F31;
- T.C40-= K40*F00 + K41*F01;
- T.C41-= K40*F10 + K41*F11;
- T.C42-= K40*F20 + K41*F21;
- T.C43-= K40*F30 + K41*F31;
- T.C44-= K40*F40 + K41*F41;
-}
-
-
-
-
-inline void L1FilterFirst( fvec *T, L1Cov &C, fvec &Chi2, fvec &NDF, L1XYMeasurementInfo &info, fvec &x, fvec &y, fvec &w )
-{
- cnst ZERO = 0.0, ONE = 1.;
-
- fvec w1 = !w;
- C.C00= w&info.C00 + w1&INF;
- C.C10= w&info.C10 + w1&INF; C.C11= w&info.C11 + w1&INF;
- C.C20= ZERO; C.C21= ZERO; C.C22= INF;
- C.C30= ZERO; C.C31= ZERO; C.C32= ZERO; C.C33= INF;
- C.C40= ZERO; C.C41= ZERO; C.C42= ZERO; C.C43= ZERO; C.C44= INF;
-
- T[0] = w&x + w1&T[0];
- T[1] = w&y + w1&T[1];
- NDF = -3.0;
- Chi2 = ZERO;
-}
-*/
#undef cnst
diff --git a/reco/L1/L1Algo/L1Station.cxx b/reco/L1/L1Algo/L1Station.cxx
index d8f46fa761bc0d84cff51a6df9b3b9f1cb6f3768..8e32bc91551f61f4382a513cbb72001acbe40240 100644
--- a/reco/L1/L1Algo/L1Station.cxx
+++ b/reco/L1/L1Algo/L1Station.cxx
@@ -82,7 +82,6 @@ void L1Station::CheckConsistency() const
backInfo.CheckConsistency();
xInfo.CheckConsistency();
yInfo.CheckConsistency();
- XYInfo.CheckConsistency();
// Check consistency of coordinate transformations
for (float x = -1.f; x < 1.1f; x += 0.2) {
@@ -134,8 +133,6 @@ std::string L1Station::ToString(int verbosityLevel, int indentLevel) const
aStream << frontInfo.ToString(indentLevel + 2) << '\n';
aStream << indent << indentChar << "Back:\n";
aStream << backInfo.ToString(indentLevel + 2) << '\n';
- aStream << indent << indentChar << "XY cov matrix:\n";
- aStream << XYInfo.ToString(indentLevel + 2) << '\n';
aStream << indent << indentChar << "X layer:\n";
aStream << xInfo.ToString(indentLevel + 2) << '\n';
aStream << indent << indentChar << "Y layer:\n";
diff --git a/reco/L1/L1Algo/L1Station.h b/reco/L1/L1Algo/L1Station.h
index 035d99974c403fafee133fb60e2f33b7f2b35978..ddc3cb76e4d05eb621e951bd4722cf0880b59a3e 100644
--- a/reco/L1/L1Algo/L1Station.h
+++ b/reco/L1/L1Algo/L1Station.h
@@ -35,7 +35,6 @@ public:
L1UMeasurementInfo backInfo {};
L1UMeasurementInfo xInfo {}; ///< x axis in front,back coordinates
L1UMeasurementInfo yInfo {}; ///< y axis in front,back coordinates
- L1XYMeasurementInfo XYInfo {};
/// Serialization function
friend class boost::serialization::access;
@@ -57,7 +56,6 @@ public:
ar& backInfo;
ar& xInfo;
ar& yInfo;
- ar& XYInfo;
}
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index 79028053cbb8ccc61f8c59c93f1cf69bd0511362..77329ba97af85b59994c93cbbce21126aa4dc6b4 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -95,9 +95,8 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
T.z = z0;
T.chi2 = fvec(0.);
T.NDF = fvec(2.);
- T.C00 = sta0.XYInfo.C00;
- T.C10 = sta0.XYInfo.C10;
- T.C11 = sta0.XYInfo.C11;
+
+ std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2);
T.C20 = T.C21 = fvec(0.);
T.C30 = T.C31 = T.C32 = fvec(0.);
@@ -139,8 +138,8 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
fvec u = hit.u;
fvec v = hit.v;
auto [x, y] = sta.ConvUVtoXY<fvec>(u, v);
- L1Filter(T, sta.frontInfo, u);
- L1Filter(T, sta.backInfo, v);
+ L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
+ L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
fldB0 = fldB1;
fldB1 = fldB2;
fldZ0 = fldZ1;
@@ -222,11 +221,11 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
// T.tx = (x1-x0)*dzi;
// T.ty = (y1-y0)*dzi;
// qp0 = 0;
- T.qp = qp0;
- T.z = z0;
- T.C00 = sta0.XYInfo.C00;
- T.C10 = sta0.XYInfo.C10;
- T.C11 = sta0.XYInfo.C11;
+ T.qp = qp0;
+ T.z = z0;
+
+ std::tie(T.C00, T.C10, T.C11) = sta0.FormXYCovarianceMatrix(hit0.du2, hit0.dv2);
+
T.C20 = T.C21 = 0;
T.C30 = T.C31 = T.C32 = 0;
T.C40 = T.C41 = T.C42 = T.C43 = 0;
@@ -263,8 +262,8 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
fit.L1AddMaterial(T, fParameters.GetMaterialThickness(i, T.x, T.y), qp0, ONE);
// if (ista == fNstationsBeforePipe) fit.L1AddPipeMaterial( T, qp0, 1);
- L1Filter(T, sta.frontInfo, u);
- L1Filter(T, sta.backInfo, v);
+ L1Filter(T, sta.frontInfo, u, hit.du2, fvec::One());
+ L1Filter(T, sta.backInfo, v, hit.dv2, fvec::One());
fldB0 = fldB1;
fldB1 = fldB2;
@@ -332,7 +331,6 @@ void L1Algo::L1KFTrackFitter()
L1Track* t[fvec::size()] {nullptr};
const L1Station* sta = fParameters.GetStations().begin();
- L1Station staFirst, staLast; // FIXME (?): Probably, we should replace these variables with references (S.Zharko)
// Spatial-time position of a hit vs. station and track in the portion
// NOTE: u- and v-axes are axes, orthogonal to front and back strips of the station, respectively.
@@ -447,30 +445,24 @@ 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];
- staFirst.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
- staFirst.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
- staFirst.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[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];
+ 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];
- staLast.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
- staLast.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
- staLast.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[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];
}
}
@@ -506,7 +498,7 @@ void L1Algo::L1KFTrackFitter()
time_last = iif(w_time[ista] > fvec::Zero(), time_last, fvec::Zero());
dt2_last = iif(w_time[ista] > fvec::Zero(), dt2_last, fvec(100.));
- FilterFirst(fit, x_last, y_last, time_last, dt2_last, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
+ FilterFirst(fit, x_last, y_last, time_last, dt2_last, d_xx_lst, d_yy_lst, d_xy_lst);
fldZ1 = z[ista];
@@ -561,12 +553,10 @@ void L1Algo::L1KFTrackFitter()
L1UMeasurementInfo vtxInfoX;
vtxInfoX.cos_phi = 1.;
vtxInfoX.sin_phi = 0.;
- vtxInfoX.sigma2 = fvec(1.e-4);
L1UMeasurementInfo vtxInfoY;
vtxInfoY.cos_phi = 0;
vtxInfoY.sin_phi = 1.;
- vtxInfoY.sigma2 = fvec(1.e-4);
if (kGlobal == fTrackingMode) {
fvec qp00 = tr.qp;
@@ -661,7 +651,7 @@ void L1Algo::L1KFTrackFitter()
ista = 0;
- FilterFirst(fit, x_first, y_first, time_first, dt2_first, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
+ FilterFirst(fit, x_first, y_first, time_first, dt2_first, d_xx_fst, d_yy_fst, d_xy_fst);
qp01 = tr.qp;
@@ -896,75 +886,12 @@ void L1Algo::GuessVec(L1TrackParFit& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy,
t.fTr.z = z0;
}
-void L1Algo::FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st)
-{
- // initialize covariance matrix
- track.C00 = st.XYInfo.C00;
- track.C10 = st.XYInfo.C10;
- track.C11 = st.XYInfo.C11;
- track.C20 = ZERO;
- track.C21 = ZERO;
- track.C22 = vINF;
- track.C30 = ZERO;
- track.C31 = ZERO;
- track.C32 = ZERO;
- track.C33 = vINF;
- track.C40 = ZERO;
- track.C41 = ZERO;
- track.C42 = ZERO;
- track.C43 = ZERO;
- track.C44 = ONE;
- track.C50 = ZERO;
- track.C51 = ZERO;
- track.C52 = ZERO;
- track.C53 = ZERO;
- track.C55 = 2.6f * 2.6f;
-
- track.x = x;
- track.y = y;
- track.NDF = -3.0;
- track.chi2 = ZERO;
-}
-
-void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, L1Station& st)
-{
- // initialize covariance matrix
- track.fTr.C00 = st.XYInfo.C00;
- track.fTr.C10 = st.XYInfo.C10;
- track.fTr.C11 = st.XYInfo.C11;
- track.fTr.C20 = ZERO;
- track.fTr.C21 = ZERO;
- track.fTr.C22 = vINF;
- track.fTr.C30 = ZERO;
- track.fTr.C31 = ZERO;
- track.fTr.C32 = ZERO;
- track.fTr.C33 = vINF;
- track.fTr.C40 = ZERO;
- track.fTr.C41 = ZERO;
- track.fTr.C42 = ZERO;
- track.fTr.C43 = ZERO;
- track.fTr.C44 = ONE;
- track.fTr.C50 = ZERO;
- track.fTr.C51 = ZERO;
- track.fTr.C52 = ZERO;
- track.fTr.C53 = ZERO;
- track.fTr.C54 = ZERO;
- track.fTr.C55 = 2.6f * 2.6f;
-
- track.fTr.x = x;
- track.fTr.y = y;
- track.fTr.t = t;
- track.fTr.NDF = -3.0;
- track.fTr.chi2 = ZERO;
-}
-
-void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, L1Station& st, fvec& /*d_xx*/,
- fvec& /*d_yy*/, fvec& /*d_xy*/)
+void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec& dt2, fvec& d_xx, fvec& d_yy, fvec& d_xy)
{
- track.fTr.C00 = st.XYInfo.C00;
- track.fTr.C10 = st.XYInfo.C10;
- track.fTr.C11 = st.XYInfo.C11;
+ track.fTr.C00 = d_xx;
+ track.fTr.C10 = d_xy;
+ track.fTr.C11 = d_yy;
track.fTr.C20 = ZERO;
track.fTr.C21 = ZERO;
track.fTr.C22 = vINF;
@@ -990,70 +917,3 @@ void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& t, fvec&
track.fTr.NDF = -3.0;
track.fTr.chi2 = ZERO;
}
-
-
-void L1Algo::FilterFirst(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt2, L1Station& st)
-{
- track.fTr.C00 = st.XYInfo.C00;
- track.fTr.C10 = st.XYInfo.C10;
- track.fTr.C11 = st.XYInfo.C11;
- track.fTr.C20 = ZERO;
- track.fTr.C21 = ZERO;
- track.fTr.C22 = vINF;
- track.fTr.C30 = ZERO;
- track.fTr.C31 = ZERO;
- track.fTr.C32 = ZERO;
- track.fTr.C33 = vINF;
- track.fTr.C40 = ZERO;
- track.fTr.C41 = ZERO;
- track.fTr.C42 = ZERO;
- track.fTr.C43 = ZERO;
- track.fTr.C44 = ONE;
- track.fTr.C50 = ZERO;
- track.fTr.C51 = ZERO;
- track.fTr.C52 = ZERO;
- track.fTr.C53 = ZERO;
- track.fTr.C54 = ZERO;
- track.fTr.C55 = dt2;
-
- track.fTr.x = x;
- track.fTr.y = y;
- track.fTr.NDF = -3.0;
- track.fTr.chi2 = ZERO;
-}
-
-void L1Algo::FilterFirstL(L1TrackParFit& track, fvec& x, fvec& y, fvec& /*t*/, fvec& dt2, L1Station& /*st*/, fvec& d_xx,
- fvec& d_yy, fvec& d_xy)
-{
- // initialize covariance matrix
- track.fTr.C00 = d_xx;
- track.fTr.C10 = d_xy;
- track.fTr.C11 = d_yy;
- // track.fTr.C00= st.XYInfo.C00;
- // track.fTr.C10= st.XYInfo.C10;
- // track.fTr.C11= st.XYInfo.C11;
- track.fTr.C20 = ZERO;
- track.fTr.C21 = ZERO;
- track.fTr.C22 = vINF;
- track.fTr.C30 = ZERO;
- track.fTr.C31 = ZERO;
- track.fTr.C32 = ZERO;
- track.fTr.C33 = vINF;
- track.fTr.C40 = ZERO;
- track.fTr.C41 = ZERO;
- track.fTr.C42 = ZERO;
- track.fTr.C43 = ZERO;
- track.fTr.C44 = ONE;
- track.fTr.C50 = ZERO;
- track.fTr.C51 = ZERO;
- track.fTr.C52 = ZERO;
- track.fTr.C53 = ZERO;
- track.fTr.C54 = ZERO;
- track.fTr.C55 = dt2;
-
- track.fTr.x = x;
- track.fTr.y = y;
- // track.fTr.t = t;
- track.fTr.NDF = -3.0; // TODO: Why -3.0? (S.Zharko)
- track.fTr.chi2 = ZERO;
-}
diff --git a/reco/L1/L1Algo/L1TrackParFit.cxx b/reco/L1/L1Algo/L1TrackParFit.cxx
index d19f8c014908471be39ee61226f92a0416a2a0ee..be04031d16afd6620e7999c105cd67fb9faa1b6f 100644
--- a/reco/L1/L1Algo/L1TrackParFit.cxx
+++ b/reco/L1/L1Algo/L1TrackParFit.cxx
@@ -77,7 +77,7 @@ void L1TrackParFit::Filter(const L1UMeasurementInfo& info, const fvec& u, const
fTr.C55 -= K5 * F5;
}
-void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w)
+void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec du2, fvec w)
{
fvec wi, zeta, zetawi, HCH;
fvec F0, F1, F2, F3, F4, F5;
@@ -97,12 +97,12 @@ void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w)
F5 = info.cos_phi * fTr.C50 + info.sin_phi * fTr.C51;
#if 0 // use mask
- const fmask mask = (HCH < info.sigma2 * 16.);
- wi = w/( (mask & info.sigma2) +HCH );
+ const fmask mask = (HCH < du2 * 16.);
+ wi = w/( (mask & du2) +HCH );
zetawi = zeta *wi;
fTr.chi2 += mask & (zeta * zetawi);
#else
- wi = w / (info.sigma2 + HCH);
+ wi = w / (du2 + HCH);
zetawi = zeta * wi;
fTr.chi2 += zeta * zetawi;
#endif // 0
diff --git a/reco/L1/L1Algo/L1TrackParFit.h b/reco/L1/L1Algo/L1TrackParFit.h
index 20caaa3290c7aed5867ff0b84a6606c1f562a962..e304521b13b07e8b2da87f4687845defc6f7487b 100644
--- a/reco/L1/L1Algo/L1TrackParFit.h
+++ b/reco/L1/L1Algo/L1TrackParFit.h
@@ -49,7 +49,7 @@ public:
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 t, fvec dt2, fvec w, fvec timeInfo);
- void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w = 1.);
+ void FilterNoP(L1UMeasurementInfo& info, fvec u, fvec du2, fvec w);
void Extrapolate(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);
void ExtrapolateStep(fvec z_out, fvec qp0, const L1FieldRegion& F, const fvec& w);
diff --git a/reco/L1/L1Algo/L1UMeasurementInfo.cxx b/reco/L1/L1Algo/L1UMeasurementInfo.cxx
index 3d86cffba4eca85e7fbdbdb263fb6ad177a9af0e..a8e45ad248f0bf402d6467422ca487d49a934885 100644
--- a/reco/L1/L1Algo/L1UMeasurementInfo.cxx
+++ b/reco/L1/L1Algo/L1UMeasurementInfo.cxx
@@ -17,7 +17,6 @@ std::string L1UMeasurementInfo::ToString(int indentLevel) const
std::string indent(indentLevel, indentChar);
aStream << indent << "cos(phi): " << std::setw(12) << std::setfill(' ') << cos_phi[0] << '\n';
aStream << indent << "sin(phi): " << std::setw(12) << std::setfill(' ') << sin_phi[0] << '\n';
- aStream << indent << "sigma2 [cm2]: " << std::setw(12) << std::setfill(' ') << sigma2[0];
return aStream.str();
}
@@ -25,15 +24,14 @@ std::string L1UMeasurementInfo::ToString(int indentLevel) const
//
void L1UMeasurementInfo::CheckConsistency() const
{
- /* (i) Checks for the horizontal equality of SIMD vector elements */
+ // (i) Checks for the horizontal equality of SIMD vector elements
L1Utils::CheckSimdVectorEquality(cos_phi, "L1UMeasurementsInfo::cos_phi");
L1Utils::CheckSimdVectorEquality(sin_phi, "L1UMeasurementsInfo::sin_phi");
- L1Utils::CheckSimdVectorEquality(sigma2, "L1UMeasurementsInfo::sigma2");
-
- /*(ii) Sigma2 possible values*/
- if (sigma2[0] < 0) {
+ /*
+ if (fabs( cos_phi[0]*cos_phi[0]+sin_phi[0]*sin_phi[0]-1.) >1.e-8) {
std::stringstream msg;
- msg << "L1UMeasurementsInfo: illegal value of sigma2: " << sigma2[0] << " [cm2] (sigma2 >= 0 excepted)";
+ msg << "L1UMeasurementsInfo: cos_phi and sin_phi do not match: " << cos_phi[0] << " "<<sin_phi[0];
throw std::logic_error(msg.str());
}
+ */
}
diff --git a/reco/L1/L1Algo/L1UMeasurementInfo.h b/reco/L1/L1Algo/L1UMeasurementInfo.h
index 63a82ac90fabaf9c23f753c1f42c027fa33053c8..5405719c6036efe7ca197fc98a1a125d761658f4 100644
--- a/reco/L1/L1Algo/L1UMeasurementInfo.h
+++ b/reco/L1/L1Algo/L1UMeasurementInfo.h
@@ -17,7 +17,6 @@ class L1UMeasurementInfo {
public:
fvec cos_phi {L1NaN::SetNaN<decltype(cos_phi)>()};
fvec sin_phi {L1NaN::SetNaN<decltype(sin_phi)>()};
- fvec sigma2 {L1NaN::SetNaN<decltype(sigma2)>()};
/// String representation of class contents
/// \param indentLevel number of indent characters in the output
@@ -27,7 +26,7 @@ public:
void CheckConsistency() const;
/// Checks, if the fields are NaN
- bool IsNaN() const { return L1NaN::IsNaN(cos_phi) || L1NaN::IsNaN(sin_phi) || L1NaN::IsNaN(sigma2); }
+ bool IsNaN() const { return L1NaN::IsNaN(cos_phi) || L1NaN::IsNaN(sin_phi); }
/// Serialization function
friend class boost::serialization::access;
@@ -36,7 +35,6 @@ public:
{
ar& cos_phi;
ar& sin_phi;
- ar& sigma2;
}
} _fvecalignment;
diff --git a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
index c55858e65d4c814ffb45e7dcf376c0a816def1eb..c6aefbad11454f271c7442e9231b1c46e837c43d 100644
--- a/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
+++ b/reco/L1/ParticleFinder/CbmL1PFFitter.cxx
@@ -98,11 +98,11 @@ inline void CbmL1PFFitter::PFFieldRegion::getL1FieldRegion(L1FieldRegion& fld, i
inline CbmL1PFFitter::PFFieldRegion::PFFieldRegion(const L1FieldRegion& fld, int i) { setFromL1FieldRegion(fld, i); }
-void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, L1Station& st)
+void FilterFirst(L1TrackPar& track, fvec& x, fvec& y, fvec& dxx, fvec& dxy, fvec& dyy)
{
- track.C00 = st.XYInfo.C00;
- track.C10 = st.XYInfo.C10;
- track.C11 = st.XYInfo.C11;
+ track.C00 = dxx;
+ track.C10 = dxy;
+ track.C11 = dyy;
track.C20 = ZERO;
track.C21 = ZERO;
track.C22 = vINF;
@@ -147,15 +147,20 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
int ista;
const L1Station* sta = CbmL1::Instance()->fpAlgo->GetParameters()->GetStations().begin();
- L1Station staFirst, staLast;
- fvec* x = new fvec[nHits];
- fvec* u = new fvec[nHits];
- fvec* v = new fvec[nHits];
- fvec* y = new fvec[nHits];
- fvec* z = new fvec[nHits];
- fvec* w = new fvec[nHits];
+
+ fvec* x = new fvec[nHits];
+ fvec* u = new fvec[nHits];
+ fvec* v = new fvec[nHits];
+ fvec* du2 = new fvec[nHits];
+ fvec* dv2 = new fvec[nHits];
+ fvec* y = new fvec[nHits];
+ fvec* z = new fvec[nHits];
+ fvec* w = new fvec[nHits];
// fvec y_temp;
fvec x_first, y_first, x_last, y_last;
+ fvec fstC00, fstC10, fstC11;
+ fvec lstC00, lstC10, lstC11;
+
fvec z0, z1, z2, dz, z_start, z_end;
L1FieldValue* fB = new L1FieldValue[nHits];
L1FieldValue fB_temp _fvecalignment;
@@ -228,6 +233,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
ista =
CbmL1::Instance()->fpAlgo->GetParameters()->GetStationIndexActive(hit->GetStationNr(), L1DetectorID::kMvd);
if (ista == -1) continue;
+ du2[ista][iVec] = hit->GetDx() * hit->GetDx();
+ dv2[ista][iVec] = hit->GetDy() * hit->GetDy();
}
else {
if (!listStsHits) continue;
@@ -242,6 +249,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
ista = CbmL1::Instance()->fpAlgo->GetParameters()->GetStationIndexActive(
CbmStsSetup::Instance()->GetStationNumber(hit->GetAddress()), L1DetectorID::kSts);
if (ista == -1) continue;
+ du2[ista][iVec] = hit->GetDu() * hit->GetDu();
+ dv2[ista][iVec] = hit->GetDv() * hit->GetDv();
}
w[ista][iVec] = 1.f;
@@ -250,31 +259,40 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
u[ista][iVec] = posx * sta[ista].frontInfo.cos_phi[0] + posy * sta[ista].frontInfo.sin_phi[0];
v[ista][iVec] = posx * sta[ista].backInfo.cos_phi[0] + posy * sta[ista].backInfo.sin_phi[0];
z[ista][iVec] = posz;
+
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
fB[ista].z[iVec] = fB_temp.z[iVec];
if (i == 0) {
- z_start[iVec] = posz;
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
- staFirst.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
- staFirst.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
- staFirst.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec];
+ z_start[iVec] = posz;
+ x_first[iVec] = x[ista][iVec];
+ y_first[iVec] = y[ista][iVec];
+
+ auto [dxx, dxy, dyy] = sta[ista].FormXYCovarianceMatrix((fscal) du2[ista][iVec], (fscal) dv2[ista][iVec]);
+
+ fstC00[iVec] = dxx;
+ fstC10[iVec] = dxy;
+ fstC11[iVec] = dyy;
}
if (i == nHitsTrack - 1) {
- z_end[iVec] = posz;
- x_last[iVec] = x[ista][iVec];
- y_last[iVec] = y[ista][iVec];
- staLast.XYInfo.C00[iVec] = sta[ista].XYInfo.C00[iVec];
- staLast.XYInfo.C10[iVec] = sta[ista].XYInfo.C10[iVec];
- staLast.XYInfo.C11[iVec] = sta[ista].XYInfo.C11[iVec];
+ z_end[iVec] = posz;
+ x_last[iVec] = x[ista][iVec];
+ y_last[iVec] = y[ista][iVec];
+
+ auto [dxx, dxy, dyy] = sta[ista].FormXYCovarianceMatrix((fscal) du2[ista][iVec], (fscal) dv2[ista][iVec]);
+
+ lstC00[iVec] = dxx;
+ lstC10[iVec] = dxy;
+ lstC11[iVec] = dyy;
}
}
}
// fit forward
i = 0;
- FilterFirst(T, x_first, y_first, staFirst);
+
+ FilterFirst(T, x_first, y_first, fstC00, fstC10, fstC11);
+
fvec qp0 = T.qp;
z1 = z[i];
sta[i].fieldSlice.GetFieldValue(T.x, T.y, b1);
@@ -303,8 +321,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0,
-1, wIn);
- L1Filter(T, sta[i].frontInfo, u[i], w1);
- L1Filter(T, sta[i].backInfo, v[i], w1);
+ L1Filter(T, sta[i].frontInfo, u[i], du2[i], w1);
+ L1Filter(T, sta[i].backInfo, v[i], dv2[i], w1);
b2 = b1;
z2 = z1;
@@ -345,7 +363,7 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
i = nHits - 1;
- FilterFirst(T, x_last, y_last, staLast);
+ FilterFirst(T, x_last, y_last, lstC00, lstC10, lstC11);
z1 = z[i];
sta[i].fieldSlice.GetFieldValue(T.x, T.y, b1);
@@ -375,8 +393,8 @@ void CbmL1PFFitter::Fit(vector<CbmStsTrack>& Tracks, vector<int>& pidHypo)
fit.L1AddMaterial(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, wIn);
fit.EnergyLossCorrection(T, CbmL1::Instance()->fpAlgo->GetParameters()->GetMaterialThickness(i, T.x, T.y), qp0, 1,
wIn);
- L1Filter(T, sta[i].frontInfo, u[i], w1);
- L1Filter(T, sta[i].backInfo, v[i], w1);
+ L1Filter(T, sta[i].frontInfo, u[i], du2[i], w1);
+ L1Filter(T, sta[i].backInfo, v[i], dv2[i], w1);
b2 = b1;
z2 = z1;