diff --git a/macro/L1/run_reco_L1global.C b/macro/L1/run_reco_L1global.C
index 9162f7cbe34e858338ad9026c0e178e124b5b108..05cfa152d7f99b2800a957f98320174423e3cb91 100644
--- a/macro/L1/run_reco_L1global.C
+++ b/macro/L1/run_reco_L1global.C
@@ -360,7 +360,7 @@ void run_reco_L1global(TString input = "", Int_t nTimeSlices = -1, Int_t firstTi
}
l1->SetGlobalMode();
l1->SetUseHitErrors(true);
- l1->SetUseMcHit(0, 0, 1, 0);
+ //l1->SetUseMcHit(0, 0, 1, 0);
// --- Material budget file names
TString mvdGeoTag;
@@ -397,60 +397,6 @@ void run_reco_L1global(TString input = "", Int_t nTimeSlices = -1, Int_t firstTi
// ------------------------------------------------------------------------
- // ==== From here on, the time-based and the event-based reconstruction
- // ==== chains differ, since time-based version of primary vertex finding
- // ==== and global tracking are not yet available. For time-based
- // ==== reconstruction, a track-based event finder is used; no global
- // ==== tracks are produced.
-
- if (eventBased) {
-
- // ----- Primary vertex finding -------------------------------------
- CbmPrimaryVertexFinder* pvFinder = new CbmPVFinderKF();
- CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
- run->AddTask(findVertex);
- std::cout << "-I- " << myName << ": Added task " << findVertex->GetName() << std::endl;
- // ----------------------------------------------------------------------
-
-
- // --- Global track finding -----------------------------------------
- CbmLitFindGlobalTracks* finder = new CbmLitFindGlobalTracks();
- finder->SetTrackingType("branch");
- finder->SetMergerType("nearest_hit");
- run->AddTask(finder);
- std::cout << "-I- : Added task " << finder->GetName() << std::endl;
- // ----------------------------------------------------------------------
-
- // --- Particle Id in TRD -----------------------------------------
- if (useTrd) {
- CbmTrdSetTracksPidLike* trdLI = new CbmTrdSetTracksPidLike("TRDLikelihood", "TRDLikelihood");
- trdLI->SetUseMCInfo(kTRUE);
- trdLI->SetUseMomDependence(kTRUE);
- run->AddTask(trdLI);
- std::cout << "-I- : Added task " << trdLI->GetName() << std::endl;
- }
- // ------------------------------------------------------------------------
-
-
- // ----- RICH reconstruction ----------------------------------------
- if (useRich) {
- CbmRichReconstruction* richReco = new CbmRichReconstruction();
- run->AddTask(richReco);
- std::cout << "-I- : Added task " << richReco->GetName() << std::endl;
- }
- // ----------------------------------------------------------------------
-
- } //? event-based reco
-
- else {
-
- // -----Â Event building from STS tracks -----------------------------
- run->AddTask(new CbmBuildEventsFromTracksReal());
- // ----------------------------------------------------------------------
-
- } //? time-based reco
-
-
// ----- Parameter database --------------------------------------------
std::cout << std::endl << std::endl;
std::cout << "-I- " << myName << ": Set runtime DB" << std::endl;
diff --git a/reco/L1/CbmL1MCTrack.h b/reco/L1/CbmL1MCTrack.h
index 3e4d4ec2a860672efb9cf79afa13b7a777807bb5..2bb0d1a13881a71b4d31ee3d805bc6321b945e60 100644
--- a/reco/L1/CbmL1MCTrack.h
+++ b/reco/L1/CbmL1MCTrack.h
@@ -61,6 +61,8 @@ public:
friend class CbmL1;
+ double pt() { return sqrt(px * px + py * py); }
+
private:
void CalculateMCCont();
void CalculateHitCont();
@@ -85,8 +87,8 @@ public:
int mother_ID = -1;
int pdg = -1;
bool isSignal {0};
- L1Vector<int> Points {"CbmL1MCTrack::Points"}; // indices of pints in L1::vMCPoints
- L1Vector<int> Hits {"CbmL1MCTrack::Hits"}; // indices of hits in algo->vHits or L1::vHits
+ L1Vector<int> Points {"CbmL1MCTrack::Points"}; // indices of pints in L1::vMCPoints
+ L1Vector<int> Hits {"CbmL1MCTrack::Hits"}; // indices of hits in algo->vHits or L1::vHits
private:
int nMCContStations = 0; // number of consecutive stations with mcPoints
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index e357801033ff1c79bf9eea14d7f144572375a340..3468efb009ddf0c772930793212699b4d6b04486 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -305,19 +305,19 @@ void CbmL1::EfficienciesPerformance()
static double L1_CATIME = 0.0;
- TL1PerfEfficiencies ntra; // efficiencies for current event
+ TL1PerfEfficiencies ntra; // efficiencies for current event
ntra.fOutDir = fHistoDir; // Setup a pointer for output directory
L1_NTRA.fOutDir = fHistoDir; // Save average efficiencies
for (vector<CbmL1Track>::iterator rtraIt = vRTracks.begin(); rtraIt != vRTracks.end(); ++rtraIt) {
ntra.ghosts += rtraIt->IsGhost();
- // if(rtraIt->IsGhost()){ // Debug.
- // cout << " " << rtraIt->GetNOfHits() << " " << 1./rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " | ";
+ //if (rtraIt->IsGhost()) { // Debug.
+ //cout << " " << rtraIt->GetNOfHits() << " " << 1. / rtraIt->T[5] << " " << rtraIt->GetMaxPurity() << " | ";
// for( map<int, int>::iterator posIt = rtraIt->hitMap.begin(); posIt != rtraIt->hitMap.end(); posIt++ ){
// cout << " (" << posIt->second << " " << posIt->first << ") ";
// }
- // cout << endl;
- // }
+ //cout << endl;
+ //}
}
int sta_nhits[algo->GetNstations()], sta_nfakes[algo->GetNstations()];
@@ -457,7 +457,7 @@ void CbmL1::EfficienciesPerformance()
// cout.precision(3);
if (fVerbose) {
if (fVerbose > 1) {
- ntra.PrintEff();
+ ntra.PrintEff(true);
cout << "Number of true and fake hits in stations: " << endl;
for (int i = 0; i < algo->GetNstations(); i++) {
cout << sta_nhits[i] - sta_nfakes[i] << "+" << sta_nfakes[i] << " ";
@@ -1191,10 +1191,10 @@ void CbmL1::TrackFitPerformance()
float z[3] = {0.f, 0.f, 0.f};
int ih = 0;
for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
- const int iMCP = mc.Points[iMCPoint];
- CbmL1MCPoint& mcP = vMCPoints[iMCP];
+ const int iMCP = mc.Points[iMCPoint];
+ CbmL1MCPoint& mcP = vMCPoints[iMCP];
const L1Station& st = algo->GetParameters()->GetStation(mcP.iStation);
- z[ih] = st.z[0];
+ z[ih] = st.z[0];
if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
ih++;
@@ -1286,8 +1286,8 @@ void CbmL1::TrackFitPerformance()
}
- { // last hit
- int iMC = vHitMCRef[it->Hits.back()]; // TODO2: adapt to linking
+ { // last hit
+ int iMC = vHitMCRef[it->Hits.back()]; // TODO2: adapt to linking
if (iMC < 0) continue;
#define L1FSTPARAMEXTRAPOLATE
@@ -1302,10 +1302,10 @@ void CbmL1::TrackFitPerformance()
float z[3] = {0.f, 0.f, 0.f};
int ih = 0;
for (unsigned int iMCPoint = 0; iMCPoint < mc.Points.size(); iMCPoint++) {
- const int iMCP = mc.Points[iMCPoint];
- CbmL1MCPoint& mcP = vMCPoints[iMCP];
+ const int iMCP = mc.Points[iMCPoint];
+ CbmL1MCPoint& mcP = vMCPoints[iMCP];
const L1Station& st = algo->GetParameters()->GetStation(mcP.iStation);
- z[ih] = st.z[0];
+ z[ih] = st.z[0];
if (ih > 0 && (z[ih] - z[ih - 1]) < 0.1) continue;
st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
ih++;
@@ -1373,10 +1373,10 @@ void CbmL1::TrackFitPerformance()
float z[3];
for (unsigned int ih = 0; ih < 3; ih++) {
if (ih >= mc.Points.size()) continue; //If nofMCPoints in track < 3
- const int iMCP = mc.Points[ih];
- CbmL1MCPoint& mcP = vMCPoints[iMCP];
+ const int iMCP = mc.Points[ih];
+ CbmL1MCPoint& mcP = vMCPoints[iMCP];
const L1Station& st = algo->GetParameters()->GetStation(mcP.iStation);
- z[ih] = st.z[0];
+ z[ih] = st.z[0];
st.fieldSlice.GetFieldValue(mcP.x, mcP.y, B[ih]);
};
fld.Set(B[0], z[0], B[1], z[1], B[2], z[2]);
@@ -1441,7 +1441,7 @@ void CbmL1::TrackFitPerformance()
for (unsigned int iHit = 0; iHit < it->Hits.size(); iHit++) {
const int iStation = vHitStore[it->Hits[iHit]].iStation;
const L1Station& st = algo->GetParameters()->GetStation(iStation);
- z[ih] = st.z[0];
+ z[ih] = st.z[0];
st.fieldSlice.GetFieldValue(vHitStore[it->Hits[iHit]].x, vHitStore[it->Hits[iHit]].y, B[ih]);
ih++;
if (ih == 3) break;
@@ -1618,8 +1618,8 @@ void CbmL1::FieldApproxCheck()
L1FieldValue B_L1;
Double_t bbb, bbb_L1;
- const int M = 5; // polinom order
- const int N = (M + 1) * (M + 2) / 2;
+ const int M = 5; // polinom order
+ const int N = (M + 1) * (M + 2) / 2;
const L1Station& st = algo->GetParameters()->GetStation(ist);
for (int i = 0; i < N; i++) {
FSl.cx[i] = st.fieldSlice.cx[i][0];
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 2d85d7a05112d6bb8143db4dfd6e661610857ab5..6c92b96a671541fe63db671af3ff88edff410283 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -151,14 +151,19 @@ struct TmpHit {
/// \param point constant reference to the input MC-point
/// \param st reference to the station info object
// TODO: Probably, L1Station& st parameter should be constant. Do we really want to modify a station here? (S.Zharko)
- void SetHitFromPoint(const CbmL1MCPoint& point, const L1Station& st)
+ void SetHitFromPoint(const CbmL1MCPoint& point, const L1Station& st, bool doSmear = 1)
{
- x = 0.5 * (point.xIn + point.xOut) + gRandom->Gaus(0, dx);
- y = 0.5 * (point.yIn + point.yOut) + gRandom->Gaus(0, dy);
+ x = 0.5 * (point.xIn + point.xOut);
+ y = 0.5 * (point.yIn + point.yOut);
z = 0.5 * (point.zIn + point.zOut);
- time = point.time + gRandom->Gaus(0, dt);
+ time = point.time;
u_front = x * st.frontInfo.cos_phi[0] + y * st.frontInfo.sin_phi[0];
u_back = x * st.backInfo.cos_phi[0] + y * st.backInfo.sin_phi[0];
+ if (doSmear) {
+ x += gRandom->Gaus(0, dx);
+ y += gRandom->Gaus(0, dy);
+ time += gRandom->Gaus(0, dt);
+ }
}
};
@@ -938,26 +943,14 @@ void CbmL1::ReadEvent(L1AlgoInputData* fData_, float& TsStart, float& TsLength,
if (th.iMC < 0) continue; // skip noise hits
if (mcUsed[iMcTrd]) continue; // one hit per MC point
-
- if (L1Algo::TrackingMode::kGlobal == fTrackingMode) {
- // SGtrd2d!! don't smear
- TmpHit thsave = th;
- th.dx = 0.;
- th.dy = 0.;
- th.dt = 0.;
- th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetParameters()->GetStation(th.iStation));
- th.dx = thsave.dx;
- th.dy = thsave.dy;
- th.dt = thsave.dt;
- // SGtrd2d!! artificial errors
- th.dx = 1.1;
- th.du = 1.1;
- th.dy = 1.1;
- th.dv = 1.1;
- }
- else {
- th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetParameters()->GetStation(th.iStation));
+ bool doSmear = true;
+ if (0) { // SGtrd2d!! debug: artificial errors
+ th.dx = .1;
+ th.du = .1;
+ th.dy = .1;
+ th.dv = .1;
}
+ th.SetHitFromPoint(vMCPoints[th.iMC], algo->GetParameters()->GetStation(th.iStation), doSmear);
mcUsed[iMcTrd] = 1;
} // replace hit by MC points
@@ -1271,7 +1264,6 @@ void CbmL1::Fill_vMCTracks()
TLorentzVector vp;
MCTrack->GetStartVertex(vr);
MCTrack->Get4Momentum(vp);
-
Int_t pdg = MCTrack->GetPdgCode();
Double_t q = 0, mass = 0.;
if (pdg < 9999999 && ((TParticlePDG*) TDatabasePDG::Instance()->GetParticle(pdg))) {
@@ -1279,6 +1271,7 @@ void CbmL1::Fill_vMCTracks()
mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
}
+ //cout << "mc track " << iMCTrack << " pdg " << pdg << " z " << vr.Z() << endl;
Int_t iTrack = vMCTracks.size(); //or iMCTrack?
CbmL1MCTrack T(mass, q, vr, vp, iTrack, mother_ID, pdg);
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 44667d32ec2110a4b15407703036901b1b26d0dc..067ab935ad0823f56d81fa7157fbb6c9e3b280a1 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -90,7 +90,7 @@ void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool M
fInitManager.TransferParametersContainer(fParameters);
- LOG(info) << fParameters.ToString();
+ LOG(info) << fParameters.ToString(3);
// Get number of station
fNstations = fInitManager.GetNstationsActive();
@@ -118,33 +118,31 @@ void L1Algo::SetData(L1Vector<L1Hit>& Hits_, int nStrips_, L1Vector<unsigned cha
NHitsIsecAll = nHits;
- vDontUsedHits_A.reset(nHits);
- vDontUsedHits_B.reset(nHits);
- vDontUsedHits_Buf.reset(nHits);
- vDontUsedHitsxy_A.reset(nHits);
- vDontUsedHitsxy_buf.reset(nHits);
- vDontUsedHitsxy_B.reset(nHits);
+ vNotUsedHits_A.reset(nHits);
+ vNotUsedHits_B.reset(nHits);
+ vNotUsedHits_Buf.reset(nHits);
+ vNotUsedHitsxy_A.reset(nHits);
+ vNotUsedHitsxy_buf.reset(nHits);
+ vNotUsedHitsxy_B.reset(nHits);
RealIHit_v.reset(nHits);
RealIHit_v_buf.reset(nHits);
RealIHit_v_buf2.reset(nHits);
#ifdef _OPENMP
- fHitToBestTrackF.reset(fNstrips);
- fHitToBestTrackB.reset(fNstrips);
- for (unsigned int j = 0; j < fHitToBestTrackB.size(); j++) {
- omp_init_lock(&fHitToBestTrackB[j]);
- omp_init_lock(&fHitToBestTrackF[j]);
+ fStripToTrackLock.reset(fNstrips);
+ for (unsigned int j = 0; j < fStripToTrackLock.size(); j++) {
+ omp_init_lock(&fStripToTrackLock[j]);
}
#endif
- fStripToTrack.clear();
- fStripToTrack.reserve(fNstrips);
+ fStripToTrackF.clear();
+ fStripToTrackF.reserve(fNstrips);
fStripToTrackB.clear();
fStripToTrackB.reserve(fNstrips);
- TripForHit[0].reset(nHits);
- TripForHit[1].reset(nHits);
+ fHitFirstTriplet.reset(nHits);
+ fHitNtriplets.reset(nHits);
fDupletPortionSize.clear();
fDupletPortionSize.reserve(2 * nHits);
@@ -244,6 +242,14 @@ void L1Algo::CreateHitPoint(const L1Hit& hit, L1HitPoint& point)
}
int L1Algo::GetMcTrackIdForHit(int iHit)
+{
+ int hitId = iHit;
+ int iMcPoint = CbmL1::Instance()->GetHitMCRefs()[hitId];
+ if (iMcPoint < 0) return -1;
+ return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
+}
+
+int L1Algo::GetMcTrackIdForUnusedHit(int iHit)
{
int hitId = (*RealIHitP)[iHit];
int iMcPoint = CbmL1::Instance()->GetHitMCRefs()[hitId];
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 8206baf28431fa8183e2cda828df341e85a2543f..433458ff01b5df3d61ad41b5a792c94b4f493924 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -65,6 +65,8 @@ class L1AlgoDraw;
using std::map;
+typedef int Tindex;
+
#ifdef PULLS
#define TRIP_PERFORMANCE
class L1AlgoPulls;
@@ -77,7 +79,6 @@ class L1AlgoEfficiencyPerformance;
template<Tindex NHits>
class L1AlgoEfficiencyPerformance;
#endif
-typedef int Tindex;
using L1StationsArray_t = std::array<L1Station, L1Constants::size::kMaxNstations>;
using L1MaterialArray_t = std::array<L1Material, L1Constants::size::kMaxNstations>;
@@ -109,11 +110,11 @@ public:
/// pack station, thread and triplet indices to an unique triplet ID
static unsigned int PackTripletId(unsigned int Station, unsigned int Thread, unsigned int Triplet)
{
-#ifndef FAST_CODE
+ //SG!!#ifndef FAST_CODE
assert(Station < L1Constants::size::kMaxNstations);
assert(Thread < L1Constants::size::kMaxNthreads);
assert(Triplet < L1Constants::size::kMaxNtriplets);
-#endif
+ //#endif
constexpr unsigned int kMoveThread = L1Constants::size::kTripletBits;
constexpr unsigned int kMoveStation = L1Constants::size::kTripletBits + L1Constants::size::kThreadBits;
return (Station << kMoveStation) + (Thread << kMoveThread) + Triplet;
@@ -142,6 +143,19 @@ public:
}
+ /// pack thread and track indices to an unique track ID
+ static int PackTrackId(int Thread, int Track)
+ {
+ return PackTripletId(0, (unsigned int) Thread, (unsigned int) Track);
+ }
+
+ /// unpack the track ID to its thread index
+ static int TrackId2Thread(int ID) { return TripletId2Thread((unsigned int) ID); }
+
+ /// unpack the track ID to its track index
+ static int TrackId2Track(int ID) { return TripletId2Triplet((unsigned int) ID); }
+
+
L1Vector<L1Triplet> fTriplets[L1Constants::size::kMaxNstations][L1Constants::size::kMaxNthreads] {
{"L1Algo::fTriplets"}}; // created triplets at station + thread
@@ -241,6 +255,9 @@ public:
/// Get mc track ID for a hit (debug tool)
int GetMcTrackIdForHit(int iHit);
+ /// Get mc track ID for a hit (debug tool)
+ int GetMcTrackIdForUnusedHit(int iHit);
+
public:
int fNstrips {0}; ///> number of strips
L1Vector<L1Hit>* vHits {nullptr}; ///> hits as a combination of front-, backstrips and z-position
@@ -261,12 +278,12 @@ public:
// L1Branch* pointer;
unsigned int NHitsIsecAll {0};
- L1Vector<L1Hit> vDontUsedHits_A {"L1Algo::vDontUsedHits_A"};
- L1Vector<L1Hit> vDontUsedHits_B {"L1Algo::vDontUsedHits_B"};
- L1Vector<L1Hit> vDontUsedHits_Buf {"L1Algo::vDontUsedHits_Buf"};
- L1Vector<L1HitPoint> vDontUsedHitsxy_A {"L1Algo::vDontUsedHitsxy_A"};
- L1Vector<L1HitPoint> vDontUsedHitsxy_buf {"L1Algo::vDontUsedHitsxy_buf"};
- L1Vector<L1HitPoint> vDontUsedHitsxy_B {"L1Algo::vDontUsedHitsxy_B"};
+ L1Vector<L1Hit> vNotUsedHits_A {"L1Algo::vNotUsedHits_A"};
+ L1Vector<L1Hit> vNotUsedHits_B {"L1Algo::vNotUsedHits_B"};
+ L1Vector<L1Hit> vNotUsedHits_Buf {"L1Algo::vNotUsedHits_Buf"};
+ L1Vector<L1HitPoint> vNotUsedHitsxy_A {"L1Algo::vNotUsedHitsxy_A"};
+ L1Vector<L1HitPoint> vNotUsedHitsxy_buf {"L1Algo::vNotUsedHitsxy_buf"};
+ L1Vector<L1HitPoint> vNotUsedHitsxy_B {"L1Algo::vNotUsedHitsxy_B"};
L1Vector<L1Track> fTracks_local[L1Constants::size::kMaxNthreads] {"L1Algo::fTracks_local"};
L1Vector<L1HitIndex_t> fRecoHits_local[L1Constants::size::kMaxNthreads] {"L1Algo::fRecoHits_local"};
@@ -275,12 +292,11 @@ public:
L1Vector<L1HitIndex_t> RealIHit_v_buf2 {"L1Algo::RealIHit_v_buf2"};
#ifdef _OPENMP
- L1Vector<omp_lock_t> fHitToBestTrackF {"L1Algo::fHitToBestTrackF"};
- L1Vector<omp_lock_t> fHitToBestTrackB {"L1Algo::fHitToBestTrackB"};
+ L1Vector<omp_lock_t> fStripToTrackLock {"L1Algo::fStripToTrackLock"};
#endif
- L1Vector<int> fStripToTrack {"L1Algo::fStripToTrack"}; // front strip to track pointers
- L1Vector<int> fStripToTrackB {"L1Algo::fStripToTrackB"}; // back strip to track pointers
+ L1Vector<int> fStripToTrackF {"L1Algo::fStripToTrack"}; // strip to track pointers
+ L1Vector<int> fStripToTrackB {"L1Algo::fStripToTrack"}; // strip to track pointers
int fNThreads {0};
bool fUseHitErrors {true};
@@ -327,7 +343,8 @@ public:
L1HitIndex_t HitsUnusedStopIndexEnd[L1Constants::size::kMaxNstations + 1] {0};
- L1Vector<int> TripForHit[2] {"L1Algo::TripForHit"}; // TODO: what does '2' stand for?
+ L1Vector<int> fHitFirstTriplet {"L1Algo::fHitFirstTriplet"}; /// link hit -> first triplet { hit, *, *}
+ L1Vector<int> fHitNtriplets {"L1Algo::fHitNtriplets"}; /// link hit ->n triplets { hit, *, *}
// fvec u_front[Portion/fvecLen], u_back[Portion/fvecLen];
@@ -487,7 +504,7 @@ private:
Tindex start_lh, Tindex n1_l, L1HitPoint* Hits_l,
// output
fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* HitTime_l, fvec* HitTimeEr, fvec* Event_l,
- fvec* d_x, fvec* d_y, fvec* d_xy, fvec* d_u, fvec* d_v);
+ fvec* d_u, fvec* d_v);
/// Get the field approximation. Add the target to parameters estimation. Propagate to middle station.
void findSingletsStep1( // input
@@ -495,7 +512,7 @@ private:
fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr,
// output
- L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_x, fvec* d_y, fvec* d_xy, fvec* d_u, fvec* d_v);
+ L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_u, fvec* d_v);
/// Find the doublets. Reformat data in the portion of doublets.
void findDoubletsStep0( // input
@@ -550,12 +567,7 @@ private:
Tindex n3, int istal, int istam, int istar, nsL1::vector<L1TrackPar>::TSimd& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3,
// output
- Tindex& nstaltriplets, L1Vector<L1HitIndex_t>* hitsn_3 = 0, L1Vector<L1HitIndex_t>* hitsr_5 = 0
-
- // #ifdef XXX
- // ,unsigned int &stat_n_trip
- // #endif
- );
+ Tindex& nstaltriplets);
/// Find neighbours of triplets. Calculate level of triplets.
diff --git a/reco/L1/L1Algo/L1Branch.h b/reco/L1/L1Algo/L1Branch.h
index 3ed84aa7e7e049b6927cabc212ed793bb9e3b8a1..07ff6c13845a2b4ad909279d1ac57d3705675de3 100644
--- a/reco/L1/L1Algo/L1Branch.h
+++ b/reco/L1/L1Algo/L1Branch.h
@@ -23,7 +23,8 @@ struct L1Branch {
char Momentum {0};
char NHits {0};
fscal chi2 {0.};
- int CandIndex {0};
+ int fID {0};
+ bool fIsAlive {0};
L1Vector<L1HitIndex_t> fHits {"L1Branch::fHits"};
// static bool compareCand(const L1Branch *a, const L1Branch *b){
@@ -68,8 +69,7 @@ struct L1Branch {
if (a.ista != b.ista) return (a.ista < b.ista);
- else
- return (a.chi2 < b.chi2);
+ return (a.chi2 < b.chi2);
}
@@ -91,19 +91,21 @@ struct L1Branch {
// return (a.Quality > b.Quality );
//}
- void Set(unsigned char iStation, unsigned char Length, float Chi2, float Qp)
+ void Set(unsigned char iStation, unsigned char Length, float Chi2, float Qp, int ID)
{
NHits = Length;
ista = iStation;
// iN = ( (static_cast<unsigned char>( Chi2 ))<<3 ) + (Level%8);
//unsigned short int ista_l = 16-iStation;
- float tmp = sqrt(Chi2) / 3.5 * 255;
- if (tmp > 255) tmp = 255;
+ // float tmp = sqrt(Chi2) / 3.5 * 255;
+ // if (tmp > 255) tmp = 255;
// unsigned short int chi_2 = 255 - static_cast<unsigned char>( tmp );
// Quality = (Length<<12) + (ista_l<<8) + chi_2;
Momentum = 1.0 / fabs(Qp);
// chi2 = chi_2;
chi2 = Chi2;
+ fID = ID;
+ fIsAlive = true;
}
// void SetLength( unsigned char Length ){
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 39473695abc4f901fb1b10dd591f57d6ceb66895..c991197116230c17866a60d55e1dc40684e5671e 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -78,7 +78,7 @@ inline void L1Algo::findSingletsStep0( // input
// output
fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, L1HitIndex_t* hitsl, fvec* HitTime_l, fvec* HitTimeEr,
// comment unused parameters, FU, 18.01.21
- fvec* /*Event_l*/, fvec* /*d_x*/, fvec* /*d_y*/, fvec* /*d_xy*/, fvec* d_u, fvec* d_v)
+ fvec* /*Event_l*/, fvec* d_u, fvec* d_v)
{
/// Prepare the portion of data of left hits of a triplet:
@@ -125,9 +125,7 @@ inline void L1Algo::findSingletsStep1( /// input 1st stage of singlet search
fvec* u_front_l, fvec* u_back_l, fvec* zPos_l, fvec* HitTime_l, fvec* HitTimeEr,
// output
// L1TrackPar *T_1,
- L1TrackPar* T_1, L1FieldRegion* fld_1,
- // comment unused parameters, FU, 18.01.21
- fvec* /*d_x*/, fvec* /*d_y*/, fvec* /*d_xy*/, fvec* d_u, fvec* d_v)
+ L1TrackPar* T_1, L1FieldRegion* fld_1, fvec* d_u, fvec* d_v)
{
/// Get the field approximation. Add the target to parameters estimation.
@@ -409,7 +407,7 @@ inline void L1Algo::findDoubletsStep0(
if (fParameters.DevIsMatchDoubletsViaMc()) { // trd2d
int indL = HitsUnusedStartIndex[iStaL] + hitsl_1[i1];
int indM = HitsUnusedStartIndex[iStaM] + imh;
- if (GetMcTrackIdForHit(indL) != GetMcTrackIdForHit(indM)) { continue; }
+ if (GetMcTrackIdForUnusedHit(indL) != GetMcTrackIdForUnusedHit(indM)) { continue; }
}
// check y-boundaries
@@ -1079,8 +1077,7 @@ inline void L1Algo::findTripletsStep3( // input
Tindex n3, int istal, int istam, int istar, nsL1::vector<L1TrackPar>::TSimd& T_3, L1Vector<L1HitIndex_t>& hitsl_3,
L1Vector<L1HitIndex_t>& hitsm_3, L1Vector<L1HitIndex_t>& hitsr_3,
// output
- Tindex& nstaltriplets, L1Vector<L1HitIndex_t>* /*hitsn_3*/, L1Vector<L1HitIndex_t>* /*hitsr_5*/
-)
+ Tindex& nstaltriplets)
{
/// Select good triplets and save them into fTriplets. Find neighbouring triplets at the next station.
@@ -1112,8 +1109,8 @@ inline void L1Algo::findTripletsStep3( // input
unsigned int Location = PackTripletId(istal, Thread, fTriplets[istal][Thread].size());
if (ihitl_priv == 0 || ihitl_priv != hitsl_3[i3]) {
- TripForHit[0][ihitl] = Location;
- TripForHit[1][ihitl] = Location;
+ fHitFirstTriplet[ihitl] = Location;
+ fHitNtriplets[ihitl] = 0;
}
ihitl_priv = hitsl_3[i3];
@@ -1153,13 +1150,13 @@ inline void L1Algo::findTripletsStep3( // input
++nstaltriplets;
- TripForHit[1][ihitl] = Location + 1;
+ fHitNtriplets[ihitl]++;
if (istal > (fNstations - 4)) continue;
- unsigned int nNeighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
+ unsigned int nNeighbours = fHitNtriplets[ihitm];
- unsigned int neighLocation = TripForHit[0][ihitm];
+ unsigned int neighLocation = fHitFirstTriplet[ihitm];
unsigned int neighStation = TripletId2Station(neighLocation);
unsigned int neighThread = TripletId2Thread(neighLocation);
unsigned int neighTriplet = TripletId2Triplet(neighLocation);
@@ -1234,9 +1231,9 @@ inline void L1Algo::f5( // input
L1Vector<unsigned int> neighCands("L1CATrackFinder::neighCands"); // save neighbour candidates
neighCands.reserve(8); // ~average is 1-2 for central, up to 5
- unsigned int nNeighbours = TripForHit[1][ihitm] - TripForHit[0][ihitm];
+ unsigned int nNeighbours = fHitNtriplets[ihitm];
- unsigned int neighLocation = TripForHit[0][ihitm];
+ unsigned int neighLocation = fHitFirstTriplet[ihitm];
unsigned int neighStation = TripletId2Station(neighLocation);
unsigned int neighThread = TripletId2Thread(neighLocation);
unsigned int neighTriplet = TripletId2Triplet(neighLocation);
@@ -1308,7 +1305,6 @@ inline void L1Algo::DupletsStaPort(
L1HitPoint* vHits_m = &((*vHitPointsUnused)[0]) + HitsUnusedStartIndex[istam];
fvec u_front[Portion / fvecLen], u_back[Portion / fvecLen];
- fvec dx0[Portion / fvecLen], dy0[Portion / fvecLen], dxy0[Portion / fvecLen];
fvec dv0[Portion / fvecLen], du0[Portion / fvecLen];
fvec zPos[Portion / fvecLen];
fvec HitTime[Portion / fvecLen];
@@ -1321,7 +1317,7 @@ inline void L1Algo::DupletsStaPort(
findSingletsStep0( // input
(ip - portionStopIndex_[istal + 1]) * Portion, n1, vHits_l,
// output
- u_front, u_back, zPos, hitsl_1, HitTime, HitTimeEr, Event, dx0, dy0, dxy0, du0, dv0);
+ u_front, u_back, zPos, hitsl_1, HitTime, HitTimeEr, Event, du0, dv0);
for (Tindex i = 0; i < n1; ++i)
L1_ASSERT(hitsl_1[i] < HitsUnusedStopIndex[istal] - HitsUnusedStartIndex[istal],
@@ -1335,7 +1331,7 @@ inline void L1Algo::DupletsStaPort(
u_front, u_back, zPos, HitTime, HitTimeEr,
// output
- T_1, fld_1, dx0, dy0, dxy0, du0, dv0);
+ T_1, fld_1, du0, dv0);
/// Find the doublets. Reformat data in the portion of doublets.
@@ -1636,11 +1632,11 @@ void L1Algo::CATrackFinder()
RealIHitP = &RealIHit_v;
RealIHitPBuf = &RealIHit_v_buf;
- vHitsUnused = &vDontUsedHits_B; /// array of hits used on current iteration
- L1Vector<L1Hit>* vHitsUnused_buf = &vDontUsedHits_A; /// buffer for copy
+ vHitsUnused = &vNotUsedHits_B; /// array of hits used on current iteration
+ L1Vector<L1Hit>* vHitsUnused_buf = &vNotUsedHits_A; /// buffer for copy
- vHitPointsUnused = &vDontUsedHitsxy_B; /// array of info for hits used on current iteration
- L1Vector<L1HitPoint>* vHitPointsUnused_buf = &vDontUsedHitsxy_A;
+ vHitPointsUnused = &vNotUsedHitsxy_B; /// array of info for hits used on current iteration
+ L1Vector<L1HitPoint>* vHitPointsUnused_buf = &vNotUsedHitsxy_A;
NHitsIsecAll = 0;
fTracks.clear();
@@ -1649,11 +1645,11 @@ void L1Algo::CATrackFinder()
fRecoHits.reserve(2 * vHits->size());
fTracks.reserve(2 * vHits->size() / fNstations);
- int nDontUsedHits = 0;
+ int nNotUsedHits = 0;
- // #pragma omp parallel for reduction(+:nDontUsedHits)
+ // #pragma omp parallel for reduction(+:nNotUsedHits)
for (int ista = 0; ista < fNstations; ++ista) {
- nDontUsedHits += (HitsStopIndex[ista] - HitsStartIndex[ista]);
+ nNotUsedHits += (HitsStopIndex[ista] - HitsStartIndex[ista]);
HitsUnusedStartIndex[ista] = HitsStartIndex[ista];
HitsUnusedStopIndex[ista] = HitsStopIndex[ista];
}
@@ -1675,23 +1671,23 @@ void L1Algo::CATrackFinder()
c_timerG.Start();
#endif // XXX
- float yStep = 1.5 / sqrt(nDontUsedHits); // empirics. 0.01*sqrt(2374) ~= 0.5
+ float yStep = 1.5 / sqrt(nNotUsedHits); // empirics. 0.01*sqrt(2374) ~= 0.5
- // float yStep = 0.5 / sqrt(nDontUsedHits); // empirics. 0.01*sqrt(2374) ~= 0.5
+ // float yStep = 0.5 / sqrt(nNotUsedHits); // empirics. 0.01*sqrt(2374) ~= 0.5
if (yStep > 0.3) yStep = 0.3;
float xStep = yStep * 3;
// float xStep = yStep * 3;
// yStep = 0.0078;
- // const float hitDensity = sqrt( nDontUsedHits );
+ // const float hitDensity = sqrt( nNotUsedHits );
// float yStep = 0.7*4/hitDensity; // empirics. 0.01*sqrt(2374) ~= 0.5
// if (yStep > 0.3)
// yStep = 1.25;
// xStep = 2.05;
- vHitsUnused = &vDontUsedHits_Buf;
+ vHitsUnused = &vNotUsedHits_Buf;
for (int iS = 0; iS < fNstations; ++iS) {
@@ -1700,7 +1696,7 @@ void L1Algo::CATrackFinder()
int start = HitsUnusedStartIndex[iS];
int nhits = HitsUnusedStopIndex[iS] - start;
if (nhits > 0) {
- vGridTime[iS].StoreHits(nhits, &((*vHits)[start]), iS, *this, start, &(vDontUsedHits_Buf[start]),
+ vGridTime[iS].StoreHits(nhits, &((*vHits)[start]), iS, *this, start, &(vNotUsedHits_Buf[start]),
&((*vHits)[start]), &(RealIHit_v[start]));
}
else { // to avoid out-of-range crash in array[start]
@@ -1721,12 +1717,12 @@ void L1Algo::CATrackFinder()
#pragma omp parallel for schedule(dynamic, 5)
#endif
for (L1HitIndex_t ih = HitsStartIndex[ista]; ih < HitsStopIndex[ista]; ++ih) {
- CreateHitPoint(vDontUsedHits_Buf[ih], vDontUsedHitsxy_B[ih]);
+ CreateHitPoint(vNotUsedHits_Buf[ih], vNotUsedHitsxy_B[ih]);
}
}
#ifdef COUNTERS
- stat_nStartHits += nDontUsedHits;
+ stat_nStartHits += nNotUsedHits;
#endif
#ifdef XXX
@@ -1753,7 +1749,7 @@ void L1Algo::CATrackFinder()
#endif
L1ASSERT(0, fNFindIterations == (int) fParameters.GetCAIterations().size());
- isec = 0; // TODO: temporary! (S.Zharko)
+ isec = 0; // TODO: temporary! (S.Zharko)
for (const auto& caIteration : fParameters.GetCAIterations()) // all finder
{
@@ -1786,13 +1782,13 @@ void L1Algo::CATrackFinder()
for (int ist = 0; ist < fNstations; ++ist) {
for (L1HitIndex_t ih = HitsUnusedStartIndex[ist]; ih < HitsUnusedStopIndex[ist]; ++ih) {
//SG!!
- TripForHit[0][ih] = 0;
- TripForHit[1][ih] = 0;
+ fHitFirstTriplet[ih] = 0;
+ fHitNtriplets[ih] = 0;
}
}
/*
- TripForHit[0].assign(HitsUnusedStopIndex[fNstations-1],0);
- TripForHit[1].assign(HitsUnusedStopIndex[fNstations-1],0);
+ fHitFirstTriplet.assign(HitsUnusedStopIndex[fNstations-1],0);
+ fHitNtriplets.assign(HitsUnusedStopIndex[fNstations-1],0);
*/
{
// #pragma omp task
@@ -1802,11 +1798,11 @@ void L1Algo::CATrackFinder()
fFirstCAstation = caIteration.GetFirstStationIndex();
fDoubletChi2Cut = caIteration.GetDoubletChi2Cut(); //11.3449 * 2.f / 3.f; // prob = 0.1
fTripletChi2Cut = caIteration.GetTripletChi2Cut(); //21.1075; // prob = 0.01%
- fPickGather = caIteration.GetPickGather(); //3.0;
- fPickNeighbour = caIteration.GetPickNeighbour(); //5.0;
- fMaxInvMom = caIteration.GetMaxInvMom(); //1.0 / 0.5; // max considered q/p
- fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
- fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
+ fPickGather = caIteration.GetPickGather(); //3.0;
+ fPickNeighbour = caIteration.GetPickNeighbour(); //5.0;
+ fMaxInvMom = caIteration.GetMaxInvMom(); //1.0 / 0.5; // max considered q/p
+ fMaxSlopePV = caIteration.GetMaxSlopePV(); //1.1;
+ fMaxSlope = caIteration.GetMaxSlope(); //2.748; // corresponds to 70 grad
// define the target
fTargX = fParameters.GetTargetPositionX();
@@ -2058,24 +2054,29 @@ void L1Algo::CATrackFinder()
// collect consequtive: the longest tracks, shorter, more shorter and so on
- for (int ilev = fNstations - 3; ilev >= min_level; ilev--) {
- // choose length in triplets number - ilev - the maximum possible triplet level among all triplets in the searched candidate
+ for (int firstTripletLevel = fNstations - 3; firstTripletLevel >= min_level; firstTripletLevel--) {
+ // choose length in triplets number - firstTripletLevel - the maximum possible triplet level among all triplets in the searched candidate
TStopwatch Time;
// how many levels to check
- int nlevel = (fNstations - 2) - ilev + 1;
+ int nlevel = (fNstations - 2) - firstTripletLevel + 1;
- const unsigned char min_best_l = (ilev > min_level) ? ilev + 2 : min_level + 3; // loose maximum
+ const unsigned char min_best_l =
+ (firstTripletLevel > min_level) ? firstTripletLevel + 2 : min_level + 3; // loose maximum
for (int i = 0; i < fNThreads; ++i) {
fTrackCandidates[i].clear();
}
- fStripToTrack.reset(fNstrips, -1);
+ fStripToTrackF.reset(fNstrips, -1);
fStripToTrackB.reset(fNstrips, -1);
- for (int istaF = fFirstCAstation; istaF <= fNstations - 3 - ilev; ++istaF) {
+ //== Loop over triplets with the required level, find and store track candidates
+
+ for (int istaF = fFirstCAstation; istaF <= fNstations - 3 - firstTripletLevel; ++istaF) {
+
+ if (--nlevel == 0) break; //TODO: SG: this is not needed
#ifdef _OPENMP
#pragma omp parallel for firstprivate(curr_tr, new_tr, best_tr, curr_chi2, best_chi2, best_L, curr_L, \
@@ -2115,12 +2116,13 @@ void L1Algo::CATrackFinder()
if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
- if ((ilev == 0) && (GetFStation((*fStripFlag)[(*vHitsUnused)[first_trip.GetLHit()].f]) != 0))
+ if ((firstTripletLevel == 0)
+ && (GetFStation((*fStripFlag)[(*vHitsUnused)[first_trip.GetLHit()].f]) != 0))
continue; // ghost supression // collect only MAPS tracks-triplets CHECK!!!
}
}
- if (first_trip.GetLevel() < ilev)
- continue; // try only triplets, which can start track with ilev+3 length. w\o it have more ghosts, but efficiency either
+ if (first_trip.GetLevel() < firstTripletLevel)
+ continue; // try only triplets, which can start track with firstTripletLevel+3 length. w\o it have more ghosts, but efficiency either
}
@@ -2144,7 +2146,7 @@ void L1Algo::CATrackFinder()
new_tr); /// reqursive func to build a tree of possible track-candidates and choose the best
// if ( best_L < min_best_l ) continue;
- if (best_L < ilev + 2) continue; // lose maximum one hit
+ if (best_L < firstTripletLevel + 2) continue; // lose maximum one hit
if (best_L < min_level + 3) continue; // should find all hits for min_level
@@ -2162,11 +2164,11 @@ void L1Algo::CATrackFinder()
}
}
#endif
- best_tr.Set(istaF, best_L, best_chi2, first_trip.GetQp());
+ best_tr.Set(istaF, best_L, best_chi2, first_trip.GetQp(),
+ PackTrackId(thread_num, fTrackCandidates[thread_num].size()));
fTrackCandidates[thread_num].push_back(best_tr);
L1Branch& tr = fTrackCandidates[thread_num].back();
- tr.CandIndex = fTrackCandidates[thread_num].size() - 1 + thread_num * 10000000;
bool check = 1;
@@ -2174,57 +2176,57 @@ void L1Algo::CATrackFinder()
phitIt != tr.fHits.end(); ++phitIt) {
const L1Hit& h = (*vHits)[*phitIt];
#ifdef _OPENMP
- omp_set_lock(&fHitToBestTrackB[h.b]);
+ omp_set_lock(&fStripToTrackLock[h.b]);
#endif
int& strip1 = (fStripToTrackB)[h.b];
if (strip1 != -1) {
- int thread = strip1 / 10000000;
- int num = (strip1 - thread * 10000000);
+ int thread = TrackId2Thread(strip1);
+ int num = TrackId2Track(strip1);
if (L1Branch::compareCand(tr, fTrackCandidates[thread][num])) {
- fTrackCandidates[thread][num].CandIndex = -1;
- strip1 = tr.CandIndex;
+ fTrackCandidates[thread][num].fID = -1;
+ strip1 = tr.fID;
}
else {
check = 0;
#ifdef _OPENMP
- omp_unset_lock(&fHitToBestTrackB[h.b]);
+ omp_unset_lock(&fStripToTrackLock[h.b]);
#endif
break;
}
}
else
- strip1 = tr.CandIndex;
+ strip1 = tr.fID;
#ifdef _OPENMP
- omp_unset_lock(&fHitToBestTrackB[h.b]);
+ omp_unset_lock(&fStripToTrackLock[h.b]);
#endif
if (check) {
#ifdef _OPENMP
- omp_set_lock(&fHitToBestTrackF[h.f]);
+ omp_set_lock(&fStripToTrackLock[h.f]);
#endif
- int& strip2 = (fStripToTrack)[h.f];
+ int& strip2 = (fStripToTrackF)[h.f];
if (strip2 != -1) {
- int thread = strip2 / 10000000;
- int num = (strip2 - thread * 10000000);
+ int thread = TrackId2Thread(strip2);
+ int num = TrackId2Track(strip2);
if (L1Branch::compareCand(tr, fTrackCandidates[thread][num])) {
- fTrackCandidates[thread][num].CandIndex = -1;
- strip2 = tr.CandIndex;
+ fTrackCandidates[thread][num].fID = -1;
+ strip2 = tr.fID;
}
else {
check = 0;
#ifdef _OPENMP
- omp_unset_lock(&fHitToBestTrackF[h.f]);
+ omp_unset_lock(&fStripToTrackLock[h.f]);
#endif
break;
}
}
else
- strip2 = tr.CandIndex;
+ strip2 = tr.fID;
#ifdef _OPENMP
- omp_unset_lock(&fHitToBestTrackF[h.f]);
+ omp_unset_lock(&fStripToTrackLock[h.f]);
#endif
}
}
@@ -2233,8 +2235,6 @@ void L1Algo::CATrackFinder()
} // iThread
} // istaF
- if (--nlevel == 0) break;
-
for (int i = 0; i < fNThreads; ++i) {
fTracks_local[i].clear();
fRecoHits_local[i].clear();
@@ -2251,11 +2251,11 @@ void L1Algo::CATrackFinder()
bool check = 1;
- if (fTrackCandidates[i][iCandidate].CandIndex != -1) {
+ if (fTrackCandidates[i][iCandidate].fID != -1) {
for (L1Vector<L1HitIndex_t>::iterator phIt = tr.fHits.begin(); /// used strips are marked
phIt != tr.fHits.end(); ++phIt) {
const L1Hit& h = (((*vHits))[*phIt]);
- if (((fStripToTrackB)[h.b] != tr.CandIndex) || ((fStripToTrack)[h.f] != tr.CandIndex)) {
+ if (((fStripToTrackB)[h.b] != tr.fID) || ((fStripToTrackF)[h.f] != tr.fID)) {
check = 0;
break;
}
diff --git a/reco/L1/L1Algo/L1Filtration.h b/reco/L1/L1Algo/L1Filtration.h
index bc5bc5da2ae1744b762b24e2249d8c18ca468bcc..2310f934fcfb94636a63abb953711e971d002148 100644
--- a/reco/L1/L1Algo/L1Filtration.h
+++ b/reco/L1/L1Algo/L1Filtration.h
@@ -346,7 +346,7 @@ inline void L1FilterVtx(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo
inline void L1FilterXY(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo& info)
{
- cnst TWO = 2.;
+ cnst TWO = 2.f;
fvec zeta0, zeta1, S00, S10, S11, si;
fvec F00, F10, F20, F30, F40, F01, F11, F21, F31, F41;
@@ -371,13 +371,13 @@ inline void L1FilterXY(L1TrackPar& T, fvec x, fvec y, const L1XYMeasurementInfo&
S10 = F10 + info.C10;
S11 = F11 + info.C11;
- si = 1. / (S00 * S11 - S10 * S10);
+ si = 1.f / (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.chi2 += zeta0 * zeta0 * S00 + 2.f * zeta0 * zeta1 * S10 + zeta1 * zeta1 * S11;
T.NDF += TWO;
K00 = F00 * S00 + F01 * S10;
@@ -414,7 +414,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[] )
{
diff --git a/reco/L1/L1Algo/L1InitManager.h b/reco/L1/L1Algo/L1InitManager.h
index 86dad6330c4dbc7592635058cb88599e496a9f53..c77951f947db6178f71ccba99b6e778f0685e20c 100644
--- a/reco/L1/L1Algo/L1InitManager.h
+++ b/reco/L1/L1Algo/L1InitManager.h
@@ -222,6 +222,11 @@ public:
///
void SetTrackingLevel(int trackingLevel);
+ /// Sets upper-bound cut on max number of doublets per one singlet
+ void SetMaxDoubletsPerSinglet(unsigned int value) { fParameters.fMaxDoubletsPerSinglet = value; }
+
+ /// Sets upper-bound cut on max number of triplets per one doublet
+ void SetMaxTripletPerDoublets(unsigned int value) { fParameters.fMaxTripletPerDoublets = value; }
/// Transfers L1Parameters object to the destination
/// \param destination Reference to the destination of the L1 object
@@ -240,6 +245,9 @@ public:
/// Flag to match doublets using MC information
void DevSetIsMatchDoubletsViaMc(bool value = true) { fParameters.fDevIsMatchDoubletsViaMc = value; }
+ /// Flag to match triplets using Mc information
+ void DevSetIsMatchTripletsViaMc(bool value = true) { fParameters.fDevIsMatchTripletsViaMc = value; }
+
private:
/// Checker for L1CAIteration container initialization (sets EInitKey::kCAIterations)
/// \return true If all L1CAIteration objects were initialized properly
diff --git a/reco/L1/L1Algo/L1Parameters.cxx b/reco/L1/L1Algo/L1Parameters.cxx
index d11d43fd331b6bc86e9b3deed383fb2fe117035a..7036097f360173d74e0939976a26f0aab74a71cb 100644
--- a/reco/L1/L1Algo/L1Parameters.cxx
+++ b/reco/L1/L1Algo/L1Parameters.cxx
@@ -40,6 +40,7 @@ L1Parameters::L1Parameters(const L1Parameters& other) noexcept
, fDevIsIgnoreHitSearchAreas(other.fDevIsIgnoreHitSearchAreas)
, fDevIsFitSingletsFromTarget(other.fDevIsFitSingletsFromTarget)
, fDevIsMatchDoubletsViaMc(other.fDevIsMatchDoubletsViaMc)
+ , fDevIsMatchTripletsViaMc(other.fDevIsMatchTripletsViaMc)
{
}
@@ -84,6 +85,7 @@ void L1Parameters::Swap(L1Parameters& other) noexcept
std::swap(fDevIsIgnoreHitSearchAreas, other.fDevIsIgnoreHitSearchAreas);
std::swap(fDevIsFitSingletsFromTarget, other.fDevIsFitSingletsFromTarget);
std::swap(fDevIsMatchDoubletsViaMc, other.fDevIsMatchDoubletsViaMc);
+ std::swap(fDevIsMatchTripletsViaMc, other.fDevIsMatchTripletsViaMc);
}
//------------------------------------------------------------------------------------------------------------------------------------
diff --git a/reco/L1/L1Algo/L1Parameters.h b/reco/L1/L1Algo/L1Parameters.h
index d8476c0b16bacd0a0530ba2de73ac1129eb892c9..5b7e36d0ce32296be3ea42838164d34f993c08ac 100644
--- a/reco/L1/L1Algo/L1Parameters.h
+++ b/reco/L1/L1Algo/L1Parameters.h
@@ -177,6 +177,9 @@ public:
/// Flag to match doublets using MC information
bool DevIsMatchDoubletsViaMc() const { return fDevIsMatchDoubletsViaMc; }
+ /// Flag to match triplets using Mc information
+ bool DevIsMatchTripletsViaMc() const { return fDevIsMatchTripletsViaMc; }
+
private:
unsigned int fMaxDoubletsPerSinglet {150}; ///< Upper-bound cut on max number of doublets per one singlet
unsigned int fMaxTripletPerDoublets {15}; ///< Upper-bound cut on max number of triplets per one doublet
@@ -227,6 +230,8 @@ private:
bool fDevIsMatchDoubletsViaMc {false}; ///< Flag to match doublets using MC information
+ bool fDevIsMatchTripletsViaMc {false}; ///< Flag to match triplets using Mc information
+
/********************************
** Friend classes declaration **
********************************/
diff --git a/reco/L1/L1Algo/L1TrackPar.h b/reco/L1/L1Algo/L1TrackPar.h
index 6a97ab80b91914c3fb8ed4086bfb7603521fc9ab..3361d53ae2e4d0adcab69e986217409fe9dc4b74 100644
--- a/reco/L1/L1Algo/L1TrackPar.h
+++ b/reco/L1/L1Algo/L1TrackPar.h
@@ -81,7 +81,14 @@ public:
void SetOneEntry(const int i0, const L1TrackPar& T1, const int i1);
- void Print(int i = -1);
+ void Print(int i = -1) const;
+
+ void PrintCorrelations(int i = -1) const;
+
+ bool IsEntryConsistent(bool printWhenWrong, int i) const;
+
+ bool IsConsistent(bool printWhenWrong, int nFilled) const;
+
// void L1Extrapolate
// (
// // L1TrackPar &T, // input track parameters (x,y,tx,ty,Q/p) and cov.matrix
@@ -95,32 +102,95 @@ public:
// =============================================================================================
-inline void L1TrackPar::Print(int i)
+inline void L1TrackPar::Print(int i) const
{
- std::cout.setf(std::ios::scientific, std::ios::floatfield);
+ // std::cout.setf(std::ios::scientific, std::ios::floatfield);
if (i == -1) {
std::cout << "T = " << std::endl;
- std::cout << x << std::endl;
- std::cout << y << std::endl;
- std::cout << tx << std::endl;
- std::cout << ty << std::endl;
- std::cout << qp << std::endl;
- std::cout << z << std::endl;
+ std::cout << " x " << x << std::endl;
+ std::cout << " y " << y << std::endl;
+ std::cout << " tx " << tx << std::endl;
+ std::cout << " ty " << ty << std::endl;
+ std::cout << " qp " << qp << std::endl;
+ std::cout << " t " << t << std::endl;
+ std::cout << " z " << z << std::endl;
+ std::cout << "C = " << std::endl;
+ std::cout << " c00 " << C00 << std::endl;
+ std::cout << " c11 " << C11 << std::endl;
+ std::cout << " c22 " << C22 << std::endl;
+ std::cout << " c33 " << C33 << std::endl;
+ std::cout << " c44 " << C44 << std::endl;
+ std::cout << " c55 " << C55 << std::endl;
}
else {
std::cout << "T = ";
- std::cout << x[i] << " ";
- std::cout << y[i] << " ";
- std::cout << tx[i] << " ";
- std::cout << ty[i] << " ";
- std::cout << qp[i] << " ";
- std::cout << z[i] << std::endl;
+ std::cout << " x " << x[i];
+ std::cout << " y " << y[i];
+ std::cout << " tx " << tx[i];
+ std::cout << " ty " << ty[i];
+ std::cout << " qp " << qp[i];
+ std::cout << " t " << t[i];
+ std::cout << " z " << z[i] << std::endl;
std::cout << "C = ";
- std::cout << C00[i] << " ";
- std::cout << C11[i] << " ";
- std::cout << C22[i] << " ";
- std::cout << C33[i] << " ";
- std::cout << C44[i] << std::endl;
+ std::cout << " c00 " << C00[i];
+ std::cout << " c11 " << C11[i];
+ std::cout << " c22 " << C22[i];
+ std::cout << " c33 " << C33[i];
+ std::cout << " c44 " << C44[i] << std::endl;
+ std::cout << " c55 " << C55[i] << std::endl;
+ }
+ PrintCorrelations(i);
+}
+
+inline void L1TrackPar::PrintCorrelations(int i) const
+{
+ // std::cout.setf(std::ios::scientific, std::ios::floatfield);
+
+ if (i == -1) {
+ fvec s0 = sqrt(C00);
+ fvec s1 = sqrt(C11);
+ fvec s2 = sqrt(C22);
+ fvec s3 = sqrt(C33);
+ fvec s4 = sqrt(C44);
+ fvec s5 = sqrt(C55);
+
+ std::cout << "K = " << std::endl;
+ std::cout << " k10 " << C10 / s1 / s0 << std::endl;
+
+ std::cout << "\n k20 " << C20 / s2 / s0 << std::endl;
+ std::cout << " k21 " << C21 / s2 / s1 << std::endl;
+
+ std::cout << "\n k30 " << C30 / s3 / s0 << std::endl;
+ std::cout << " k31 " << C31 / s3 / s1 << std::endl;
+ std::cout << " k32 " << C32 / s3 / s2 << std::endl;
+
+ std::cout << "\n k40 " << C40 / s4 / s0 << std::endl;
+ std::cout << " k41 " << C41 / s4 / s1 << std::endl;
+ std::cout << " k42 " << C42 / s4 / s2 << std::endl;
+ std::cout << " k43 " << C43 / s4 / s3 << std::endl;
+
+ std::cout << "\n k50 " << C50 / s5 / s0 << std::endl;
+ std::cout << " k51 " << C51 / s5 / s1 << std::endl;
+ std::cout << " k52 " << C52 / s5 / s2 << std::endl;
+ std::cout << " k53 " << C53 / s5 / s3 << std::endl;
+ std::cout << " k54 " << C54 / s5 / s4 << std::endl;
+ }
+ else {
+ float s0 = sqrt(C00[i]);
+ float s1 = sqrt(C11[i]);
+ float s2 = sqrt(C22[i]);
+ float s3 = sqrt(C33[i]);
+ float s4 = sqrt(C44[i]);
+ float s5 = sqrt(C55[i]);
+
+ std::cout << "K = " << std::endl;
+ std::cout << " " << C10[i] / s1 / s0 << std::endl;
+ std::cout << " " << C20[i] / s2 / s0 << " " << C21[i] / s2 / s1 << std::endl;
+ std::cout << " " << C30[i] / s3 / s0 << " " << C31[i] / s3 / s1 << " " << C32[i] / s3 / s2 << std::endl;
+ std::cout << " " << C40[i] / s4 / s0 << " " << C41[i] / s4 / s1 << " " << C42[i] / s4 / s2 << " "
+ << C43[i] / s4 / s3 << std::endl;
+ std::cout << " " << C50[i] / s5 / s0 << " " << C51[i] / s5 / s1 << " " << C52[i] / s5 / s2 << " "
+ << C53[i] / s5 / s3 << " " << C54[i] / s5 / s4 << std::endl;
}
}
@@ -159,4 +229,71 @@ inline void L1TrackPar::SetOneEntry(const int i0, const L1TrackPar& T1, const in
NDF[i0] = T1.NDF[i1];
} // SetOneEntry
+inline bool L1TrackPar::IsEntryConsistent(bool printWhenWrong, int i) const
+{
+ bool ok = true;
+
+ const fvec* v = &x;
+ const fvec* v1 = &NDF;
+ for (; v < v1; v++) {
+ ok = ok && std::isfinite((*v)[i]);
+ }
+ // verify diagonal elements
+
+ ok = ok && (C00[i] > 0.f);
+ ok = ok && (C11[i] > 0.f);
+ ok = ok && (C22[i] > 0.f);
+ ok = ok && (C33[i] > 0.f);
+ ok = ok && (C44[i] > 0.f);
+ ok = ok && (C55[i] > 0.f);
+
+ // verify non-diagonal elements
+ ok = ok && (C10[i] * C10[i] <= C11[i] * C00[i]);
+
+ ok = ok && (C20[i] * C20[i] <= C22[i] * C00[i]);
+ ok = ok && (C21[i] * C21[i] <= C22[i] * C11[i]);
+
+ ok = ok && (C30[i] * C30[i] <= C33[i] * C00[i]);
+ ok = ok && (C31[i] * C31[i] <= C33[i] * C11[i]);
+ ok = ok && (C32[i] * C32[i] <= C33[i] * C22[i]);
+
+ ok = ok && (C40[i] * C40[i] <= C44[i] * C00[i]);
+ ok = ok && (C41[i] * C41[i] <= C44[i] * C11[i]);
+ ok = ok && (C42[i] * C42[i] <= C44[i] * C22[i]);
+ ok = ok && (C43[i] * C43[i] <= C44[i] * C33[i]);
+
+ ok = ok && (C50[i] * C50[i] <= C55[i] * C00[i]);
+ ok = ok && (C51[i] * C51[i] <= C55[i] * C11[i]);
+ ok = ok && (C52[i] * C52[i] <= C55[i] * C22[i]);
+ ok = ok && (C53[i] * C53[i] <= C55[i] * C33[i]);
+ ok = ok && (C54[i] * C54[i] <= C55[i] * C44[i]);
+
+ if (!ok && printWhenWrong) {
+ std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
+ Print(i);
+ }
+ return ok;
+}
+
+inline bool L1TrackPar::IsConsistent(bool printWhenWrong, int nFilled) const
+{
+ assert(nFilled <= fvecLen);
+ bool ok = true;
+ if (nFilled < 0) { nFilled = fvecLen; }
+ for (int i = 0; i < nFilled; ++i) {
+ ok = ok && IsEntryConsistent(printWhenWrong, i);
+ }
+
+ if (!ok && printWhenWrong) {
+ std::cout << "L1TrackPar parameters are not consistent: " << std::endl;
+ if (nFilled == fvecLen) { std::cout << " All vector elements are filled " << std::endl; }
+ else {
+ std::cout << " Only first " << nFilled << " vector elements are filled " << std::endl;
+ }
+ Print(-1);
+ }
+ return ok;
+}
+
+
#endif
diff --git a/reco/L1/L1Algo/L1Triplet.h b/reco/L1/L1Algo/L1Triplet.h
index 77a90f14150c0f22e0c95cbcad2cca665065e2be..988c06fda59bc8aafa13dc995c8ea490bb48cf9c 100644
--- a/reco/L1/L1Algo/L1Triplet.h
+++ b/reco/L1/L1Algo/L1Triplet.h
@@ -50,8 +50,8 @@ public:
L1HitIndex_t GetMHit() const { return fHitM; }
L1HitIndex_t GetRHit() const { return fHitR; }
- void SetNNeighbours(char n) { fNneighbours = n; }
- char GetNNeighbours() const { return fNneighbours; }
+ void SetNNeighbours(int n) { fNneighbours = n; }
+ int GetNNeighbours() const { return fNneighbours; }
void SetFNeighbour(unsigned int n) { fFirstNeighbour = n; }
unsigned int GetFNeighbour() const { return fFirstNeighbour; }
@@ -86,7 +86,7 @@ private:
L1HitIndex_t fHitL = 0; // left hit index (16b) in vHits array
L1HitIndex_t fHitM = 0; // middle hit index (16b)
L1HitIndex_t fHitR = 0; // right hit index (16b)
- char fNneighbours = 0; // n of neighbouring triplets
+ int fNneighbours = 0; // n of neighbouring triplets
unsigned char fLevel = 0; // triplet level
// == its possible position on the longest track candidate it belongs to.
// level 0 = rightmost triplet of a track candidate
diff --git a/reco/L1/vectors/P4_F32vec4.h b/reco/L1/vectors/P4_F32vec4.h
index aaaca2fe7bb9fc766b13c2def800adad6adfc916..765a33d6872c9f236cdb96226937a64c1c7da922 100644
--- a/reco/L1/vectors/P4_F32vec4.h
+++ b/reco/L1/vectors/P4_F32vec4.h
@@ -239,6 +239,12 @@ public:
/// Checks, if all of the bands is NaN
bool IsNanAll() const { return std::isnan(v[0]) && std::isnan(v[1]) && std::isnan(v[2]) && std::isnan(v[3]); }
+ /// Checks, if all of the bands are finite
+ bool IsFiniteAll() const
+ {
+ return std::isfinite(v[0]) && std::isfinite(v[1]) && std::isfinite(v[2]) && std::isfinite(v[3]);
+ }
+
} __attribute__((aligned(16)));