diff --git a/core/base/CbmTrackingDetectorInterfaceBase.h b/core/base/CbmTrackingDetectorInterfaceBase.h
index 48807029ec3622fd23141ec0c30c8ad5bb02891d..33cfb21e19676dbeec1777ea3b6600b832be07f5 100644
--- a/core/base/CbmTrackingDetectorInterfaceBase.h
+++ b/core/base/CbmTrackingDetectorInterfaceBase.h
@@ -12,6 +12,9 @@
#ifndef CbmTrackingDetectorInterfaceBase_h
#define CbmTrackingDetectorInterfaceBase_h 1
+#include "CbmHit.h"
+#include "CbmPixelHit.h"
+
#include <array>
#include <string>
#include <tuple>
@@ -20,8 +23,6 @@
#include <cmath>
-class CbmHit;
-
/// @class CbmTrackingDetectorInterfaceBase
/// @brief Abstract class, which should be inherited by every detecting subsystem tracking interface class
///
@@ -101,6 +102,15 @@ public:
/// @return Flag: true - station provides time measurements, false - station does not provide time measurements
virtual bool IsTimeInfoProvided(int stationId) const = 0;
+ /// @brief Gets x,y,t ranges of a CbmPixelHit
+ /// @param hit A hit
+ /// @return range X, Y, T
+ std::tuple<double, double, double> GetHitRanges(const CbmPixelHit& hit) const
+ {
+ // by default assume gaussian distributions of errors
+ return std::tuple(3.5 * hit.GetDx(), 3.5 * hit.GetDy(), 3.5 * hit.GetTimeError());
+ }
+
/// Technical flag: true - all downcasts are done with dynamic_cast followed by checks for nullptr (increases
/// computing time, better for debug), false - all downcasts are done with static_cast without sanity checks
/// (decreases computing time, but can cause bugs)
diff --git a/core/detectors/trd/CbmTrdTrackingInterface.cxx b/core/detectors/trd/CbmTrdTrackingInterface.cxx
index 1d2d4e27db8c491be426bb254d7df0805330a81c..1d5ff0d3820d85f4bc075edc32ae3c0879bc1b10 100644
--- a/core/detectors/trd/CbmTrdTrackingInterface.cxx
+++ b/core/detectors/trd/CbmTrdTrackingInterface.cxx
@@ -11,6 +11,8 @@
#include "CbmTrdTrackingInterface.h"
+#include "CbmTrdHit.h"
+
#include "FairDetector.h"
#include "FairRunAna.h"
#include <Logger.h>
@@ -71,6 +73,34 @@ std::tuple<double, double> CbmTrdTrackingInterface::GetStereoAnglesSensor(int ad
return std::tuple(0., TMath::Pi() / 2.);
}
+std::tuple<double, double, double> CbmTrdTrackingInterface::GetHitRanges(const CbmPixelHit& hit) const
+{
+ const CbmTrdHit* trdHit = dynamic_cast<const CbmTrdHit*>(&hit);
+ if (!trdHit) { LOG(fatal) << "CbmTrdTrackingInterface::GetHitRange: input hit is not of type CbmTrdHit"; }
+ const CbmTrdParModDigi* par = dynamic_cast<const CbmTrdParModDigi*>(fTrdDigiPar->GetModulePar(trdHit->GetAddress()));
+
+ auto [rangeX, rangeY, rangeT] = CbmTrackingDetectorInterfaceBase::GetHitRanges(hit);
+
+ if (trdHit->GetClassType() == 1) {
+ // TRD 2D hit
+ // TODO: replace with more sophisticated values once the errors are correct
+ rangeX = 1.7;
+ rangeY = 3.5;
+ rangeT = 200.;
+ }
+ else { // TRD 1D hit
+ if ((par->GetOrientation() == 1) || (par->GetOrientation() == 3)) {
+ // Y coordinate is uniformly distributed
+ rangeY = sqrt(3.) * trdHit->GetDy();
+ }
+ else { // X coordinate is uniformly distributed
+ rangeX = sqrt(3.) * trdHit->GetDx();
+ }
+ }
+
+ return std::tuple(rangeX, rangeY, rangeT);
+}
+
//-------------------------------------------------------------------------------------------------------------------------------------
//
InitStatus CbmTrdTrackingInterface::Init()
diff --git a/core/detectors/trd/CbmTrdTrackingInterface.h b/core/detectors/trd/CbmTrdTrackingInterface.h
index bd5828265c0212b2c49c1cbb1247621a170b3323..c1215c303f40ee5a3a5805b521731e0bc4c9a31f 100644
--- a/core/detectors/trd/CbmTrdTrackingInterface.h
+++ b/core/detectors/trd/CbmTrdTrackingInterface.h
@@ -113,6 +113,11 @@ public:
/// @return Flag: true - station provides time measurements, false - station does not provide time measurements
bool IsTimeInfoProvided(int /*stationId*/) const { return true; }
+ /// @brief Gets x,y,t ranges of a CbmTrdHit
+ /// @param hit A hit
+ /// @return range X, Y, T
+ std::tuple<double, double, double> GetHitRanges(const CbmPixelHit& hit) const;
+
/// @brief FairTask: Init method
InitStatus Init();
diff --git a/reco/L1/CbmCaTimeSliceReader.cxx b/reco/L1/CbmCaTimeSliceReader.cxx
index b8c99253113c01ea16e37a9e47c34bb74e56b6ff..b47900dccecb67d16432900453f7bfbe6ec49fbe 100644
--- a/reco/L1/CbmCaTimeSliceReader.cxx
+++ b/reco/L1/CbmCaTimeSliceReader.cxx
@@ -396,10 +396,15 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
aHit.y = hitRecord.fY;
aHit.z = hitRecord.fZ;
aHit.t = hitRecord.fT;
- aHit.dx = hitRecord.fDx;
+ aHit.dx2 = hitRecord.fDx2;
aHit.dxy = hitRecord.fDxy;
- aHit.dy = hitRecord.fDy;
- aHit.dt2 = hitRecord.fDt * hitRecord.fDt;
+ aHit.dy2 = hitRecord.fDy2;
+ aHit.dt2 = hitRecord.fDt2;
+
+ aHit.rangeX = hitRecord.fRangeX;
+ aHit.rangeY = hitRecord.fRangeY;
+ aHit.rangeT = hitRecord.fRangeT;
+
aHit.ID = static_cast<int>(fpIODataManager->GetNofHits());
aHit.iSt = hitRecord.fStaId;
fpIODataManager->PushBackHit(aHit, hitRecord.fDataStream);
@@ -417,8 +422,8 @@ void TimeSliceReader::StoreHitRecord(const HitRecord& hitRecord)
aHitQa.iStation = hitRecord.fStaId;
aHitQa.x = hitRecord.fX;
aHitQa.y = hitRecord.fY;
- aHitQa.dx = hitRecord.fDx;
- aHitQa.dy = hitRecord.fDy;
+ aHitQa.dx = sqrt(hitRecord.fDx2);
+ aHitQa.dy = sqrt(hitRecord.fDy2);
aHitQa.dxy = hitRecord.fDxy;
aHitQa.time = hitRecord.fT;
fpvQaHits->push_back(aHitQa);
diff --git a/reco/L1/CbmCaTimeSliceReader.h b/reco/L1/CbmCaTimeSliceReader.h
index 81eb53e4c29bcce0f42b8c400d253c79384329f1..966abdefefbb644850d529a47e13680662149f58 100644
--- a/reco/L1/CbmCaTimeSliceReader.h
+++ b/reco/L1/CbmCaTimeSliceReader.h
@@ -160,14 +160,14 @@ namespace cbm::ca
/// Stores recorded hit information into registered hit containers
void StoreHitRecord(const HitRecord& hitRecord);
- bool fbReadTracks = true; ///< flag to read reconstructed tracks from reco.root
+ bool fbReadTracks = true; ///< flag to read reconstructed tracks from reco.root
/// @brief Pointers to the tracking detector interfaces for each subsystem
DetIdArr_t<const CbmTrackingDetectorInterfaceBase*> fvpDetInterface = {{nullptr}};
// Input data branches
- CbmTimeSlice* fpBrTimeSlice = nullptr; ///< Pointer to the TS object
- DetIdArr_t<TClonesArray*> fvpBrHits = {{nullptr}}; ///< Input branch for hits
+ CbmTimeSlice* fpBrTimeSlice = nullptr; ///< Pointer to the TS object
+ DetIdArr_t<TClonesArray*> fvpBrHits = {{nullptr}}; ///< Input branch for hits
// Branches for reconstructed tracks. The input at the moment (as for 27.02.2023) depends on the selected
// tracking mode. For simulations in CBM, the CA tracking is used only in STS + MVD detectors. In this case
@@ -271,15 +271,18 @@ int cbm::ca::TimeSliceReader::ReadHitsForDetector()
if (iStActive == -1) { continue; } // Cut off inactive stations
// Fill out data common for all the detectors
- hitRecord.fStaId = iStActive;
- hitRecord.fX = pPixelHit->GetX();
- hitRecord.fY = pPixelHit->GetY();
- hitRecord.fZ = pPixelHit->GetZ();
- hitRecord.fDx = pPixelHit->GetDx();
- hitRecord.fDy = pPixelHit->GetDy();
- hitRecord.fDxy = pPixelHit->GetDxy();
- hitRecord.fT = pPixelHit->GetTime();
- hitRecord.fDt = pPixelHit->GetTimeError();
+ hitRecord.fStaId = iStActive;
+ hitRecord.fX = pPixelHit->GetX();
+ hitRecord.fY = pPixelHit->GetY();
+ hitRecord.fZ = pPixelHit->GetZ();
+ hitRecord.fDx2 = pPixelHit->GetDx() * pPixelHit->GetDx();
+ hitRecord.fDy2 = pPixelHit->GetDy() * pPixelHit->GetDy();
+ hitRecord.fDxy = pPixelHit->GetDxy();
+ hitRecord.fT = pPixelHit->GetTime();
+ hitRecord.fDt2 = pPixelHit->GetTimeError() * pPixelHit->GetTimeError();
+
+ std::tie(hitRecord.fRangeX, hitRecord.fRangeY, hitRecord.fRangeT) = pDetInterface->GetHitRanges(*pPixelHit);
+
hitRecord.fDet = static_cast<int>(DetID);
hitRecord.fDataStream = (static_cast<int64_t>(hitRecord.fDet) << 60) | pPixelHit->GetAddress();
hitRecord.fExtId = iHext;
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 46ade0005137adbf622c500e526d146894a1b184..d8248eb8440d86d2452c542c4246aef128ab056a 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -74,12 +74,16 @@ struct TmpHit {
double x; ///< position of hit along x axis [cm]
double y; ///< position of hit along y axis [cm]
double z; ///< position of hit along z axis [cm]
- double dx; ///< hit position uncertainty along x axis [cm]
- double dy; ///< hit position uncertainty along y axis [cm]
+ double dx2; ///< hit position uncertainty along x axis [cm]
+ double dy2; ///< hit position uncertainty along y axis [cm]
double dxy; ///< hit position covariance along x and y axes [cm2]
int iMC; ///< index of MCPoint in the fvMCPoints array
double time = 0.; ///< time of the hit [ns]
- double dt = 1.e4; ///< time error of the hit [ns]
+ double dt2 = 1.e8; ///< time error of the hit [ns]
+ double rangeX; ///< hit range in X [cm]
+ double rangeY; ///< hit range in Y [cm]
+ double rangeT; ///< hit range in T [ns]
+
int Det;
};
@@ -455,9 +459,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
h->Position(pos);
h->PositionError(err);
- th.dx = h->GetDx();
- th.dy = h->GetDy();
-
+ th.dx2 = h->GetDx() * h->GetDx();
+ th.dy2 = h->GetDy() * h->GetDy();
th.dxy = h->GetDxy();
th.x = pos.X();
@@ -465,8 +468,10 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.z = pos.Z();
// Get time
- th.time = h->GetTime(); // currently ignored by the tracking
- th.dt = h->GetTimeError(); // currently ignored by the tracking
+ th.time = h->GetTime(); // currently ignored by the tracking
+ th.dt2 = h->GetTimeError() * h->GetTimeError(); // currently ignored by the tracking
+
+ std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmMvdTrackingInterface::Instance()->GetHitRanges(*h);
}
th.Det = 0;
th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMvd>(hitIndex) : -1;
@@ -517,7 +522,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
//Get time
th.time = h->GetTime();
- th.dt = h->GetTimeError();
+ th.dt2 = h->GetTimeError() * h->GetTimeError();
TVector3 pos, err;
h->Position(pos);
@@ -527,10 +532,13 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.y = pos.Y();
th.z = pos.Z();
- th.dx = h->GetDx();
- th.dy = h->GetDy();
+ th.dx2 = h->GetDx() * h->GetDx();
+ th.dy2 = h->GetDy() * h->GetDy();
th.dxy = h->GetDxy();
+
+ std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmStsTrackingInterface::Instance()->GetHitRanges(*h);
}
+
th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kSts>(hitIndex) : -1;
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
tmpHits.push_back(th);
@@ -565,8 +573,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
//Get time
th.time = h->GetTime();
- th.dt = h->GetTimeError();
-
+ th.dt2 = h->GetTimeError() * h->GetTimeError();
// th.iSector = 0;
th.iStripF = firstDetStrip + hitIndex;
@@ -577,9 +584,11 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.y = h->GetY();
th.z = h->GetZ();
- th.dx = h->GetDx();
- th.dy = h->GetDy();
+ th.dx2 = h->GetDx() * h->GetDx();
+ th.dy2 = h->GetDy() * h->GetDy();
th.dxy = 0; /// FIXME: ???
+
+ std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmMuchTrackingInterface::Instance()->GetHitRanges(*h);
}
th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kMuch>(th.ExtIndex) : -1;
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
@@ -627,7 +636,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
// if (h->GetPlaneId()==3) continue;
th.time = h->GetTime();
- th.dt = h->GetTimeError();
+ th.dt2 = h->GetTimeError() * h->GetTimeError();
// th.iSector = 0;
th.iStripF = NStrips;
@@ -642,10 +651,12 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.y = pos.Y();
th.z = pos.Z();
- th.dx = fabs(h->GetDx());
- th.dy = fabs(h->GetDy());
+ th.dx2 = h->GetDx() * h->GetDx();
+ th.dy2 = h->GetDy() * h->GetDy();
th.dxy = 0;
+ std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmTrdTrackingInterface::Instance()->GetHitRanges(*h);
+
th.iMC = fPerformance ? MatchHitWithMc<L1DetectorID::kTrd>(hitIndex) : -1;
th.fDataStream = (static_cast<int64_t>(th.Det) << 60) | h->GetAddress();
@@ -685,10 +696,10 @@ void CbmL1::ReadEvent(CbmEvent* event)
if (th.iStation < 0) continue;
th.time = h->GetTime();
- th.dt = h->GetTimeError();
+ th.dt2 = h->GetTimeError() * h->GetTimeError();
- th.dx = h->GetDx();
- th.dy = h->GetDy();
+ th.dx2 = h->GetDx() * h->GetDx();
+ th.dy2 = h->GetDy() * h->GetDy();
th.dxy = h->GetDxy();
// th.iSector = 0;
@@ -704,6 +715,7 @@ void CbmL1::ReadEvent(CbmEvent* event)
th.y = pos.Y();
th.z = pos.Z();
+ std::tie(th.rangeX, th.rangeY, th.rangeT) = CbmTofTrackingInterface::Instance()->GetHitRanges(*h);
//TODO: is it still needed here? (S.Zharko)
if (L1Algo::TrackingMode::kMcbm == fTrackingMode && th.z > 400) continue;
@@ -751,8 +763,8 @@ void CbmL1::ReadEvent(CbmEvent* event)
s.iStation = th.iStation;
s.x = th.x;
s.y = th.y;
- s.dx = th.dx;
- s.dy = th.dy;
+ s.dx = sqrt(th.dx2);
+ s.dy = sqrt(th.dy2);
s.dxy = th.dxy;
s.time = th.time;
@@ -767,10 +779,15 @@ void CbmL1::ReadEvent(CbmEvent* event)
h.y = th.y;
h.z = th.z;
h.t = th.time;
- h.dx = th.dx;
- h.dy = th.dy;
+ h.dx2 = th.dx2;
+ h.dy2 = th.dy2;
h.dxy = th.dxy;
- h.dt2 = th.dt * th.dt;
+ h.dt2 = th.dt2;
+
+ h.rangeX = th.rangeX;
+ h.rangeY = th.rangeY;
+ h.rangeT = th.rangeT;
+
h.ID = iHit;
h.iSt = th.iStation;
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 875c9d617fd1d3254eea14e6e617f6c27adc065d..d54967897f05dc1f54f011d4063b32907a478c99 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -285,9 +285,9 @@ public:
L1Grid vGrid[L1Constants::size::kMaxNstations]; ///<
L1Grid vGridTime[L1Constants::size::kMaxNstations]; ///<
- fscal fMaxDx[L1Constants::size::kMaxNstations];
- fscal fMaxDy[L1Constants::size::kMaxNstations];
- fscal fMaxDt[L1Constants::size::kMaxNstations];
+ fscal fMaxRangeX[L1Constants::size::kMaxNstations];
+ fscal fMaxRangeY[L1Constants::size::kMaxNstations];
+ fscal fMaxRangeT[L1Constants::size::kMaxNstations];
double fCaRecoTime {0.}; // time of the track finder + fitter
diff --git a/reco/L1/L1Algo/L1BranchExtender.cxx b/reco/L1/L1Algo/L1BranchExtender.cxx
index d7d7e2fd893595b1a058814bbe265b51e1821ed8..d30eac508306e52ce54d990508c530851448d172 100644
--- a/reco/L1/L1Algo/L1BranchExtender.cxx
+++ b/reco/L1/L1Algo/L1BranchExtender.cxx
@@ -80,9 +80,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& Tout, const bool up
T.ResetErrors(1., 1., .1, .1, 1., (sta0.timeInfo ? hit0.dt2 : 1.e6), 1.e6);
T.NDF = fvec(2.);
- T.C00 = hit0.dx * hit0.dx;
+ T.C00 = hit0.dx2;
T.C10 = hit0.dxy;
- T.C11 = hit0.dy * hit0.dy;
+ T.C11 = hit0.dy2;
L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
L1FieldRegion fld _fvecalignment;
@@ -203,8 +203,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool upstream, co
const fscal iz = 1.f / (tr.z[0] - fParameters.GetTargetPositionZ()[0]);
L1HitAreaTime area(vGridTime[ista], tr.x[0] * iz, tr.y[0] * iz,
- (sqrt(fPickGather * tr.C00) + fMaxDx[ista] + fMaxDZ * abs(tr.tx))[0] * iz,
- (sqrt(fPickGather * tr.C11) + fMaxDy[ista] + fMaxDZ * abs(tr.ty))[0] * iz, tr.t[0],
+ (sqrt(fPickGather * tr.C00) + fMaxRangeX[ista] + fMaxDZ * abs(tr.tx))[0] * iz,
+ (sqrt(fPickGather * tr.C11) + fMaxRangeY[ista] + fMaxDZ * abs(tr.ty))[0] * iz, tr.t[0],
sqrt(tr.C55[0]));
for (L1HitIndex_t ih = -1; true;) { // loop over the hits in the area
@@ -226,7 +226,7 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool upstream, co
if (sta.timeInfo && tr.nTimeMeasurements[0] > 0) {
fscal dt = hit.t - tr.t[0];
- if (dt * dt > (tr.C55[0] + hit.dt2) * 25) continue;
+ if (fabs(dt) > sqrt(25. * tr.C55[0]) + hit.rangeT) continue;
}
//if (GetFUsed((*fStripFlag)[hit.f] | (*fStripFlag)[hit.b])) continue; // if used
@@ -250,8 +250,8 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& Tout, const bool upstream, co
fscal d2 = d_x * d_x + d_y * d_y;
if (d2 > r2_best) continue;
- fscal dxm_est2 = (pickGather2 * (abs(C00 + fMaxDx[ista] * fMaxDx[ista])))[0];
- if (d_x * d_x > dxm_est2) continue;
+ fscal dxm_est = sqrt(pickGather2 * C00)[0] + fMaxRangeX[ista];
+ if (fabs(d_x) > dxm_est) continue;
r2_best = d2;
iHit_best = globalInd;
diff --git a/reco/L1/L1Algo/L1CaTrackFinder.cxx b/reco/L1/L1Algo/L1CaTrackFinder.cxx
index 589b6b3f5efc789b8b9cbebbe32c113dbecb2290..4d6e49acfae802be47495ed9519886bfe8001674 100644
--- a/reco/L1/L1Algo/L1CaTrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinder.cxx
@@ -85,7 +85,7 @@ void L1Algo::CaTrackFinder()
fscal timeOfFlightMax =
1.5 * l * sqrt(1. + L1TrackPar::kProtonMass * L1TrackPar::kProtonMass / minProtonMomentum / minProtonMomentum)
* L1TrackPar::kClightNsInv;
- fscal dt = 3.5 * sqrt(h.dt2);
+ fscal dt = h.rangeT;
L1HitTimeInfo& info = fHitTimeInfo[caHitId];
info.fEventTimeMin = h.t - dt - timeOfFlightMax;
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index ddbb447b2910005edb9c4c145bdb365594b33cbd..3b5e1574f2a4906480103acbbc790e0971e93785 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -184,9 +184,9 @@ void L1Algo::CaTrackFinderSlice()
L1HitIndex_t nGridHitsFilled = 0;
for (int iS = 0; iS < fParameters.GetNstationsActive(); ++iS) {
- fMaxDx[iS] = 0.;
- fMaxDy[iS] = 0.;
- fMaxDt[iS] = 0.;
+ fMaxRangeX[iS] = 0.;
+ fMaxRangeY[iS] = 0.;
+ fMaxRangeT[iS] = 0.;
bool timeUninitialised = 1;
fscal lasttime = 0;
fscal starttime = 0;
@@ -200,11 +200,9 @@ void L1Algo::CaTrackFinderSlice()
for (L1HitIndex_t ih = 0; ih < fSliceHitIds[iS].size(); ++ih) {
const L1Hit& h = fInputData.GetHit(fSliceHitIds[iS][ih]);
- fscal dt = sqrt(h.dt2);
-
- if (fMaxDx[iS] < h.dx) { fMaxDx[iS] = h.dx; }
- if (fMaxDy[iS] < h.dy) { fMaxDy[iS] = h.dy; }
- if (fMaxDt[iS] < dt) { fMaxDt[iS] = dt; }
+ if (fMaxRangeX[iS] < h.rangeX) { fMaxRangeX[iS] = h.rangeX; }
+ if (fMaxRangeY[iS] < h.rangeY) { fMaxRangeY[iS] = h.rangeY; }
+ if (fMaxRangeT[iS] < h.rangeT) { fMaxRangeT[iS] = h.rangeT; }
auto [x, y] = GetHitCoorOnGrid(h);
if (x < gridMinX) { gridMinX = x; }
diff --git a/reco/L1/L1Algo/L1Fit.cxx b/reco/L1/L1Algo/L1Fit.cxx
index b8ef95582e130e78f823074c72634b6f28ace13b..8895cc4062c9fe3c3f6e6685789afa29b78d18ec 100644
--- a/reco/L1/L1Algo/L1Fit.cxx
+++ b/reco/L1/L1Algo/L1Fit.cxx
@@ -310,9 +310,9 @@ void L1Fit::FilterXY(const L1XYMeasurementInfo& info, cnst& x, cnst& y)
void L1Fit::FilterHit(const L1Station& sta, const L1Hit& hit)
{
L1XYMeasurementInfo info;
- info.C00 = hit.dx * hit.dx;
+ info.C00 = hit.dx2;
info.C10 = hit.dxy;
- info.C11 = hit.dy * hit.dy;
+ info.C11 = hit.dy2;
FilterXY(info, hit.x, hit.y);
FilterTime(hit.t, hit.dt2, sta.timeInfo);
}
diff --git a/reco/L1/L1Algo/L1Hit.h b/reco/L1/L1Algo/L1Hit.h
index b3c375a84ad3717f9aec1d9323dd07d6e7483400..88999fa034b80bf58472a61986a8a0d5c71f1060 100644
--- a/reco/L1/L1Algo/L1Hit.h
+++ b/reco/L1/L1Algo/L1Hit.h
@@ -34,16 +34,20 @@ public:
/// tracking detectors (MVD, MuCh, TRD, TOF) f == b and corresponds to the index of the hit. Indexes f and b
/// do not intersect between different detector stations.
- float x {0.f}; ///< measured X coordinate [cm]
- float y {0.f}; ///< measured Y coordinate [cm]
- float z = 0.f; ///< fixed Z coordinate [cm]
- float t = 0.f; ///< measured time [ns]
- float dx {0.f}; ///< uncertainty of X coordinate [cm]
- float dy {0.f}; ///< uncertainty of Y coordinate [cm]
- float dxy {0.f}; ///< X/Y covariance [cm2]
- float dt2 = 0.f; ///< measured uncertainty of time [ns]
- int ID = 0; ///< index of hit before hit sorting
- int iSt = 0; ///< index of station in the active stations array
+ float x {0.f}; ///< measured X coordinate [cm]
+ float y {0.f}; ///< measured Y coordinate [cm]
+ float z = 0.f; ///< fixed Z coordinate [cm]
+ float t = 0.f; ///< measured time [ns]
+ float dx2 {0.f}; ///< rms^2 of uncertainty of X coordinate [cm]
+ float dy2 {0.f}; ///< rms^2 of uncertainty of Y coordinate [cm]
+ float dxy {0.f}; ///< X/Y covariance [cm2]
+ float dt2 {0.f}; ///< measured uncertainty of time [ns]
+ float rangeX {0.f}; ///< +-range of uncertainty of X coordinate [cm]
+ float rangeY {0.f}; ///< +-range of uncertainty of Y coordinate [cm]
+ float rangeT {0.f}; ///< +-range of uncertainty of time [ns]
+
+ int ID = 0; ///< index of hit before hit sorting
+ int iSt = 0; ///< index of station in the active stations array
private:
/// Serialization method, used to save L1Hit objects into binary or text file in a defined order
@@ -56,10 +60,13 @@ private:
ar& y;
ar& z;
ar& t;
- ar& dx;
+ ar& dx2;
ar& dxy;
- ar& dy;
+ ar& dy2;
ar& dt2;
+ ar& rangeX;
+ ar& rangeY;
+ ar& rangeT;
ar& ID;
ar& iSt;
}
diff --git a/reco/L1/L1Algo/L1HitPoint.h b/reco/L1/L1Algo/L1HitPoint.h
index 5ded6ff21997b52ed6c606b8f505f50ebbdb3de7..8cf9d874ced37eca2041386625f1dfba179fe7ba 100644
--- a/reco/L1/L1Algo/L1HitPoint.h
+++ b/reco/L1/L1Algo/L1HitPoint.h
@@ -17,14 +17,17 @@ struct L1HitPoint {
void Set(const L1Hit& hit)
{
- x = hit.x;
- y = hit.y;
- z = hit.z;
- t = hit.t;
- dx = hit.dx;
- dxy = hit.dxy;
- dy = hit.dy;
- dt2 = hit.dt2;
+ x = hit.x;
+ y = hit.y;
+ z = hit.z;
+ t = hit.t;
+ dx2 = hit.dx2;
+ dxy = hit.dxy;
+ dy2 = hit.dy2;
+ dt2 = hit.dt2;
+ rangeX = hit.rangeX;
+ rangeY = hit.rangeY;
+ rangeT = hit.rangeT;
fIsSuppressed = 0;
}
@@ -33,19 +36,20 @@ struct L1HitPoint {
fscal T() const { return t; }
fscal X() const { return x; }
fscal Y() const { return y; }
- fscal dX() const { return dx; }
+ fscal dX2() const { return dx2; }
fscal dXY() const { return dxy; }
- fscal dY() const { return dy; }
- fscal dX2() const { return dx * dx; }
- fscal dY2() const { return dy * dy; }
+ fscal dY2() const { return dy2; }
fscal dT2() const { return dt2; }
+ fscal RangeX() const { return rangeX; }
+ fscal RangeY() const { return rangeY; }
+ fscal RangeT() const { return rangeT; }
bool IsSuppressed() const { return fIsSuppressed; }
void SetIsSuppresed(bool val) { fIsSuppressed = val; }
private:
- fscal x {0.}, y {0.}, z {0.}, t {0.}, dx {0.}, dxy {0.}, dy {0.}, dt2 {0.};
+ fscal x {0.}, y {0.}, z {0.}, t {0.}, dx2 {0.}, dxy {0.}, dy2 {0.}, dt2 {0.}, rangeX {0.}, rangeY {0.}, rangeT {0.};
// fIsSuppressed flag is used to suppress duplicated hits at the module overlaps
// TODO: collect those hits on the track instead of suppressing them
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index d5058c45a8fb61a1313ec1daff5649be302f0ae7..c3c0a0b730139aad9da6b5d00c516c0baf3aca3d 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -145,9 +145,9 @@ void L1Algo::L1KFTrackFitter()
if (!sta[ista].timeInfo) { dt2[ista][iVec] = 1.e4; }
z[ista][iVec] = hit.z;
sta[ista].fieldSlice.GetFieldValue(x[ista], y[ista], fB_temp);
- cov_xy[ista].C00[iVec] = hit.dx * hit.dx;
+ cov_xy[ista].C00[iVec] = hit.dx2;
cov_xy[ista].C10[iVec] = hit.dxy;
- cov_xy[ista].C11[iVec] = hit.dy * hit.dy;
+ cov_xy[ista].C11[iVec] = hit.dy2;
fB[ista].x[iVec] = fB_temp.x[iVec];
fB[ista].y[iVec] = fB_temp.y[iVec];
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx
index 29e6bff76a0a868e7b12885239f409a1f6e4ad97..ad90a88bf57f33aa7b5c476acf3a74f10d30d199 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.cxx
+++ b/reco/L1/L1Algo/L1TripletConstructor.cxx
@@ -241,7 +241,9 @@ void L1TripletConstructor::FitDoublets()
fFit.SetQp0(fFit.Tr().qp);
// check if there are other hits close to the doublet on the same station
- {
+ if (L1Algo::kMcbm != fAlgo->fTrackingMode) {
+ // TODO: SG: adjust cuts, put them to parameters
+
fscal tx = T2.tx[0];
fscal ty = T2.ty[0];
fscal tt = T2.vi[0] * sqrt(1. + tx * tx + ty * ty); // dt/dl * dl/dz
@@ -254,14 +256,14 @@ void L1TripletConstructor::FitDoublets()
if ((staM().timeInfo) && (T2.nTimeMeasurements[0] > 0)) {
fscal dt = T2.t[0] + tt * dz - hitm.T();
- if (dt * dt > 30. * (T2.C55[0] + hitm1.dT2())) { continue; }
+ if (fabs(dt) > sqrt(30. * T2.C55[0]) + hitm1.RangeT()) { continue; }
}
fscal dx = T2.x[0] + tx * dz - hitm1.X();
- if (dx * dx > 20. * (T2.C00[0] + hitm1.dX2())) { continue; }
+ if (fabs(dx) > 1.3 * (3.5 * sqrt(T2.C00[0]) + hitm1.RangeX())) { continue; }
fscal dy = T2.y[0] + ty * dz - hitm1.Y();
- if (dy * dy > 30. * (T2.C11[0] + hitm1.dY2())) { continue; }
+ if (fabs(dy) > 1.6 * (3.5 * sqrt(T2.C11[0]) + hitm1.RangeY())) { continue; }
if (fAlgo->fParameters.DevIsSuppressOverlapHitsViaMc()) {
int indL = fAlgo->fGridHitStartIndex[fIstaL] + fIhitL;
@@ -422,9 +424,9 @@ void L1TripletConstructor::FitTriplets()
y[ih] = hit.y;
z[ih] = hit.z;
t[ih] = hit.t;
- cov[ih].C00 = hit.dx * hit.dx;
+ cov[ih].C00 = hit.dx2;
cov[ih].C10 = hit.dxy;
- cov[ih].C11 = hit.dy * hit.dy;
+ cov[ih].C11 = hit.dy2;
dt2[ih] = hit.dt2;
};
@@ -612,7 +614,8 @@ void L1TripletConstructor::CollectHits(const L1TrackPar& Tr, const int iSta, con
fFit.SetTrack(Tr);
L1TrackPar& T = fFit.Tr();
-
+ //std::cout << T.chi2[0] << std::endl;
+ T.chi2 = 0.;
// if make it bigger the found hits will be rejected later because of the chi2 cut.
const fvec Pick_m22 = (fvec(chi2Cut) - T.chi2);
@@ -621,9 +624,9 @@ void L1TripletConstructor::CollectHits(const L1TrackPar& Tr, const int iSta, con
const fscal time = T.t[0];
L1HitAreaTime areaTime(fAlgo->vGridTime[iSta], T.x[0] * iz, T.y[0] * iz,
- (sqrt(Pick_m22 * T.C00) + fAlgo->fMaxDx[iSta] + fAlgo->fMaxDZ * abs(T.tx))[0] * iz,
- (sqrt(Pick_m22 * T.C11) + fAlgo->fMaxDy[iSta] + fAlgo->fMaxDZ * abs(T.ty))[0] * iz, time,
- sqrt(timeError2) * 5);
+ (sqrt(Pick_m22 * T.C00) + fAlgo->fMaxRangeX[iSta] + fAlgo->fMaxDZ * abs(T.tx))[0] * iz,
+ (sqrt(Pick_m22 * T.C11) + fAlgo->fMaxRangeY[iSta] + fAlgo->fMaxDZ * abs(T.ty))[0] * iz, time,
+ sqrt(timeError2) * 5 + fAlgo->fMaxRangeT[iSta]);
for (L1HitIndex_t ih = 0; (ih < nHits) && ((int) collectedHits.size() < maxNhits);
ih++) { // loop over all station hits
@@ -642,7 +645,7 @@ void L1TripletConstructor::CollectHits(const L1TrackPar& Tr, const int iSta, con
// check y-boundaries
//TODO: move hardcoded cuts to parameters
if ((sta.timeInfo) && (T.nTimeMeasurements[0] > 0)) {
- if (fabs(time - hit.T()) > sqrt(timeError2 + hit.dT2()) * 5) continue;
+ if (fabs(time - hit.T()) > 1.4 * (3.5 * sqrt(timeError2) + hit.RangeT())) continue;
if (fabs(time - hit.T()) > 30) continue;
}
@@ -652,24 +655,24 @@ void L1TripletConstructor::CollectHits(const L1TrackPar& Tr, const int iSta, con
fvec y, C11;
fFit.ExtrapolateYC11Line(z, y, C11);
- fscal dy_est2 = Pick_m22[0] * fabs(C11[0] + hit.dY2());
+ fscal dy_est = sqrt(Pick_m22[0] * C11[0]) + hit.RangeY();
fscal xm = hit.X();
fscal ym = hit.Y();
const fscal dY = ym - y[0];
- if (dY * dY > dy_est2) continue;
+ if (fabs(dY) > dy_est) continue;
// check x-boundaries
fvec x, C00;
fFit.ExtrapolateXC00Line(z, x, C00);
- fscal dx_est2 = Pick_m22[0] * fabs(C00[0] + hit.dX2());
+ fscal dx_est = sqrt(Pick_m22[0] * C00[0]) + hit.RangeX();
const fscal dX = xm - x[0];
- if (dX * dX > dx_est2) continue;
+ if (fabs(dX) > dx_est) continue;
// check chi2
diff --git a/reco/L1/catools/CaToolsHitRecord.h b/reco/L1/catools/CaToolsHitRecord.h
index e7cba94cb9720b18ec04e134123e01d9780c7a94..c5d91fd958deb25cc7e7aa5f2a8444e067054c97 100644
--- a/reco/L1/catools/CaToolsHitRecord.h
+++ b/reco/L1/catools/CaToolsHitRecord.h
@@ -24,10 +24,13 @@ namespace ca::tools
double fY = -1.; ///< y component of hit position [cm]
double fZ = -1.; ///< z component of hit position [cm]
double fT = -1.; ///< time of hit [ns]
- double fDx = -1.; ///< error of x component of hit position [cm]
- double fDy = -1.; ///< error of y component of hit position [cm]
+ double fDx2 = -1.; ///< error of x component of hit position [cm]
+ double fDy2 = -1.; ///< error of y component of hit position [cm]
double fDxy = -1.; ///< correlation between x and y components [cm]
- double fDt = -1.; ///< time error of hit [ns]
+ double fDt2 = -1.; ///< time error of hit [ns]
+ double fRangeX = -1.; ///< range of x [cm]
+ double fRangeY = -1.; ///< range of y [cm]
+ double fRangeT = -1.; ///< range of t [ns]
int64_t fDataStream = -1; ///< Global index of detector module
int fPointId = -1; ///< index of MC point
int fExtId = -2; ///< external index of hit
diff --git a/reco/L1/qa/CbmCaInputQaTrd.cxx b/reco/L1/qa/CbmCaInputQaTrd.cxx
index bf33aedf723cf5c09906dd20beafa4cd12fe80b5..f3a8e536fd40f6dc0bb3131db7e277afe38a873b 100644
--- a/reco/L1/qa/CbmCaInputQaTrd.cxx
+++ b/reco/L1/qa/CbmCaInputQaTrd.cxx
@@ -46,13 +46,13 @@ void CbmCaInputQaTrd::DefineParameters()
SetRange(fRHitDy, 0.0000, 5.00); // [cm]
SetRange(fRHitDu, 0.0000, 5.00); // [cm]
SetRange(fRHitDv, 0.0000, 5.00); // [cm]
- SetRange(fRHitDt, 0.0000, 10.00); // [ns]
+ SetRange(fRHitDt, 0.0000, 100.00); // [ns]
// Residuals
SetRange(fRResX, -10.00, 10.00);
SetRange(fRResY, -10.00, 10.00);
SetRange(fRResU, -2.00, 2.00);
SetRange(fRResV, -10.00, 10.00);
- SetRange(fRResT, -0.50, 0.50);
+ SetRange(fRResT, -200.0, 200.0);
// QA result selection criteria
SetRange(fEffRange, 10.0, 30.0); ///< Range for hit efficiency approximation
fResMeanThrsh = 0.50; ///< Maximum allowed deviation of residual mean from zero [sigma]