From a03ca6843a0afee940fd3381c22983bb63cc1f6c Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Tue, 8 Nov 2022 23:28:08 +0000
Subject: [PATCH] L1: remove obsolete code
---
reco/L1/CMakeLists.txt | 7 -
reco/L1/L1Algo/L1Algo.h | 1 -
reco/L1/L1Algo/L1TrackFitter.cxx | 628 ------------------
reco/L1/OffLineInterface/CbmL1MuchFinder.cxx | 300 ---------
reco/L1/OffLineInterface/CbmL1MuchFinder.h | 70 --
.../L1/OffLineInterface/CbmL1MuchFinderQa.cxx | 422 ------------
reco/L1/OffLineInterface/CbmL1MuchFinderQa.h | 72 --
reco/L1/OffLineInterface/CbmL1MuchHit.cxx | 36 -
reco/L1/OffLineInterface/CbmL1MuchHit.h | 34 -
reco/L1/OffLineInterface/CbmL1MuchTrack.cxx | 17 -
reco/L1/OffLineInterface/CbmL1MuchTrack.h | 49 --
reco/L1/OffLineInterface/CbmL1SttHit.cxx | 56 --
reco/L1/OffLineInterface/CbmL1SttHit.h | 36 -
reco/L1/OffLineInterface/CbmL1SttTrack.cxx | 25 -
reco/L1/OffLineInterface/CbmL1SttTrack.h | 50 --
.../OffLineInterface/CbmL1SttTrackFinder.cxx | 411 ------------
.../L1/OffLineInterface/CbmL1SttTrackFinder.h | 70 --
17 files changed, 2284 deletions(-)
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchFinder.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchFinder.h
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchFinderQa.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchFinderQa.h
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchHit.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchHit.h
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchTrack.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1MuchTrack.h
delete mode 100644 reco/L1/OffLineInterface/CbmL1SttHit.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1SttHit.h
delete mode 100644 reco/L1/OffLineInterface/CbmL1SttTrack.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1SttTrack.h
delete mode 100644 reco/L1/OffLineInterface/CbmL1SttTrackFinder.cxx
delete mode 100644 reco/L1/OffLineInterface/CbmL1SttTrackFinder.h
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 8057a67b64..7519440585 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -35,13 +35,6 @@ set(SRCS
OffLineInterface/CbmL1StsTrackFinder.cxx
OffLineInterface/CbmL1GlobalTrackFinder.cxx
OffLineInterface/CbmL1GlobalFindTracksEvents.cxx
- #OffLineInterface / CbmL1MuchFinder.cxx
- #OffLineInterface / CbmL1MuchHit.cxx
- #OffLineInterface / CbmL1MuchTrack.cxx
- #OffLineInterface / CbmL1MuchFinderQa.cxx
- #OffLineInterface / CbmL1SttHit.cxx
- #OffLineInterface / CbmL1SttTrackFinder.cxx
- #OffLineInterface / CbmL1SttTrack.cxx
L1Algo/L1Algo.cxx
L1Algo/L1CATrackFinder.cxx
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index f00d23276e..4079694853 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -379,7 +379,6 @@ public:
void KFTrackFitter_simple(); // version, which use procedured used during the reconstruction
void L1KFTrackFitter(); // version from SIMD-KF benchmark
- void L1KFTrackFitterMuch();
float GetMaxInvMom() const { return fMaxInvMom[0]; }
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index cfeb883c74..60c1b48d72 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -757,634 +757,6 @@ void L1Algo::L1KFTrackFitter()
}
}
-void L1Algo::L1KFTrackFitterMuch()
-{
- // cout << " Start L1 Track Fitter " << endl;
- int start_hit = 0; // for interation in fRecoHits[]
-
- // static L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- // static L1FieldRegion fld _fvecalignment;
- L1FieldValue fldB0, fldB1, fldB2 _fvecalignment;
- L1FieldRegion fld _fvecalignment;
-
-
- L1FieldValue fldB01, fldB11, fldB21 _fvecalignment;
- L1FieldRegion fld1 _fvecalignment;
-
- const int nHits = fParameters.GetNstationsActive();
- int iVec = 0, i = 0;
- int nTracks_SIMD = fvec::size();
- L1TrackPar T; // fitting parametr coresponding to current track
-
- L1TrackParFit T1; // fitting parametr coresponding to current track
- L1TrackPar& tr = T1.fTr;
- T1.SetParticleMass(GetDefaultParticleMass());
-
- L1Fit fit;
- fit.SetParticleMass(GetDefaultParticleMass());
-
- L1Track* t[fvec::size()] {nullptr};
-
- const L1Station* sta = fParameters.GetStations().begin();
- L1Station staFirst, staLast;
-
-
- fvec u[L1Constants::size::kMaxNstations];
- fvec v[L1Constants::size::kMaxNstations];
- fvec du2[L1Constants::size::kMaxNstations];
- fvec dv2[L1Constants::size::kMaxNstations];
-
- fvec x[L1Constants::size::kMaxNstations];
- fvec y[L1Constants::size::kMaxNstations];
- fvec d_xx[L1Constants::size::kMaxNstations];
- fvec d_yy[L1Constants::size::kMaxNstations];
- fvec d_xy[L1Constants::size::kMaxNstations];
-
- fvec z[L1Constants::size::kMaxNstations];
-
- fvec time[L1Constants::size::kMaxNstations];
- fvec dt2[L1Constants::size::kMaxNstations];
-
- fvec x_first;
- fvec y_first;
- fvec d_xx_fst;
- fvec d_yy_fst;
- fvec d_xy_fst;
-
- fvec time_first;
- fvec dt2_fst;
-
- fvec x_last;
- fvec y_last;
- fvec d_xx_lst;
- fvec d_yy_lst;
- fvec d_xy_lst;
-
- fvec time_last;
- fvec dt2_lst;
- fvec dz; /// !!!
-
- fvec Sy[L1Constants::size::kMaxNstations];
- fvec w[L1Constants::size::kMaxNstations];
- fvec y_temp;
- fvec x_temp;
- fvec fldZ0;
- fvec fldZ1;
- fvec fldZ2;
- fvec z_start;
- fvec z_end;
-
-
- L1FieldValue fB[L1Constants::size::kMaxNstations], fB_temp _fvecalignment;
-
- int iSta[L1Constants::size::kMaxNstations]; /// !!!
- fvec ZSta[L1Constants::size::kMaxNstations];
- for (int iHit = 0; iHit < nHits; iHit++) {
- ZSta[iHit] = sta[iHit].z;
- }
-
- unsigned short N_vTracks = fTracks.size();
-
- for (unsigned short itrack = 0; itrack < N_vTracks; itrack += fvec::size()) {
- if (N_vTracks - itrack < static_cast<unsigned short>(fvec::size())) nTracks_SIMD = N_vTracks - itrack;
-
- for (i = 0; i < nTracks_SIMD; i++)
- t[i] = &fTracks[itrack + i]; // current track
-
- // get hits of current track
- for (i = 0; i < nHits; i++) {
- w[i] = ZERO;
- z[i] = ZSta[i];
- }
-
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- for (i = 0; i < fParameters.GetNstationsActive(); i++) {
- d_xx[i][iVec] = 0;
- d_yy[i][iVec] = 0;
- }
- }
-
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- int nHitsTrack = t[iVec]->NHits;
- int nHitsTrackField = 0;
- for (i = 0; i < nHitsTrack; i++) {
- const L1Hit& hit = fInputData.GetHit(fRecoHits[start_hit++]);
- const int ista = hit.iSt;
- if (ista < fNfieldStations) nHitsTrackField++;
- iSta[i] = ista;
- w[ista][iVec] = 1.;
-
- d_xx[i][iVec] = 0;
- d_yy[i][iVec] = 0;
-
- u[ista][iVec] = hit.u;
- v[ista][iVec] = hit.v;
- std::tie(x_temp, y_temp) = sta[ista].ConvUVtoXY<fvec>(u[ista], v[ista]);
- x[ista][iVec] = x_temp[iVec];
- y[ista][iVec] = y_temp[iVec];
- time[ista][iVec] = hit.t;
- dt2[ista][iVec] = hit.dt2;
- du2[ista][iVec] = hit.du2;
- dv2[ista][iVec] = hit.dv2;
- std::tie(d_xx[ista], d_xy[ista], d_yy[ista]) = sta[ista].FormXYCovarianceMatrix(du2[ista], dv2[ista]);
-
- // mom[ista][iVec] = hit.p;
- z[ista][iVec] = hit.z;
- 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] = z[ista][iVec];
- x_first[iVec] = x[ista][iVec];
- y_first[iVec] = y[ista][iVec];
- time_first[iVec] = time[ista][iVec];
- dt2_fst[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];
- }
- else if (i == nHitsTrack - 1) {
- z_end[iVec] = z[ista][iVec];
- x_last[iVec] = x[ista][iVec];
- y_last[iVec] = y[ista][iVec];
- time_last[iVec] = time[ista][iVec];
- dt2_lst[iVec] = dt2[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];
- 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];
- }
- }
-
- // fscal prevZ = z_end[iVec];
- fscal prevZ = z_start[iVec];
- fscal cursy = 0., curSy = 0.;
- // for(i = nHitsTrack - 1; i >= 0; i-- ){
- for (i = 0; i < nHitsTrackField; i++) {
- const int ista = iSta[i];
- const L1Station& st = fParameters.GetStation(ista);
- const fscal curZ = z[ista][iVec];
- fscal dZ = curZ - prevZ;
- fscal By = st.fieldSlice.cy[0][0];
- curSy += dZ * cursy + dZ * dZ * By / 2.;
- cursy += dZ * By;
- Sy[ista][iVec] = curSy;
- prevZ = curZ;
- }
- }
-
- //fit backward
-
- int nHitsField = 0;
-
- for (i = 0; i < nHits; i++)
- if (iSta[i] < fNfieldStations) nHitsField++;
-
- GuessVec(T, x, y, z, Sy, w, nHitsField, &z_start);
-
-
- GuessVec(T1, x, y, z, Sy, w, nHitsField, &z_start);
-
-
- for (int iter = 0; iter < 2; iter++) { // 1.5 iterations
-
- fvec qp0 = T.qp;
-
- fvec qp01 = tr.qp;
-
- i = 0;
-
- FilterFirst(T, x_first, y_first, staFirst);
-
- FilterFirst(T1, x_first, y_first, time_first, dt2_fst, staFirst, d_xx_fst, d_yy_fst, d_xy_fst);
-
- fldZ1 = z[i];
-
- sta[i].fieldSlice.GetFieldValue(T.x, T.y, fldB1);
-
- fldB1.Combine(fB[i], w[i]);
-
- fldZ2 = z[i + 2];
- dz = fldZ2 - fldZ1;
- sta[i].fieldSlice.GetFieldValue(T.x + T.tx * dz, T.y + T.ty * dz, fldB2);
- fldB2.Combine(fB[i + 2], w[i + 2]);
-
- fld.Set(fldB2, fldZ2, fldB1, fldZ1, fldB0, fldZ0);
-
- for (++i; i < nHits; i++) {
- fmask initialised = (z[i] <= z_end) & (z_start < z[i]);
- fvec w1 = iif(initialised, w[i], fvec::Zero());
- fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
-
- fldZ0 = z[i];
- dz = (fldZ1 - fldZ0);
-
-
- if (i < 8) {
-
- sta[i].fieldSlice.GetFieldValue(T.x - T.tx * dz, T.y - T.ty * dz, fldB0);
- fldB0.Combine(fB[i], w[i]);
- fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
-
- L1Extrapolate(T, z[i], qp0, fld, &w1);
-
- T1.Extrapolate(z[i], qp01, fld, w1);
-
- //if (i == fNstationsBeforePipe) {
- //fit.L1AddPipeMaterial(T, qp0, wIn);
- //fit.EnergyLossCorrection(T, fit.PipeRadThick, qp0, fvec(-1.f), wIn);
-
- //T1.AddPipeMaterial(qp01, wIn);
- //T1.EnergyLossCorrection(fit.PipeRadThick, qp01, fvec(-1.f), wIn);
- //}
-
- fldB2 = fldB1;
- fldZ2 = fldZ1;
- fldB1 = fldB0;
- fldZ1 = fldZ0;
-
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(-1.f), wIn);
-
- T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick,
- 1);
- }
- else {
- // L1AddMaterial( T, sta[i].materialInfo, qp0, wIn );
- }
-
- L1UMeasurementInfo info = sta[i].frontInfo;
- info.sigma2 = d_xx[i];
- L1Filter(T, info, u[i], w1);
- T1.Filter(info, u[i], d_xx[i], w1);
-
- info = sta[i].backInfo;
- info.sigma2 = d_yy[i];
-
- L1Filter(T, info, v[i], w1);
- T1.Filter(info, v[i], d_yy[i], w1);
- T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
- }
-
- if (i >= 8) {
-
- fvec z_last = z[i];
-
- // T1.ExtrapolateLine( tr.z + 10, &w1);
- // L1ExtrapolateLine( T, T.z + 10);
-
- fvec d_z = tr.z - z_last;
-
- const fvec st = fvec(10);
-
- fvec nofSteps = (abs(d_z) / 10);
-
- int max_steps = 0;
-
- for (int j = 0; j < 4; j++) {
- nofSteps[j] = int(fabs(d_z[j]) / 10);
- if (max_steps < nofSteps[j]) max_steps = nofSteps[j];
- }
-
-
- fvec nofSteps1 = fvec::Zero();
-
- fvec stepSize = iif(nofSteps < fvec::One(), d_z, st) * w1;
- fvec z_cur = tr.z;
- fvec w2 = w1;
-
-
- for (int iStep = 0; iStep < max_steps + 1; iStep++) {
-
- const fmask maskLastStep = (nofSteps == nofSteps1);
- z_cur = iif(maskLastStep, z_last, tr.z + stepSize);
-
- // fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01));
- // T1.ExtrapolateLine1( z, &w2, v_mc);
-
- T1.ExtrapolateLine(z_cur, &w2);
- nofSteps1 += iif(!maskLastStep, fvec::One(), fvec::Zero());
- w2 = iif(!maskLastStep, w2, fvec::Zero());
-
- // T1.ExtrapolateLine( z_last, &w1);
- // L1ExtrapolateLine( T, z_last);
- // TODO: Unify energy loss corrections (S.Zharko)
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- if (i == 11 || i == 14 || i == 17)
- T1.EnergyLossCorrectionIron(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
- fvec(-1.f), wIn);
- if (i == 8)
- T1.EnergyLossCorrectionCarbon(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
- fvec(-1.f), wIn);
- if (i == 9 || i == 10 || i == 12 || i == 13 || i == 15 || i == 16 || i == 18 || i == 19)
- T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
- fvec(-1.f), wIn);
-
- T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, wIn,
- sta[i].materialInfo.thick / (nofSteps + fvec(1)), 1);
-
- wIn = iif(!maskLastStep, wIn, fvec::Zero());
- }
- else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
- }
- }
-
- // T1.ExtrapolateLine( z_last, &w1);
- //
- L1UMeasurementInfo info = sta[i].frontInfo;
- info.sigma2 = d_xx[i];
- L1Filter(T, info, u[i], w1);
- T1.Filter(info, u[i], d_xx[i], w1);
-
- info = sta[i].backInfo;
- info.sigma2 = d_yy[i];
- L1Filter(T, info, v[i], w1);
- T1.Filter(info, v[i], d_yy[i], w1);
- T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
- }
- }
- // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
-
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
-
- t[iVec]->TLast[0] = tr.x[iVec];
- t[iVec]->TLast[1] = tr.y[iVec];
- t[iVec]->TLast[2] = tr.tx[iVec];
- t[iVec]->TLast[3] = tr.ty[iVec];
- t[iVec]->TLast[4] = tr.qp[iVec];
- t[iVec]->TLast[5] = tr.z[iVec];
- t[iVec]->TLast[6] = tr.t[iVec];
-
- t[iVec]->CLast[0] = tr.C00[iVec];
- t[iVec]->CLast[1] = tr.C10[iVec];
- t[iVec]->CLast[2] = tr.C11[iVec];
- t[iVec]->CLast[3] = tr.C20[iVec];
- t[iVec]->CLast[4] = tr.C21[iVec];
- t[iVec]->CLast[5] = tr.C22[iVec];
- t[iVec]->CLast[6] = tr.C30[iVec];
- t[iVec]->CLast[7] = tr.C31[iVec];
- t[iVec]->CLast[8] = tr.C32[iVec];
- t[iVec]->CLast[9] = tr.C33[iVec];
- t[iVec]->CLast[10] = tr.C40[iVec];
- t[iVec]->CLast[11] = tr.C41[iVec];
- t[iVec]->CLast[12] = tr.C42[iVec];
- t[iVec]->CLast[13] = tr.C43[iVec];
- t[iVec]->CLast[14] = tr.C44[iVec];
- t[iVec]->CLast[15] = tr.C50[iVec];
- t[iVec]->CLast[16] = tr.C51[iVec];
- t[iVec]->CLast[17] = tr.C52[iVec];
- t[iVec]->CLast[18] = tr.C53[iVec];
- t[iVec]->CLast[19] = tr.C54[iVec];
- t[iVec]->CLast[20] = tr.C55[iVec];
-
-
- t[iVec]->chi2 = tr.chi2[iVec];
- t[iVec]->NDF = (int) tr.NDF[iVec];
- }
-
- if (iter == 1) continue; // only 1.5 iterations
- // fit backward
-
- i = nHits - 1; //0;
-
- FilterFirst(T, x_last, y_last, staLast);
-
- FilterFirstL(T1, x_last, y_last, time_last, dt2_lst, staLast, d_xx_lst, d_yy_lst, d_xy_lst);
-
- qp0 = T.qp;
- qp01 = tr.qp;
-
-
- // fvec initialised = fvec(z[i] < z_end) & fvec(z_start <= z[i]);
- // fvec wIn = (ONE & (initialised));
-
- // T1.EnergyLossCorrectionAl( fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(1.f), wIn);
- //
- // T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, ONE, sta[i].materialInfo.thick, 0);
-
- for (--i; i >= 0; i--) {
-
- fmask initialised = (z[i] < z_end) & (z_start <= z[i]);
- fvec w1 = iif(initialised, w[i], fvec::Zero());
- fvec wIn = iif(initialised, fvec::One(), fvec::Zero());
-
- if (i >= fNfieldStations - 1) {
-
- fvec z_last = z[i];
-
- // T1.ExtrapolateLine( tr.z + 10, &w1);
- // L1ExtrapolateLine( T, T.z + 10);
-
-
- fvec d_z = tr.z - z_last;
- const fvec st = fvec(10);
- fvec nofSteps = (abs(d_z) / 10);
-
- int max_steps = 0;
-
- for (size_t j = 0; j < fvec::size(); j++) {
- nofSteps[j] = int(fabs(d_z[j]) / 10); //*w1[i];
- if (max_steps < nofSteps[j]) max_steps = nofSteps[j];
- }
-
- fvec nofSteps1 = fvec(0);
- fvec stepSize = wIn * iif((nofSteps < fvec::One()), d_z, st);
- fvec z_cur = tr.z;
- fvec w2 = wIn;
-
- for (int iStep = 0; iStep < max_steps + 1; iStep++) {
-
- const fmask maskLastStep = (nofSteps == nofSteps1);
- z_cur = iif(maskLastStep, z_last, tr.z - stepSize);
-
- // fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01));
- // T1.ExtrapolateLine1( z_cur, &w2, v_mc);
-
- T1.ExtrapolateLine(z_cur, &w2);
- nofSteps1 += iif(!maskLastStep, fvec::One(), fvec::Zero());
-
- // TODO: Unify the selection of energy loss correction! (S.Zharko)
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- if (i == 11 || i == 14 || i == 17)
- T1.EnergyLossCorrectionIron(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
- fvec(1.f), w2);
- if (i == 8)
- T1.EnergyLossCorrectionCarbon(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
- fvec(1.f), w2);
- if (i == 9 || i == 10 || i == 12 || i == 13 || i == 15 || i == 16 || i == 18)
- T1.EnergyLossCorrectionAl(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + 1), qp01,
- fvec(1.f), w2);
-
- T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y) / (nofSteps + fvec(1)), qp01, w2,
- sta[i].materialInfo.thick / (nofSteps + fvec(1)), 0);
-
- w2 = iif(!maskLastStep, w2, fvec::Zero());
- }
- else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, w2);
- }
- }
-
- //fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01));
- // T1.ExtrapolateLine1( z_last, &wIn, v_mc);
-
- T1.ExtrapolateLine(z_last, &wIn);
- L1ExtrapolateLine(T, z_last);
-
-
- L1UMeasurementInfo info = sta[i].frontInfo;
- info.sigma2 = d_xx[i];
-
- L1Filter(T, info, u[i], w1);
- T1.Filter(info, u[i], d_xx[i], w1);
-
- info = sta[i].backInfo;
- info.sigma2 = d_yy[i];
-
- L1Filter(T, info, v[i], w1);
- T1.Filter(info, v[i], d_yy[i], w1);
- T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
- }
-
- if (i < fNfieldStations - 1) {
-
- if (i == fNfieldStations - 2) { // next-to-last field station
-
- T1.ExtrapolateLine(z[i + 1]);
-
- fldZ1 = z[i + 1];
- sta[i + 1].fieldSlice.GetFieldValue(tr.x, tr.y, fldB1);
- fldB1.Combine(fB[i + 1], w[i + 1]);
-
- fldZ2 = z[i - 1];
- dz = fldZ2 - fldZ1;
-
- sta[i + 1].fieldSlice.GetFieldValue(tr.x + tr.tx * dz, tr.y + tr.ty * dz, fldB2);
- fldB2.Combine(fB[i - 1], w[i - 1]);
- }
-
-
- fldZ0 = z[i];
- dz = (fldZ1 - fldZ0);
- sta[i].fieldSlice.GetFieldValue(tr.x - tr.tx * dz, tr.y - tr.ty * dz, fldB0);
- fldB0.Combine(fB[i], w[i]);
-
- fld.Set(fldB0, fldZ0, fldB1, fldZ1, fldB2, fldZ2);
-
- // fvec v_mc = fabs(1/qp01)/sqrt(mass2+fabs(1/qp01)*fabs(1/qp01));
-
- T1.Extrapolate(z[i], qp01, fld, w1);
- // L1Extrapolate( T, z[i], qp0, fld,&w1 );
-
-
- if constexpr (L1Constants::control::kIfUseRadLengthTable) {
- T1.EnergyLossCorrection(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, fvec(1.f), wIn);
- T1.AddThickMaterial(fParameters.GetMaterialThickness(i, tr.x, tr.y), qp01, wIn, sta[i].materialInfo.thick,
- 0);
- }
- else {
- fit.L1AddMaterial(T, sta[i].materialInfo, qp0, wIn);
- }
-
- L1UMeasurementInfo info = sta[i].frontInfo;
- info.sigma2 = du2[i];
- T1.Filter(info, u[i], info.sigma2, w1);
-
- info = sta[i].backInfo;
- info.sigma2 = dv2[i];
- T1.Filter(info, v[i], info.sigma2, w1);
- T1.FilterTime(time[i], dt2[i], w1, sta[i].timeInfo);
-
-
- fldB2 = fldB1;
- fldZ2 = fldZ1;
- fldB1 = fldB0;
- fldZ1 = fldZ0;
- }
- }
- // fit.L1AddHalfMaterial( T, sta[i].materialInfo, qp0 );
-
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
-
- t[iVec]->TFirst[0] = tr.x[iVec];
- t[iVec]->TFirst[1] = tr.y[iVec];
- t[iVec]->TFirst[2] = tr.tx[iVec];
- t[iVec]->TFirst[3] = tr.ty[iVec];
- t[iVec]->TFirst[4] = tr.qp[iVec];
- t[iVec]->TFirst[5] = tr.z[iVec];
- t[iVec]->TFirst[6] = tr.t[iVec];
-
- t[iVec]->CFirst[0] = tr.C00[iVec];
- t[iVec]->CFirst[1] = tr.C10[iVec];
- t[iVec]->CFirst[2] = tr.C11[iVec];
- t[iVec]->CFirst[3] = tr.C20[iVec];
- t[iVec]->CFirst[4] = tr.C21[iVec];
- t[iVec]->CFirst[5] = tr.C22[iVec];
- t[iVec]->CFirst[6] = tr.C30[iVec];
- t[iVec]->CFirst[7] = tr.C31[iVec];
- t[iVec]->CFirst[8] = tr.C32[iVec];
- t[iVec]->CFirst[9] = tr.C33[iVec];
- t[iVec]->CFirst[10] = tr.C40[iVec];
- t[iVec]->CFirst[11] = tr.C41[iVec];
- t[iVec]->CFirst[12] = tr.C42[iVec];
- t[iVec]->CFirst[13] = tr.C43[iVec];
- t[iVec]->CFirst[14] = tr.C44[iVec];
- t[iVec]->CFirst[15] = tr.C50[iVec];
- t[iVec]->CFirst[16] = tr.C51[iVec];
- t[iVec]->CFirst[17] = tr.C52[iVec];
- t[iVec]->CFirst[18] = tr.C53[iVec];
- t[iVec]->CFirst[19] = tr.C54[iVec];
- t[iVec]->CFirst[20] = tr.C55[iVec];
-
- t[iVec]->chi2 = tr.chi2[iVec];
- t[iVec]->NDF = (int) tr.NDF[iVec];
- }
-
- // extrapolate to the PV region
- L1TrackParFit T1pv = T1;
- T1pv.Extrapolate(fParameters.GetTargetPositionZ(), T1pv.fTr.qp, fld, fvec(1.));
- for (iVec = 0; iVec < nTracks_SIMD; iVec++) {
- t[iVec]->Tpv[0] = T1pv.fTr.x[iVec];
- t[iVec]->Tpv[1] = T1pv.fTr.y[iVec];
- t[iVec]->Tpv[2] = T1pv.fTr.tx[iVec];
- t[iVec]->Tpv[3] = T1pv.fTr.ty[iVec];
- t[iVec]->Tpv[4] = T1pv.fTr.qp[iVec];
- t[iVec]->Tpv[5] = T1pv.fTr.z[iVec];
- t[iVec]->Tpv[6] = T1pv.fTr.t[iVec];
-
- t[iVec]->Cpv[0] = T1pv.fTr.C00[iVec];
- t[iVec]->Cpv[1] = T1pv.fTr.C10[iVec];
- t[iVec]->Cpv[2] = T1pv.fTr.C11[iVec];
- t[iVec]->Cpv[3] = T1pv.fTr.C20[iVec];
- t[iVec]->Cpv[4] = T1pv.fTr.C21[iVec];
- t[iVec]->Cpv[5] = T1pv.fTr.C22[iVec];
- t[iVec]->Cpv[6] = T1pv.fTr.C30[iVec];
- t[iVec]->Cpv[7] = T1pv.fTr.C31[iVec];
- t[iVec]->Cpv[8] = T1pv.fTr.C32[iVec];
- t[iVec]->Cpv[9] = T1pv.fTr.C33[iVec];
- t[iVec]->Cpv[10] = T1pv.fTr.C40[iVec];
- t[iVec]->Cpv[11] = T1pv.fTr.C41[iVec];
- t[iVec]->Cpv[12] = T1pv.fTr.C42[iVec];
- t[iVec]->Cpv[13] = T1pv.fTr.C43[iVec];
- t[iVec]->Cpv[14] = T1pv.fTr.C44[iVec];
- t[iVec]->Cpv[15] = T1pv.fTr.C50[iVec];
- t[iVec]->Cpv[16] = T1pv.fTr.C51[iVec];
- t[iVec]->Cpv[17] = T1pv.fTr.C52[iVec];
- t[iVec]->Cpv[18] = T1pv.fTr.C53[iVec];
- t[iVec]->Cpv[19] = T1pv.fTr.C54[iVec];
- t[iVec]->Cpv[20] = T1pv.fTr.C55[iVec];
- }
- } // iter
- }
-}
-
void L1Algo::GuessVec(L1TrackPar& t, fvec* xV, fvec* yV, fvec* zV, fvec* Sy, fvec* wV, int NHits, fvec* zCur)
// gives nice initial approximation for x,y,tx,ty - almost same as KF fit. qp - is shifted by 4%, resid_ual - ~3.5% (KF fit resid_ual - 1%).
diff --git a/reco/L1/OffLineInterface/CbmL1MuchFinder.cxx b/reco/L1/OffLineInterface/CbmL1MuchFinder.cxx
deleted file mode 100644
index c47706649e..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchFinder.cxx
+++ /dev/null
@@ -1,300 +0,0 @@
-/* Copyright (C) 2006-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov, Denis Bertini [committer] */
-
-#include "CbmL1MuchFinder.h"
-
-#include "CbmKF.h"
-#include "CbmKFHit.h"
-#include "CbmKFMaterial.h"
-#include "CbmKFMath.h"
-#include "CbmKFPixelMeasurement.h"
-#include "CbmKFTrackInterface.h"
-#include "CbmL1MuchHit.h"
-#include "CbmL1MuchTrack.h"
-#include "CbmMCTrack.h"
-#include "CbmMuchHit.h"
-#include "CbmMuchPoint.h"
-#include "CbmMuchTrack.h"
-#include "CbmStsKFTrackFitter.h"
-#include "CbmStsTrack.h"
-#include "CbmStsTrackMatch.h"
-#include "CbmVertex.h"
-
-#include "FairRootManager.h"
-
-#include "TClonesArray.h"
-#include "TFile.h"
-#include "TLorentzVector.h"
-#include "TVector3.h"
-
-#include <iostream>
-#include <vector>
-
-#include <cmath>
-
-using std::cout;
-using std::endl;
-using std::fabs;
-using std::vector;
-
-ClassImp(CbmL1MuchFinder);
-
-CbmL1MuchFinder::CbmL1MuchFinder(const char* name, Int_t iVerbose) : FairTask(name, iVerbose)
-{
- fTrackCollection = new TClonesArray("CbmMuchTrack", 100);
- histodir = 0;
-}
-
-
-CbmL1MuchFinder::~CbmL1MuchFinder() {}
-
-InitStatus CbmL1MuchFinder::Init() { return ReInit(); }
-
-InitStatus CbmL1MuchFinder::ReInit()
-{
- fMuchPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("MuchPoint");
- fMuchHits = (TClonesArray*) FairRootManager::Instance()->GetObject("MuchHit");
- fStsTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("StsTrack");
- fMCTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
- fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("StsTrackMatch");
- //fPrimVtx = (CbmVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
- // Get pointer to PrimaryVertex object from IOManager if it exists
- // The old name for the object is "PrimaryVertex" the new one
- // "PrimaryVertex." Check first for the new name
- fPrimVtx = dynamic_cast<CbmVertex*>(FairRootManager::Instance()->GetObject("PrimaryVertex."));
- if (nullptr == fPrimVtx) {
- fPrimVtx = dynamic_cast<CbmVertex*>(FairRootManager::Instance()->GetObject("PrimaryVertex"));
- }
- if (nullptr == fPrimVtx) {
- Error("CbmL1MuchFinder::ReInit", "vertex not found!");
- return kERROR;
- }
-
- fStsFitter.Init();
-
- FairRootManager::Instance()->Register("MuchTrack", "Much", fTrackCollection, IsOutputBranchPersistent("MuchTrack"));
-
- return kSUCCESS;
-}
-
-void CbmL1MuchFinder::SetParContainers() {}
-
-void CbmL1MuchFinder::Finish() { Write(); }
-
-void CbmL1MuchFinder::Exec(Option_t* /*option*/)
-{
- const int MaxBranches = 50;
-
- static bool first_call_murec = 1;
- fTrackCollection->Clear();
- static int EventCounter = 0;
- EventCounter++;
- cout << " MuRec event " << EventCounter << endl;
-
- int MuNStations = CbmKF::Instance()->MuchStation2MCIDMap.size();
-
- if (first_call_murec) {
- first_call_murec = 0;
- TDirectory* curdir = gDirectory;
- histodir = gDirectory->mkdir("MuRec");
- histodir->cd();
- fhNBranches = new TH1F("NBranches", "N Branches", MaxBranches, 0, MaxBranches);
- curdir->cd();
- }
-
- int NHits = fMuchHits->GetEntriesFast();
- vector<CbmL1MuchHit> vMuchHits;
-
- for (int ih = 0; ih < NHits; ih++) {
- CbmMuchHit* h = (CbmMuchHit*) fMuchHits->At(ih);
- CbmL1MuchHit m(h, ih);
- vMuchHits.push_back(m);
- }
-
- vector<CbmL1MuchTrack> vTracks;
-
- int NStsTracks = fStsTracks->GetEntriesFast();
-
- CbmL1MuchTrack Branches[MaxBranches];
-
- for (int itr = 0; itr < NStsTracks; itr++) {
-
- CbmStsTrack* sts = (CbmStsTrack*) fStsTracks->At(itr);
- if (sts->GetNStsHits() + sts->GetNMvdHits() < 4) continue;
-
- // MC
- if (0 && fSTSTrackMatch && fMCTracks) {
- CbmStsTrackMatch* m = (CbmStsTrackMatch*) fSTSTrackMatch->At(itr);
- if (!m) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- if (mcTrackID < 0) continue;
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if (!mcTrack) continue;
- if (abs(mcTrack->GetPdgCode()) != 13) continue;
- }
-
- fStsFitter.DoFit(sts, 13); // refit with muon mass
-
- int NBranches = 1;
-
- Branches[0].SetStsTrack(sts);
- Branches[0].StsID = itr;
- Branches[0].NHits = 0;
- Branches[0].NMissed = 0;
- Branches[0].NMissedStations = 0;
- Branches[0].ok = 1;
- Branches[0].stopped = 0;
- Branches[0].vHits.clear();
- //cout<<"Sts track N "<<itr<<" with initial mom="<<1./Branches[0].T[4]<<endl;
- for (int ist = 0; ist < MuNStations; ist++) {
-
- int NBranchesOld = NBranches;
-
- for (int ibr = 0; ibr < NBranchesOld; ibr++) {
- CbmL1MuchTrack& tr = Branches[ibr];
- if (tr.stopped) continue;
- //if( ist%3 ==0 ) cout<<" | ";
- double Zdet = CbmKF::Instance()->vMuchDetectors[ist].ZReference;
- tr.Extrapolate(Zdet);
- if (fabs(tr.T[4]) > 100.) tr.stopped = 1;
- if (1. < 0.5 * fabs(tr.T[4])) tr.stopped = 1; // 0.5 GeV, stop transport
- //if( tr.stopped ) cout<<"Sts track N "<<itr<<" stopped at Mu station "<<ist
- //<<" with mom="<<1./tr.T[4]<<endl;
- if (tr.stopped) continue;
- if (CbmKF::Instance()->vMuchDetectors[ist].IsOutside(tr.T[0], tr.T[1])) {
- //cout<<" out ";
- tr.NMissedStations++;
- continue;
- }
-
- vector<int> NewHits;
- for (int ih = 0; ih < NHits; ih++) {
- CbmL1MuchHit& h = vMuchHits[ih];
- if (h.iStation != ist) continue;
- if (0) { // !!! Cut for the time of flight (ns)
- double hl = sqrt(h.FitPoint.x * h.FitPoint.x + h.FitPoint.y * h.FitPoint.y + h.FitPoint.z * h.FitPoint.z);
- double hp = 1. / fabs(tr.T[4]);
- double texp = hl * sqrt(1. - 0.1057 * 0.1057 / (hp * hp)) / 29.9792458; //ns
- if (h.time - texp > 1.0) continue;
- }
- double dx = tr.T[0] - h.FitPoint.x;
- double dy = tr.T[1] - h.FitPoint.y;
- double c0 = tr.C[0] + h.FitPoint.V[0];
- double c1 = tr.C[1] + h.FitPoint.V[1];
- double c2 = tr.C[2] + h.FitPoint.V[2];
- double chi2 = 0.5 * (dx * dx * c0 - 2 * dx * dy * c1 + dy * dy * c2) / (c0 * c2 - c1 * c1);
- if (chi2 <= 20.) NewHits.push_back(ih);
- }
- int nnew = NewHits.size();
- for (int ih = 1; ih < nnew; ih++) {
- if (NBranches >= MaxBranches) break;
- CbmL1MuchTrack& t = Branches[NBranches++];
- t = tr;
- CbmL1MuchHit& h = vMuchHits[NewHits[ih]];
- t.vHits.push_back(&h);
- t.NHits++;
- double qp0 = t.T[4];
- h.Filter(t, 1, qp0);
- }
- if (nnew > 0) {
- CbmL1MuchHit& h = vMuchHits[NewHits[0]];
- tr.vHits.push_back(&h);
- tr.NHits++;
- double qp0 = tr.T[4];
- h.Filter(tr, 1, qp0);
- }
- else
- tr.NMissed++;
- }
- }
- int bestbr = 0;
- for (int ibr = 1; ibr < NBranches; ibr++) {
- if ((Branches[ibr].NHits > Branches[bestbr].NHits)
- || (Branches[ibr].NHits == Branches[bestbr].NHits) && (Branches[ibr].chi2 < Branches[bestbr].chi2))
- bestbr = ibr;
- }
- vTracks.push_back(Branches[bestbr]);
- // MC
- if (fSTSTrackMatch && fMCTracks) {
- CbmStsTrackMatch* m = (CbmStsTrackMatch*) fSTSTrackMatch->At(itr);
- if (!m) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- if (mcTrackID < 0) continue;
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if (!mcTrack) continue;
- if (abs(mcTrack->GetPdgCode()) == 13) fhNBranches->Fill(NBranches);
- }
- } // tracks
-
- int NTracks = vTracks.size();
- cout << "NTracks=" << NTracks << endl;
- //sort tracks and check for common muon hits
-
- vector<CbmL1MuchTrack*> vpTracks;
- for (int i = 0; i < NTracks; i++)
- vpTracks.push_back(&(vTracks[i]));
- sort(vpTracks.begin(), vpTracks.end(), CbmL1MuchTrack::Compare);
-
- int NOutTracks = fTrackCollection->GetEntriesFast();
-
- for (int it = 0; it < NTracks; it++) {
- CbmL1MuchTrack& tr = *vpTracks[it];
- if (tr.NDF <= 0 || tr.chi2 > 10. * tr.NDF) {
- tr.ok = 0;
- continue;
- }
- int nbusy = 0;
- for (int ih = 0; ih < tr.NHits; ih++)
- if (tr.vHits[ih]->busy) nbusy++;
- if (0 && nbusy > 2) {
- tr.ok = 0;
- continue;
- }
- {
- new ((*fTrackCollection)[NOutTracks]) CbmMuchTrack();
- CbmMuchTrack* MuchTrack = (CbmMuchTrack*) fTrackCollection->At(NOutTracks++);
- MuchTrack->SetChi2(tr.GetRefChi2());
- MuchTrack->SetNDF(tr.GetRefNDF());
- FairTrackParam tp;
- CbmKFMath::CopyTC2TrackParam(&tp, tr.T, tr.C);
- MuchTrack->SetMuchTrack(&tp);
- MuchTrack->SetStsTrackID(tr.StsID);
- int nh = 0;
- for (vector<CbmL1MuchHit*>::iterator i = tr.vHits.begin(); i != tr.vHits.end(); i++) {
- if (++nh > 30) break;
- MuchTrack->AddHitIndex((*i)->index);
- }
- MuchTrack->SetNMissedHits(tr.NMissed);
- MuchTrack->SetNMissedStations(tr.NMissedStations);
- }
- for (int ih = 0; ih < tr.NHits; ih++)
- tr.vHits[ih]->busy = 1;
- }
-
- if (EventCounter < 100 && EventCounter % 10 == 0 || EventCounter >= 100 && EventCounter % 100 == 0) Write();
- //cout<<"end of MuRec"<<endl;
-}
-
-
-void CbmL1MuchFinder::Write()
-{
- TFile HistoFile("MuRec.root", "RECREATE");
- writedir2current(histodir);
- HistoFile.Close();
-}
-
-void CbmL1MuchFinder::writedir2current(TObject* obj)
-{
- if (!obj->IsFolder()) obj->Write();
- else {
- TDirectory* cur = gDirectory;
- TDirectory* sub = cur->mkdir(obj->GetName());
- sub->cd();
- TList* listSub = ((TDirectory*) obj)->GetList();
- TIter it(listSub);
- while (TObject* obj1 = it())
- writedir2current(obj1);
- cur->cd();
- }
-}
diff --git a/reco/L1/OffLineInterface/CbmL1MuchFinder.h b/reco/L1/OffLineInterface/CbmL1MuchFinder.h
deleted file mode 100644
index f2c93fcf22..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchFinder.h
+++ /dev/null
@@ -1,70 +0,0 @@
-/* Copyright (C) 2006-2009 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Denis Bertini [committer], Sergey Gorbunov */
-
-#ifndef CBM_L1_MuchFinder_h
-#define CBM_L1_MuchFinder_h
-
-#include "CbmStsKFTrackFitter.h"
-
-#include "FairTask.h"
-
-#include "TH1.h"
-#include "TLorentzVector.h"
-
-
-class TClonesArray;
-
-class CbmL1MuchFinder : public FairTask {
-public:
- /** Constructor **/
- CbmL1MuchFinder(const char* name = "CbmL1MuchFinder", Int_t iVerbose = 1);
-
- /** Destructor **/
- ~CbmL1MuchFinder();
-
- /// * FairTask methods
-
- /** Intialisation at begin of run. To be implemented in the derived class.
- *@value Success If not kSUCCESS, task will be set inactive.
- **/
- InitStatus Init();
-
- /** Reinitialisation.
- *@value Success If not kSUCCESS, task will be set inactive.
- **/
- InitStatus ReInit();
-
- /** Intialise parameter containers.
- **/
- void SetParContainers();
-
- void Exec(Option_t* option);
-
- /** Action after each event. **/
- void Finish();
-
-private:
- TClonesArray* fMuchPoints; //! Much MC points
- TClonesArray* fMuchHits; //! Much Hits
- TClonesArray* fStsTracks;
- TClonesArray* fMCTracks;
- TClonesArray* fSTSTrackMatch;
- TClonesArray* fTrackCollection; //! Much tracks
-
- CbmVertex* fPrimVtx;
- CbmStsKFTrackFitter fStsFitter;
-
- TDirectory* histodir;
-
- void Write();
- void writedir2current(TObject* obj);
-
- TH1F* fhNBranches;
-
-public:
- ClassDef(CbmL1MuchFinder, 1);
-};
-
-
-#endif
diff --git a/reco/L1/OffLineInterface/CbmL1MuchFinderQa.cxx b/reco/L1/OffLineInterface/CbmL1MuchFinderQa.cxx
deleted file mode 100644
index 1050eeb015..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchFinderQa.cxx
+++ /dev/null
@@ -1,422 +0,0 @@
-/* Copyright (C) 2007-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#include "CbmL1MuchFinderQa.h"
-
-#include "CbmKF.h"
-#include "CbmKFHit.h"
-#include "CbmKFMaterial.h"
-#include "CbmKFMath.h"
-#include "CbmKFPixelMeasurement.h"
-#include "CbmKFTrackInterface.h"
-#include "CbmL1MuchHit.h"
-#include "CbmL1MuchTrack.h"
-#include "CbmMCTrack.h"
-#include "CbmMuchHit.h"
-#include "CbmMuchPoint.h"
-#include "CbmMuchTrack.h"
-#include "CbmStsKFTrackFitter.h"
-#include "CbmStsTrack.h"
-#include "CbmStsTrackMatch.h"
-#include "CbmVertex.h"
-
-#include "FairRootManager.h"
-
-#include "TClonesArray.h"
-#include "TFile.h"
-#include "TLorentzVector.h"
-#include "TProfile.h"
-#include "TVector3.h"
-
-#include <iostream>
-#include <vector>
-
-#include <cmath>
-
-using std::cout;
-using std::endl;
-using std::fabs;
-using std::vector;
-
-ClassImp(CbmL1MuchFinderQa);
-
-CbmL1MuchFinderQa::CbmL1MuchFinderQa(const char* name, Int_t iVerbose) : FairTask(name, iVerbose)
-{
- histodir = 0;
- fhPerfAll = 0;
- fhPerfSignal = 0;
-}
-
-
-CbmL1MuchFinderQa::~CbmL1MuchFinderQa() {}
-
-InitStatus CbmL1MuchFinderQa::Init() { return ReInit(); }
-
-InitStatus CbmL1MuchFinderQa::ReInit()
-{
- fMuchPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("MuchPoint");
- fMuchHits = (TClonesArray*) FairRootManager::Instance()->GetObject("MuchHit");
- fStsTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("StsTrack");
- fMuchTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MuchTrack");
- fMCTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
- //fPrimVtx = (CbmVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
- // Get pointer to PrimaryVertex object from IOManager if it exists
- // The old name for the object is "PrimaryVertex" the new one
- // "PrimaryVertex." Check first for the new name
- fPrimVtx = dynamic_cast<CbmVertex*>(FairRootManager::Instance()->GetObject("PrimaryVertex."));
- if (nullptr == fPrimVtx) {
- fPrimVtx = dynamic_cast<CbmVertex*>(FairRootManager::Instance()->GetObject("PrimaryVertex"));
- }
- if (nullptr == fPrimVtx) {
- Error("CbmL1MuchFinderQa::ReInit", "vertex not found!");
- return kERROR;
- }
-
- fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("StsTrackMatch");
- fStsFitter.Init();
-
- return kSUCCESS;
-}
-
-void CbmL1MuchFinderQa::SetParContainers() {}
-
-void CbmL1MuchFinderQa::Finish() { Write(); }
-
-void CbmL1MuchFinderQa::Exec(Option_t* /*option*/)
-{
-
- static int EventCounter = 0;
- EventCounter++;
- cout << " MuRecQa event " << EventCounter << endl;
-
- static bool first_call = 1;
-
- int MuNStations = CbmKF::Instance()->MuchStation2MCIDMap.size();
-
- if (first_call) {
- first_call = 0;
- TDirectory* curdir = gDirectory;
- histodir = gDirectory->mkdir("MuRecQa");
- histodir->cd();
-
- fhPerfSignal = new TProfile("eff_signal", "Signal Mu eff vs mom", 60, 0, 30);
- fhPerfAll = new TProfile("eff_all", "All Mu eff vs mom", 60, 0, 30);
- fhGhost = new TProfile("ghost", "ghost vs mom", 60, 0, 30);
-
- char S1[255], S2[255];
- for (int ist = 0; ist < MuNStations; ist++) {
- sprintf(S1, "station %i", ist);
- //TDirectory *d = histodir->mkdir(S1);
- //d->cd();
- sprintf(S1, "Pull_x_%i", ist);
- sprintf(S2, "Pull_x of sts track at Mu detector %i ", ist);
- histPull_dx[ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "Pull_y_%i", ist);
- sprintf(S2, "Pull_y of sts track at Mu detector %i ", ist);
- histPull_dy[ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "Pull_tx_%i", ist);
- sprintf(S2, "Pull_tx of sts track at Mu detector %i ", ist);
- histPull_tx[ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "Pull_ty_%i", ist);
- sprintf(S2, "Pull_ty of sts track at Mu detector %i ", ist);
- histPull_ty[ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "Pull_qp_%i", ist);
- sprintf(S2, "Pull_qp of sts track at Mu detector %i ", ist);
- histPull_qp[ist] = new TH1F(S1, S2, 100, -10., 10.);
-
- sprintf(S1, "chi2hit_%i", ist);
- sprintf(S2, "chi2 between track and hit at Mu detector %i ", ist);
- histChi2[ist] = new TH1F(S1, S2, 1000, 0, 100);
-
-
- sprintf(S1, "b_Pull_x_%i", ist);
- sprintf(S2, "b_Pull_x of sts track at Mu detector %i ", ist);
- histPull_dx[20 + ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "b_Pull_y_%i", ist);
- sprintf(S2, "b_Pull_y of sts track at Mu detector %i ", ist);
- histPull_dy[20 + ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "b_Pull_tx_%i", ist);
- sprintf(S2, "b_Pull_tx of sts track at Mu detector %i ", ist);
- histPull_tx[20 + ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "b_Pull_ty_%i", ist);
- sprintf(S2, "b_Pull_ty of sts track at Mu detector %i ", ist);
- histPull_ty[20 + ist] = new TH1F(S1, S2, 100, -10., 10.);
- sprintf(S1, "b_Pull_qp_%i", ist);
- sprintf(S2, "b_Pull_qp of sts track at Mu detector %i ", ist);
- histPull_qp[20 + ist] = new TH1F(S1, S2, 100, -10., 10.);
- }
- curdir->cd();
- }
-
- if (!fSTSTrackMatch) return;
- if (!fMuchHits) return;
- if (!fMuchPoints) return;
- if (!fStsTracks) return;
- if (!fMuchTracks) return;
- if (!fSTSTrackMatch) return;
- if (!fMCTracks) return;
-
- int NHits = fMuchHits->GetEntriesFast();
- int NStsTracks = fStsTracks->GetEntriesFast();
- int NMuchTracks = fMuchTracks->GetEntriesFast();
- //int NMCTracks = fMCTracks->GetEntriesFast();
-
- bool was_problem = 0;
-
- static int NMu5 = 0, NMuRec5 = 0, NMu10 = 0, NMuRec10 = 0;
- static int NMuS5 = 0, NMuSRec5 = 0, NMuS10 = 0, NMuSRec10 = 0;
- static int NGhost5 = 0, NGhostRec5 = 0, NGhost10 = 0, NGhostRec10 = 0;
-
- int IndMuTrack[NStsTracks];
-
- for (int itr = 0; itr < NStsTracks; itr++)
- IndMuTrack[itr] = -1;
-
- for (int itr = 0; itr < NMuchTracks; itr++) {
- CbmMuchTrack* Tmu = (CbmMuchTrack*) fMuchTracks->At(itr);
- if (!Tmu) continue;
- int ists = Tmu->GetStsTrackID();
- if (ists < 0) continue;
- IndMuTrack[ists] = itr;
- }
-
- for (int ists = 0; ists < NStsTracks; ists++) {
- CbmStsTrack* Tsts = (CbmStsTrack*) fStsTracks->At(ists);
- CbmStsTrackMatch* m = (CbmStsTrackMatch*) fSTSTrackMatch->At(ists);
- if (!Tsts) continue;
- if (!m) continue;
- double ratio = 1. * m->GetNofTrueHits() / (m->GetNofTrueHits() + m->GetNofWrongHits() + m->GetNofFakeHits());
- if (ratio < .7) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- if (mcTrackID < 0) continue;
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if (!mcTrack) continue;
- //if (fabs((mcTrack->GetStartVertex()).z())>1.) continue;
- //if( mcTrack->GetMotherID() !=-1 ) continue;
-
- Double_t P = mcTrack->GetP();
-
- int NMuHits = 0;
- int NMissedStations = 0;
- int NH = fMuchHits->GetEntriesFast();
- int ista_last = -1;
- for (int ih = 0; ih < NH; ih++) {
- CbmMuchHit* h = (CbmMuchHit*) fMuchHits->At(ih);
- int ist = h->GetStationNr() - 1;
- if (ist != ista_last + 1) { NMissedStations += ist - 1 - ista_last; }
- ista_last = ist;
- int ip = h->GetRefIndex();
- if (ip < 0) continue;
- CbmMuchPoint* p = (CbmMuchPoint*) fMuchPoints->At(ip);
- if (mcTrackID == p->GetTrackID()) NMuHits++;
- }
-
- if (Tsts->GetNStsHits() + Tsts->GetNMvdHits() < 4) continue;
-
- bool Muon = (abs(mcTrack->GetPdgCode()) == 13);
- bool Signal = (mcTrack->GetMotherId() < 0 && mcTrackID < 3 && Muon);
-
- if (Muon && (NMuHits < 3 || NMissedStations > 0)) continue;
-
- bool MuonFlag = 0;
- if (IndMuTrack[ists] >= 0) {
- CbmMuchTrack* Tmu = (CbmMuchTrack*) fMuchTracks->At(IndMuTrack[ists]);
- if (Tmu && Tmu->GetNHits() + Tmu->GetNMissedHits() >= 3)
- MuonFlag = Tmu->GetNMissedHits() <= 2 && Tmu->GetNMissedStations() == 0;
- }
-
- //cout<<ok<<" "<<Tmu->GetNHits()<<" "<< Tmu->GetNMissedHits()<<endl;
- if (Muon) { // muon
- if (Signal) {
- fhPerfSignal->Fill(P, MuonFlag);
- if (P < 5) {
- NMuS5++;
- NMuSRec5 += MuonFlag;
- }
- else {
- NMuS10++;
- NMuSRec10 += MuonFlag;
- }
- }
- fhPerfAll->Fill(P, MuonFlag);
- if (P < 5) {
- NMu5++;
- NMuRec5 += MuonFlag;
- }
- else {
- NMu10++;
- NMuRec10 += MuonFlag;
- }
- if (!MuonFlag) was_problem = 1;
- }
- else {
- fhGhost->Fill(P, MuonFlag);
- if (P < 5) {
- NGhost5++;
- NGhostRec5 += MuonFlag;
- }
- else {
- NGhost10++;
- NGhostRec10 += MuonFlag;
- }
- }
- }
-
- cout << " N Signal Mu Total/Reconstructed <5Gev = " << NMuS5 << "/" << NMuSRec5 << "; >=5Gev = " << NMuS10 << "/"
- << NMuSRec10 << endl;
- cout << " N Mu Total/Reconstructed <5Gev = " << NMu5 << "/" << NMuRec5 << "; >=5Gev = " << NMu10 << "/" << NMuRec10
- << endl;
- cout << " N Ghost Total/Reconstructed <5Gev = " << NGhost5 << "/" << NGhostRec5 << "; >=5Gev = " << NGhost10 << "/"
- << NGhostRec10 << endl;
-
- //if( was_problem ) cout<<"\n problem\n"<<endl;//return;
-
- vector<CbmL1MuchHit> vMuchHits;
-
- for (int ih = 0; ih < NHits; ih++) {
- CbmMuchHit* h = (CbmMuchHit*) fMuchHits->At(ih);
- CbmL1MuchHit m(h, ih);
- vMuchHits.push_back(m);
- }
-
- /*
- for( int mcTrackID=0; mcTrackID<NMCTracks; mcTrackID++ ){
- CbmMCTrack* mcTrack= (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if( !mcTrack ) continue;
- //cout<<"mc track "<<mcTrack->GetMotherID()<<"->"<<mcTrackID
- //<<" ( pdg "<<mcTrack->GetPdgCode()<<", p "<<mcTrack->GetMomentum().Mag()<<" ) : ";
- for( int ist=0; ist<MuNStations; ist++ ){
- //if( ist%3 ==0 ) cout<<" | ";
- CbmMuchPoint *mp = 0;
- for( int ih=0; ih<NHits; ih++ ){
- CbmL1MuchHit &h = vMuchHits[ih];
- if( h.iStation != ist ) continue;
- CbmMuchHit *mh = (CbmMuchHit*) fMuchHits->At(h.index);
- int ip = mh->GetRefIndex();
- if( ip<0 ) continue;
- CbmMuchPoint *pp = (CbmMuchPoint *) fMuchPoints->At(ip);
- if( mcTrackID == pp->GetTrackID() ){
- mp = pp;
- break;
- }
- }
- //if(mp ) cout<<"+ ";
- //else cout<<"- ";
- }
- //cout<<endl;
- }
- */
-
- for (int itr = 0; itr < NStsTracks; itr++) {
- CbmStsTrack* Tsts = (CbmStsTrack*) fStsTracks->At(itr);
- if (!Tsts) continue;
- if (Tsts->GetNStsHits() + Tsts->GetNMvdHits() < 4) continue;
-
- CbmStsTrackMatch* m = (CbmStsTrackMatch*) fSTSTrackMatch->At(itr);
- if (!m) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- if (mcTrackID < 0) continue;
- CbmMCTrack* mcTrack = (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if (!mcTrack) continue;
- if (abs(mcTrack->GetPdgCode()) != 13) continue;
- double ratio = 1. * m->GetNofTrueHits() / (m->GetNofTrueHits() + m->GetNofWrongHits() + m->GetNofFakeHits());
- if (ratio < .7) continue;
-
- fStsFitter.DoFit(Tsts, 13); // refit with muon mass
- CbmL1MuchTrack tr;
- tr.SetStsTrack(Tsts);
- tr.StsID = itr;
- tr.NHits = 0;
- tr.NMissed = 0;
- tr.ok = 1;
-
- for (int ist = 0; ist < MuNStations; ist++) {
- double Zdet = CbmKF::Instance()->vMuchDetectors[ist].ZReference;
- tr.Extrapolate(Zdet);
- if (fabs(tr.T[4]) > 100.) break;
- for (int ih = 0; ih < NHits; ih++) {
- CbmL1MuchHit& h = vMuchHits[ih];
- if (h.iStation != ist) continue;
- CbmMuchHit* mh = (CbmMuchHit*) fMuchHits->At(h.index);
- int ip = mh->GetRefIndex();
- if (ip < 0) continue;
- CbmMuchPoint* p = (CbmMuchPoint*) fMuchPoints->At(ip);
- if (p->GetTrackID() != mcTrackID) continue;
-
- double dx = tr.T[0] - h.FitPoint.x;
- double dy = tr.T[1] - h.FitPoint.y;
- double c0 = tr.C[0] + h.FitPoint.V[0];
- double c1 = tr.C[1] + h.FitPoint.V[1];
- double c2 = tr.C[2] + h.FitPoint.V[2];
- double chi2 = 0.5 * (dx * dx * c0 - 2 * dx * dy * c1 + dy * dy * c2) / (c0 * c2 - c1 * c1);
-
- histChi2[ist]->Fill(chi2);
- {
- TVector3 pos, mom;
- CbmL1MuchTrack tt = tr;
- p->Position(pos);
- p->Momentum(mom);
- tt.Propagate(pos.Z());
- Int_t j = ist + 20;
- double sx = sqrt(tt.C[0]);
- double sy = sqrt(tt.C[2]);
- double stx = sqrt(tt.C[5]);
- double sty = sqrt(tt.C[9]);
- double sqp = sqrt(tt.C[14]);
- double dx = (tt.T[0] - pos.X()) / sx;
- double dy = (tt.T[1] - pos.Y()) / sy;
- double dtx = (tt.T[2] - mom.X() / mom.Z()) / stx;
- double dty = (tt.T[3] - mom.Y() / mom.Z()) / sty;
- double dqp = (fabs(tt.T[4]) - 1. / mom.Mag()) / sqp;
- if (sx > 0.0) histPull_dx[j]->Fill(dx);
- if (sy > 0.0) histPull_dy[j]->Fill(dy);
- if (stx > 0.0) histPull_tx[j]->Fill(dtx);
- if (sty > 0.0) histPull_ty[j]->Fill(dty);
- if (sqp > 0.0) histPull_qp[j]->Fill(dqp);
- }
-
- double qp0 = tr.T[4];
- h.Filter(tr, 1, qp0);
- {
- TVector3 pos, mom;
- CbmL1MuchTrack tt = tr;
- p->Position(pos);
- p->Momentum(mom);
- tt.Propagate(pos.Z());
- if (tt.C[0] > 0.0) histPull_dx[ist]->Fill((tt.T[0] - pos.X()) / sqrt(tt.C[0]));
- if (tt.C[2] > 0.0) histPull_dy[ist]->Fill((tt.T[1] - pos.Y()) / sqrt(tt.C[2]));
- if (tt.C[5] > 0.0) histPull_tx[ist]->Fill((tt.T[2] - mom.X() / mom.Z()) / sqrt(tt.C[5]));
- if (tt.C[9] > 0.0) histPull_ty[ist]->Fill((tt.T[3] - mom.Y() / mom.Z()) / sqrt(tt.C[9]));
- if (tt.C[14] > 0.0) histPull_qp[ist]->Fill((fabs(tt.T[4]) - 1. / mom.Mag()) / sqrt(tt.C[14]));
- }
- }
- }
- } // tracks
-
-
- if (EventCounter < 100 && EventCounter % 10 == 0 || EventCounter >= 100 && EventCounter % 100 == 0) Write();
-}
-
-
-void CbmL1MuchFinderQa::Write()
-{
- TFile HistoFile("MuRecQa.root", "RECREATE");
- writedir2current(histodir);
- HistoFile.Close();
-}
-
-void CbmL1MuchFinderQa::writedir2current(TObject* obj)
-{
- if (!obj->IsFolder()) obj->Write();
- else {
- TDirectory* cur = gDirectory;
- TDirectory* sub = cur->mkdir(obj->GetName());
- sub->cd();
- TList* listSub = ((TDirectory*) obj)->GetList();
- TIter it(listSub);
- while (TObject* obj1 = it())
- writedir2current(obj1);
- cur->cd();
- }
-}
diff --git a/reco/L1/OffLineInterface/CbmL1MuchFinderQa.h b/reco/L1/OffLineInterface/CbmL1MuchFinderQa.h
deleted file mode 100644
index d736a23a37..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchFinderQa.h
+++ /dev/null
@@ -1,72 +0,0 @@
-/* Copyright (C) 2007-2009 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef CBM_L1_MuchFinderQa_h
-#define CBM_L1_MuchFinderQa_h
-
-#include "CbmStsKFTrackFitter.h"
-
-#include "FairTask.h"
-
-#include "TH1.h"
-#include "TLorentzVector.h"
-
-
-class TClonesArray;
-class TProfile;
-
-class CbmL1MuchFinderQa : public FairTask {
-public:
- /** Constructor **/
- CbmL1MuchFinderQa(const char* name = "CbmL1MuchFinderQa", Int_t iVerbose = 1);
-
- /** Destructor **/
- ~CbmL1MuchFinderQa();
-
- /// * FairTask methods
-
- /** Intialisation at begin of run. To be implemented in the derived class.
- *@value Success If not kSUCCESS, task will be set inactive.
- **/
- InitStatus Init();
-
- /** Reinitialisation.
- *@value Success If not kSUCCESS, task will be set inactive.
- **/
- InitStatus ReInit();
-
- /** Intialise parameter containers.
- **/
- void SetParContainers();
-
- void Exec(Option_t* option);
-
- /** Action after each event. **/
- void Finish();
-
-private:
- TClonesArray* fMuchPoints; //! Much MC points
- TClonesArray* fMuchHits; //! Much Hits
- TClonesArray* fStsTracks;
- TClonesArray* fMCTracks;
- TClonesArray* fMuchTracks; //! Much tracks
- TClonesArray* fSTSTrackMatch;
-
- CbmVertex* fPrimVtx;
- CbmStsKFTrackFitter fStsFitter;
-
- TDirectory* histodir;
-
- void Write();
- void writedir2current(TObject* obj);
-
- TProfile *fhPerfSignal, *fhPerfAll, *fhGhost;
- TH1F *histPull_dx[40], *histPull_dy[40], *histPull_tx[40], *histPull_ty[40], *histPull_qp[40], *histChi2[40];
-
-public:
- ClassDef(CbmL1MuchFinderQa, 1);
-};
-
-
-#endif
diff --git a/reco/L1/OffLineInterface/CbmL1MuchHit.cxx b/reco/L1/OffLineInterface/CbmL1MuchHit.cxx
deleted file mode 100644
index f6099225e6..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchHit.cxx
+++ /dev/null
@@ -1,36 +0,0 @@
-/* Copyright (C) 2007 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#include "CbmL1MuchHit.h"
-
-#include "CbmKF.h"
-#include "CbmKFTrackInterface.h"
-#include "CbmMuchHit.h"
-
-ClassImp(CbmL1MuchHit);
-
-void CbmL1MuchHit::Create(CbmMuchHit* h, int ind)
-{
- FitPoint.x = h->GetX();
- FitPoint.y = h->GetY();
- FitPoint.z = h->GetZ();
- FitPoint.V[0] = h->GetDx() * h->GetDx();
- FitPoint.V[1] = 0;
- FitPoint.V[2] = h->GetDy() * h->GetDy();
-
- CbmKF* KF = CbmKF::Instance();
- iStation = h->GetStationNr() - 1;
- MaterialIndex = KF->GetMaterialIndex(KF->MuchStation2MCIDMap[iStation]);
- time = h->GetTime(0);
- busy = 0;
- index = ind;
-}
-
-Int_t CbmL1MuchHit::Filter(CbmKFTrackInterface& track, Bool_t downstream, Double_t& QP0)
-{
- Bool_t err = 0;
- err = err || track.Propagate(FitPoint.z, QP0);
- err = err || FitPoint.Filter(track);
- return err;
-}
diff --git a/reco/L1/OffLineInterface/CbmL1MuchHit.h b/reco/L1/OffLineInterface/CbmL1MuchHit.h
deleted file mode 100644
index bc8e9b1204..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchHit.h
+++ /dev/null
@@ -1,34 +0,0 @@
-/* Copyright (C) 2007 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef CBM_L1_MuchHit_h
-#define CBM_L1_MuchHit_h
-
-#include "CbmKFHit.h"
-#include "CbmKFPixelMeasurement.h"
-
-class CbmMuchHit;
-class CbmKFTrackInterface;
-
-class CbmL1MuchHit : public CbmKFHit {
-public:
- CbmL1MuchHit() {}
- CbmL1MuchHit(CbmMuchHit* h, int index) { Create(h, index); }
- ~CbmL1MuchHit() {}
-
- void Create(CbmMuchHit* h, int index);
- Int_t Filter(CbmKFTrackInterface& track, Bool_t downstream, Double_t& QP0);
-
- CbmKFPixelMeasurement FitPoint;
-
- int index;
- int iStation;
- Double_t time;
- bool busy;
-
-public:
- ClassDef(CbmL1MuchHit, 1);
-};
-
-#endif
diff --git a/reco/L1/OffLineInterface/CbmL1MuchTrack.cxx b/reco/L1/OffLineInterface/CbmL1MuchTrack.cxx
deleted file mode 100644
index d659ad01b3..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchTrack.cxx
+++ /dev/null
@@ -1,17 +0,0 @@
-/* Copyright (C) 2007 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#include "CbmL1MuchTrack.h"
-
-#include "CbmKFMath.h"
-#include "CbmStsTrack.h"
-
-ClassImp(CbmL1MuchTrack);
-
-void CbmL1MuchTrack::SetStsTrack(CbmStsTrack* track)
-{
- CbmKFMath::CopyTrackParam2TC(track->GetParamLast(), T, C);
- chi2 = 0; //track->GetChi2();
- NDF = 0; //track->GetNDF();
-}
diff --git a/reco/L1/OffLineInterface/CbmL1MuchTrack.h b/reco/L1/OffLineInterface/CbmL1MuchTrack.h
deleted file mode 100644
index 12987e1723..0000000000
--- a/reco/L1/OffLineInterface/CbmL1MuchTrack.h
+++ /dev/null
@@ -1,49 +0,0 @@
-/* Copyright (C) 2007 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef CBM_L1_MuchTrack_h
-#define CBM_L1_MuchTrack_h
-
-#include "CbmKFTrackInterface.h"
-#include "CbmL1MuchHit.h"
-
-#include <vector>
-
-class CbmStsTrack;
-
-
-class CbmL1MuchTrack : public CbmKFTrackInterface {
-public:
- CbmL1MuchTrack() {}
- ~CbmL1MuchTrack() {}
-
- double* GetTrack() { return T; }
- double* GetCovMatrix() { return C; }
- double& GetRefChi2() { return chi2; }
- int& GetRefNDF() { return NDF; }
- double GetMass() { return 0.1057; }
- bool IsElectron() { return 0; }
- int GetNOfHits() { return vHits.size(); }
- CbmKFHit* GetHit(int i) { return vHits[i]; }
-
- void SetStsTrack(CbmStsTrack* track);
-
- double T[6], C[15], chi2;
- int NDF;
- std::vector<CbmL1MuchHit*> vHits;
- int NHits, NMissed, NMissedStations;
- bool ok;
- bool stopped;
- int StsID;
-
- static bool Compare(const CbmL1MuchTrack* p1, const CbmL1MuchTrack* p2)
- {
- return (p1->NHits > p2->NHits) || (p1->NHits == p2->NHits) && (p1->chi2 < p2->chi2);
- }
-
-public:
- ClassDef(CbmL1MuchTrack, 1);
-};
-
-#endif
diff --git a/reco/L1/OffLineInterface/CbmL1SttHit.cxx b/reco/L1/OffLineInterface/CbmL1SttHit.cxx
deleted file mode 100644
index cf04ec174b..0000000000
--- a/reco/L1/OffLineInterface/CbmL1SttHit.cxx
+++ /dev/null
@@ -1,56 +0,0 @@
-/* Copyright (C) 2008 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#include "CbmL1SttHit.h"
-
-#include "CbmKF.h"
-#include "CbmKFTrackInterface.h"
-#include "CbmSttHit.h"
-
-#include "TMath.h"
-
-ClassImp(CbmL1SttHit);
-
-void CbmL1SttHit::Create(CbmSttHit* h, int ind)
-{
- static const Double_t angle = 10. * TMath::DegToRad();
-
- Double_t phi = 0.; // rotation angle
- Int_t plane = h->GetPlaneID();
- Int_t irot = (plane - 1) % 6;
- irot = irot / 2;
- //if (irot == 1) phi = angle;
- //else if (irot == 2) phi = -angle;
- if (irot == 1) phi = -angle;
- else if (irot == 2)
- phi = angle;
- FitPoint.Set(h->GetZ(), h->GetU(), phi, h->GetDx() * h->GetDx());
- MaterialIndex = 0;
- time = 0.;
- iStation = h->GetPlaneID() - 1;
-
- /*
- FitPoint.x = h->GetX();
- FitPoint.y = h->GetY();
- FitPoint.z = h->GetZ();
- FitPoint.V[0] = h->GetDx()*h->GetDx();
- FitPoint.V[1] = 0;
- FitPoint.V[2] = h->GetDy()*h->GetDy();
-
- CbmKF *KF = CbmKF::Instance();
- iStation = h->GetStationNr()-1;
- MaterialIndex = KF->GetMaterialIndex(KF->MuchStation2MCIDMap[iStation]);
- time = h->GetTime(0);
- */
- busy = 0;
- index = ind;
-}
-
-Int_t CbmL1SttHit::Filter(CbmKFTrackInterface& track, Bool_t downstream, Double_t& QP0)
-{
- Bool_t err = 0;
- err = err || track.Propagate(FitPoint.z, QP0);
- err = err || FitPoint.Filter(track);
- return err;
-}
diff --git a/reco/L1/OffLineInterface/CbmL1SttHit.h b/reco/L1/OffLineInterface/CbmL1SttHit.h
deleted file mode 100644
index 2593d7efa2..0000000000
--- a/reco/L1/OffLineInterface/CbmL1SttHit.h
+++ /dev/null
@@ -1,36 +0,0 @@
-/* Copyright (C) 2008 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef CBM_L1_STTHIT_h
-#define CBM_L1_STTHIT_h
-
-#include "CbmKFHit.h"
-//#include "CbmKFPixelMeasurement.h"
-#include "CbmKFUMeasurement.h"
-
-class CbmSttHit;
-class CbmKFTrackInterface;
-
-class CbmL1SttHit : public CbmKFHit {
-public:
- CbmL1SttHit() {}
- CbmL1SttHit(CbmSttHit* h, int index) : CbmKFHit() { Create(h, index); }
- ~CbmL1SttHit() {}
-
- void Create(CbmSttHit* h, int index);
- Int_t Filter(CbmKFTrackInterface& track, Bool_t downstream, Double_t& QP0);
-
- //CbmKFPixelMeasurement FitPoint;
- CbmKFUMeasurement FitPoint;
-
- int index;
- int iStation;
- Double_t time;
- bool busy;
-
-public:
- ClassDef(CbmL1SttHit, 1);
-};
-
-#endif
diff --git a/reco/L1/OffLineInterface/CbmL1SttTrack.cxx b/reco/L1/OffLineInterface/CbmL1SttTrack.cxx
deleted file mode 100644
index f697745b89..0000000000
--- a/reco/L1/OffLineInterface/CbmL1SttTrack.cxx
+++ /dev/null
@@ -1,25 +0,0 @@
-/* Copyright (C) 2008 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#include "CbmL1SttTrack.h"
-
-#include "CbmKFMath.h"
-#include "CbmMuchTrack.h"
-#include "CbmStsTrack.h"
-
-ClassImp(CbmL1SttTrack);
-
-void CbmL1SttTrack::SetStsTrack(CbmStsTrack* track)
-{
- CbmKFMath::CopyTrackParam2TC(track->GetParamLast(), T, C);
- chi2 = 0; //track->GetChi2();
- NDF = 0; //track->GetNDF();
-}
-
-void CbmL1SttTrack::SetMuchTrack(CbmMuchTrack* track)
-{
- CbmKFMath::CopyTrackParam2TC(track->GetMuchTrack(), T, C);
- chi2 = 0; //track->GetChi2();
- NDF = 0; //track->GetNDF();
-}
diff --git a/reco/L1/OffLineInterface/CbmL1SttTrack.h b/reco/L1/OffLineInterface/CbmL1SttTrack.h
deleted file mode 100644
index 401e478e01..0000000000
--- a/reco/L1/OffLineInterface/CbmL1SttTrack.h
+++ /dev/null
@@ -1,50 +0,0 @@
-/* Copyright (C) 2008 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef CBM_L1_STTTRACK_h
-#define CBM_L1_STTTRACK_h
-
-#include "CbmKFTrackInterface.h"
-#include "CbmL1SttHit.h"
-
-#include <vector>
-
-class CbmStsTrack;
-class CbmMuchTrack;
-
-class CbmL1SttTrack : public CbmKFTrackInterface {
-public:
- CbmL1SttTrack() {}
- ~CbmL1SttTrack() {}
-
- double* GetTrack() { return T; }
- double* GetCovMatrix() { return C; }
- double& GetRefChi2() { return chi2; }
- int& GetRefNDF() { return NDF; }
- double GetMass() { return 0.1057; }
- bool IsElectron() { return 0; }
- int GetNOfHits() { return vHits.size(); }
- CbmKFHit* GetHit(int i) { return vHits[i]; }
-
- void SetStsTrack(CbmStsTrack* track);
- void SetMuchTrack(CbmMuchTrack* track);
-
- double T[6], C[15], chi2;
- int NDF;
- std::vector<CbmL1SttHit*> vHits;
- int NHits, NMissed, NMissedStations;
- bool ok;
- bool stopped;
- int StsID;
-
- static bool Compare(const CbmL1SttTrack* p1, const CbmL1SttTrack* p2)
- {
- return (p1->NHits > p2->NHits) || (p1->NHits == p2->NHits) && (p1->chi2 < p2->chi2);
- }
-
-public:
- ClassDef(CbmL1SttTrack, 1);
-};
-
-#endif
diff --git a/reco/L1/OffLineInterface/CbmL1SttTrackFinder.cxx b/reco/L1/OffLineInterface/CbmL1SttTrackFinder.cxx
deleted file mode 100644
index a4aeacdcfd..0000000000
--- a/reco/L1/OffLineInterface/CbmL1SttTrackFinder.cxx
+++ /dev/null
@@ -1,411 +0,0 @@
-/* Copyright (C) 2008-2021 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-// -------------------------------------------------------------------------
-// ----- CbmL1SttTrackFinder source file -----
-// ----- Created 8/03/08 by A. Zinchenko -----
-// -------------------------------------------------------------------------
-
-#include "CbmL1SttTrackFinder.h"
-
-#include "CbmKF.h"
-#include "CbmKFHit.h"
-#include "CbmKFMaterial.h"
-#include "CbmKFMath.h"
-#include "CbmKFPixelMeasurement.h"
-#include "CbmKFTrackInterface.h"
-#include "CbmL1SttHit.h"
-#include "CbmL1SttTrack.h"
-#include "CbmMCTrack.h"
-#include "CbmMuchHit.h"
-#include "CbmMuchPoint.h"
-#include "CbmMuchTrack.h"
-#include "CbmStsKFTrackFitter.h"
-#include "CbmStsTrack.h"
-#include "CbmStsTrackMatch.h"
-#include "CbmSttHit.h"
-#include "CbmSttPoint.h"
-#include "CbmSttTrack.h"
-#include "CbmVertex.h"
-
-#include "FairRootManager.h"
-
-#include "TClonesArray.h"
-#include "TFile.h"
-#include "TLorentzVector.h"
-#include "TVector3.h"
-
-#include <iostream>
-#include <vector>
-
-#include <cmath>
-
-using std::cout;
-using std::endl;
-using std::fabs;
-using std::vector;
-
-ClassImp(CbmL1SttTrackFinder);
-
-CbmL1SttTrackFinder::CbmL1SttTrackFinder(const char* name, Int_t iVerbose) : FairTask(name, iVerbose)
-{
- fTrackCollection = new TClonesArray("CbmSttTrack", 100);
- histodir = 0;
-}
-
-
-CbmL1SttTrackFinder::~CbmL1SttTrackFinder() {}
-
-InitStatus CbmL1SttTrackFinder::Init() { return ReInit(); }
-
-InitStatus CbmL1SttTrackFinder::ReInit()
-{
- fSttPoints = (TClonesArray*) FairRootManager::Instance()->GetObject("SttPoint");
- fSttHits = (TClonesArray*) FairRootManager::Instance()->GetObject("SttHit");
- fMuchTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MuchTrack");
- fStsTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("StsTrack");
- fMCTracks = (TClonesArray*) FairRootManager::Instance()->GetObject("MCTrack");
- fSTSTrackMatch = (TClonesArray*) FairRootManager::Instance()->GetObject("StsTrackMatch");
- // fPrimVtx = (CbmVertex *) FairRootManager::Instance() ->GetObject("PrimaryVertex");
- // Get pointer to PrimaryVertex object from IOManager if it exists
- // The old name for the object is "PrimaryVertex" the new one
- // "PrimaryVertex." Check first for the new name
- fPrimVtx = dynamic_cast<CbmVertex*>(FairRootManager::Instance()->GetObject("PrimaryVertex."));
- if (nullptr == fPrimVtx) {
- fPrimVtx = dynamic_cast<CbmVertex*>(FairRootManager::Instance()->GetObject("PrimaryVertex"));
- }
- if (nullptr == fPrimVtx) {
- Error("CbmL1SttTrackFinder::ReInit", "vertex not found!");
- return kERROR;
- }
- fStsFitter.Init();
-
- FairRootManager::Instance()->Register("SttTrack", "Stt", fTrackCollection, IsOutputBranchPersistent("SttTrack"));
-
- cout << " **************************************************" << endl;
- if (fMuchTracks) cout << " *** Using Much tracks for Stt tracking *** " << endl;
- else
- cout << " *** Using Sts tracks for Stt tracking *** " << endl;
- cout << " **************************************************" << endl;
-
- return kSUCCESS;
-}
-
-void CbmL1SttTrackFinder::SetParContainers() {}
-
-void CbmL1SttTrackFinder::Finish() { Write(); }
-
-void CbmL1SttTrackFinder::Exec(Option_t* /*option*/)
-{
- const int MaxBranches = 50;
-
- static bool first = 1;
- fTrackCollection->Clear();
- static int EventCounter = 0;
- EventCounter++;
- cout << " SttRec event " << EventCounter << endl;
-
- //int MuNStations = CbmKF::Instance()->MuchStation2MCIDMap.size();
- static const Int_t nStations = CbmKF::Instance()->SttStationIDMap.size();
- //const int nStations = 18; // to be taken elsewhere !!!
-
- if (first) {
- first = 0;
- TDirectory* curdir = gDirectory;
- histodir = gDirectory->mkdir("SttRec");
- histodir->cd();
- fhNBranches = new TH1F("NBranches", "N Branches", MaxBranches, 0, MaxBranches);
- curdir->cd();
- }
-
- int NHits = fSttHits->GetEntriesFast();
- vector<CbmL1SttHit> vSttHits;
-
- for (int ih = 0; ih < NHits; ++ih) {
- CbmSttHit* h = (CbmSttHit*) fSttHits->UncheckedAt(ih);
- CbmL1SttHit m(h, ih);
- vSttHits.push_back(m);
- }
-
- vector<CbmL1SttTrack> vTracks;
-
- Int_t nStsTracks;
- TClonesArray* seedTracks;
- if (fMuchTracks) seedTracks = fMuchTracks;
- else
- seedTracks = fStsTracks;
- nStsTracks = seedTracks->GetEntriesFast();
- cout << " Seed tracks: " << nStsTracks << endl;
-
- CbmL1SttTrack Branches[MaxBranches];
- Double_t scatAng[MaxBranches] = {0.};
-
- for (int itr = 0; itr < nStsTracks; ++itr) {
-
- Int_t nOK = 0;
- TObject* sts = seedTracks->UncheckedAt(itr);
- if (!fMuchTracks) {
- if (((CbmStsTrack*) sts)->GetNStsHits() + ((CbmStsTrack*) sts)->GetNMvdHits() < 4) continue;
- }
-
- // MC
- /*
- if( 0&&fSTSTrackMatch && fMCTracks ){
- CbmStsTrackMatch *m = (CbmStsTrackMatch*) fSTSTrackMatch->At(itr);
- if( !m ) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- if (mcTrackID<0) continue;
- CbmMCTrack* mcTrack= (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if( !mcTrack ) continue;
- if( abs(mcTrack->GetPdgCode())!=13 ) continue;
- }
- */
- // Check if track passes all the planes
- if (1 && fSTSTrackMatch && fMCTracks) {
- Int_t itr1 = itr;
- if (fMuchTracks) itr1 = ((CbmMuchTrack*) sts)->GetStsTrackID();
- CbmStsTrackMatch* m = (CbmStsTrackMatch*) fSTSTrackMatch->At(itr1);
- if (!m) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- //if (mcTrackID<0) continue;
- //CbmMCTrack* mcTrack= (CbmMCTrack*) fMCTracks->At(mcTrackID);
- //if( !mcTrack ) continue;
- //if( abs(mcTrack->GetPdgCode())!=13 ) continue;
- Int_t hitPlanes[20];
- TVector3 mom0(-1e+7), mom1;
- for (Int_t i = 0; i < nStations; ++i)
- hitPlanes[i] = -1;
- for (int ih = 0; ih < NHits; ++ih) {
- CbmL1SttHit& h = vSttHits[ih];
- CbmSttHit* hit = (CbmSttHit*) fSttHits->UncheckedAt(h.index);
- CbmSttPoint* p = (CbmSttPoint*) fSttPoints->UncheckedAt(hit->GetRefIndex());
- if (p->GetTrackID() != mcTrackID) continue;
- if (p->GetStationNo() == 1 && TMath::Sqrt(p->GetX() * p->GetX() + p->GetY() * p->GetY()) > 220) continue;
- if (hitPlanes[h.iStation] < 0) {
- hitPlanes[h.iStation] = 1;
- ++nOK;
- }
- if (mom0[0] < -1e+5) p->Momentum(mom0);
- else
- p->Momentum(mom1);
- if (itr < MaxBranches) scatAng[itr] = TMath::Max(scatAng[itr], mom1.Angle(mom0) * TMath::RadToDeg());
- }
- if (nOK < nStations) {
- //cout << " Track " << mcTrackID << " has " << nOK << " points " << endl;
- //continue;
- }
- }
-
- if (!fMuchTracks) fStsFitter.DoFit((CbmStsTrack*) sts, 13); // refit with muon mass
-
- int NBranches = 1;
-
- fMuchTracks == 0x0 ? Branches[0].SetStsTrack((CbmStsTrack*) sts) : Branches[0].SetMuchTrack((CbmMuchTrack*) sts);
- Branches[0].StsID = itr;
- Branches[0].NHits = 0;
- Branches[0].NMissed = 0;
- Branches[0].NMissedStations = 0;
- Branches[0].ok = 1;
- Branches[0].stopped = 0;
- Branches[0].vHits.clear();
- //cout<<"Sts track N "<<itr<<" with initial mom="<<1./Branches[0].T[4]<<endl;
- for (Int_t ist = 0; ist < nStations; ++ist) {
-
- int NBranchesOld = NBranches;
-
- for (Int_t ibr = 0; ibr < NBranchesOld; ++ibr) {
- CbmL1SttTrack& tr = Branches[ibr];
- if (tr.stopped) continue;
- //if( ist%3 ==0 ) cout<<" | ";
- //double Zdet = CbmKF::Instance()->vMuchDetectors[ist].ZReference;
- double Zdet = CbmKF::Instance()->vSttMaterial[ist].ZReference;
- //cout << Zdet << endl;
- //double Zdet = zPlanes[ist];
- tr.Extrapolate(Zdet);
- if (fabs(tr.T[4]) > 100.) tr.stopped = 1;
- if (1. < 0.5 * fabs(tr.T[4])) tr.stopped = 1; // 0.5 GeV, stop transport
- //if( tr.stopped ) cout<<"Sts track N "<<itr<<" stopped at Mu station "<<ist
- //<<" with mom="<<1./tr.T[4]<<endl;
- if (tr.stopped) continue;
- /*
- if( CbmKF::Instance()->vMuchDetectors[ist].IsOutside( tr.T[0], tr.T[1] ) ){
- //cout<<" out ";
- tr.NMissedStations++;
- continue;
- }
- */
-
- vector<int> NewHits;
- Int_t firstHit = -1;
- Double_t uTr = 0.;
- for (int ih = 0; ih < NHits; ++ih) {
- CbmL1SttHit& h = vSttHits[ih];
- if (h.iStation != ist) continue;
-
- if (firstHit < 0) {
- // Track coordinate transformation
- uTr = tr.T[0] * h.FitPoint.phi_c + tr.T[1] * h.FitPoint.phi_s;
- //uTr = tr.T[0] * h.FitPoint.phi_c - tr.T[1] * h.FitPoint.phi_s;
- firstHit = 1;
- }
- CbmSttHit* hit = (CbmSttHit*) fSttHits->UncheckedAt(h.index);
- //cout << tr.T[0] << " " << tr.T[1] << " " << hit->GetX() << " " << hit->GetY() << " " << uTr << " " << h.FitPoint.u << " " << h.FitPoint.phi_c << " " << h.FitPoint.phi_s << endl;
-
- /*
- if(0){ // !!! Cut for the time of flight (ns)
- double hl = sqrt(h.FitPoint.x*h.FitPoint.x+h.FitPoint.y*h.FitPoint.y+h.FitPoint.z*h.FitPoint.z);
- double hp = 1./fabs(tr.T[4]);
- double texp = hl*sqrt(1. - 0.1057*0.1057/(hp*hp))/29.9792458; //ns
- if( h.time - texp > 1.0 ) continue;
- }
- */
- /*
- double dx = tr.T[0] - h.FitPoint.x;
- double dy = tr.T[1] - h.FitPoint.y;
- double c0 = tr.C[0] + h.FitPoint.V[0];
- double c1 = tr.C[1] + h.FitPoint.V[1];
- double c2 = tr.C[2] + h.FitPoint.V[2];
- double chi2 = 0.5*(dx*dx*c0-2*dx*dy*c1+dy*dy*c2)/(c0*c2-c1*c1);
- */
- Double_t du = uTr - h.FitPoint.u;
- //Double_t c0 = tr.C[0] + h.FitPoint.sigma2;
- //Double_t chi2 = du * du / c0; // w/out correlations at the moment !!!
- Double_t w = h.FitPoint.sigma2 + h.FitPoint.phi_cc * tr.C[0] + h.FitPoint.phi_2sc * tr.C[1]
- + h.FitPoint.phi_ss * tr.C[2];
- Double_t chi21 = du * du / w;
- //cout << " chi2 " << ist << " " << du << " " << chi21 << " " << chi21 << endl;
- if (chi21 <= 20.) NewHits.push_back(ih);
- //if ( chi21 <= 100. ) NewHits.push_back( ih );
- }
- int nnew = NewHits.size();
- for (int ih = 1; ih < nnew; ++ih) {
- if (NBranches >= MaxBranches) break;
- CbmL1SttTrack& t = Branches[NBranches++];
- t = tr;
- CbmL1SttHit& h = vSttHits[NewHits[ih]];
- t.vHits.push_back(&h);
- t.NHits++;
- double qp0 = t.T[4];
- h.Filter(t, 1, qp0);
- }
- if (nnew > 0) {
- CbmL1SttHit& h = vSttHits[NewHits[0]];
- tr.vHits.push_back(&h);
- tr.NHits++;
- double qp0 = tr.T[4];
- h.Filter(tr, 1, qp0);
- }
- else
- tr.NMissed++;
- } // for( int ibr=0; ibr<NBranchesOld;
- } // for( int ist=0; ist<nStations;
- int bestbr = 0;
- for (int ibr = 1; ibr < NBranches; ++ibr) {
- if ((Branches[ibr].NHits > Branches[bestbr].NHits)
- || (Branches[ibr].NHits == Branches[bestbr].NHits) && (Branches[ibr].chi2 < Branches[bestbr].chi2))
- bestbr = ibr;
- }
- vTracks.push_back(Branches[bestbr]);
- // MC
- /*
- if( fSTSTrackMatch && fMCTracks ){
- CbmStsTrackMatch *m = (CbmStsTrackMatch*) fSTSTrackMatch->At(itr);
- if( !m ) continue;
- Int_t mcTrackID = m->GetMCTrackId();
- if (mcTrackID<0) continue;
- CbmMCTrack* mcTrack= (CbmMCTrack*) fMCTracks->At(mcTrackID);
- if( !mcTrack ) continue;
- if( abs(mcTrack->GetPdgCode())==13 ) fhNBranches->Fill(NBranches);
- }
- */
- if (nOK == nStations) fhNBranches->Fill(NBranches);
- else
- (vTracks.back()).ok = kFALSE;
- //cout << itr << " " << nOK << " " << (vTracks.back()).ok << endl;
- } // for( int itr=0; itr<nStsTracks;
-
- int NTracks = vTracks.size();
- cout << "NTracks=" << NTracks << endl;
- //sort tracks and check for common muon hits
-
- vector<CbmL1SttTrack*> vpTracks;
- for (int i = 0; i < NTracks; ++i)
- vpTracks.push_back(&(vTracks[i]));
- sort(vpTracks.begin(), vpTracks.end(), CbmL1SttTrack::Compare);
-
- int NOutTracks = fTrackCollection->GetEntriesFast();
-
- for (int it = 0; it < NTracks; ++it) {
- CbmL1SttTrack& tr = *vpTracks[it];
- if (tr.NDF <= 0 || tr.chi2 > 10. * tr.NDF) {
- //tr.ok = 0;
- //continue;
- }
- int nbusy = 0;
- for (int ih = 0; ih < tr.NHits; ++ih)
- if (tr.vHits[ih]->busy) nbusy++;
- if (0 && nbusy > 2) {
- tr.ok = 0;
- continue;
- }
- // Check if track contains hit from all doublets in all stations
- Int_t nDoubl[20] = {0};
- for (Int_t ih = 0; ih < tr.NHits; ++ih) {
- Int_t i2 = tr.vHits[ih]->iStation / 2;
- if (nDoubl[i2] == 0) nDoubl[i2]++;
- }
- Int_t nHit = 0;
- for (Int_t i = 0; i < nStations / 2; ++i)
- nHit += nDoubl[i];
- if (nHit < nStations / 2) continue;
- if (!(tr.ok)) continue;
-
- {
- new ((*fTrackCollection)[NOutTracks]) CbmSttTrack();
- CbmSttTrack* track = (CbmSttTrack*) fTrackCollection->At(NOutTracks++);
- track->SetChi2(tr.GetRefChi2());
- track->SetNDF(tr.GetRefNDF());
- FairTrackParam tp;
- CbmKFMath::CopyTC2TrackParam(&tp, tr.T, tr.C);
- track->SetSttTrack(&tp);
- track->SetStsTrackID(tr.StsID);
- int nh = 0;
- for (vector<CbmL1SttHit*>::iterator i = tr.vHits.begin(); i != tr.vHits.end(); i++) {
- if (++nh > 30) break;
- track->AddHitIndex((*i)->index);
- }
- track->SetNMissedHits(tr.NMissed);
- track->SetNMissedStations(tr.NMissedStations);
- if (tr.StsID < MaxBranches) track->SetScatAng(scatAng[tr.StsID]);
- }
- for (int ih = 0; ih < tr.NHits; ++ih)
- tr.vHits[ih]->busy = 1;
- }
-
- if (EventCounter < 100 && EventCounter % 10 == 0 || EventCounter >= 100 && EventCounter % 100 == 0) Write();
- cout << "end of SttRec " << fTrackCollection->GetEntriesFast() << endl;
-}
-
-
-void CbmL1SttTrackFinder::Write()
-{
- TFile HistoFile("SttRec.root", "RECREATE");
- writedir2current(histodir);
- HistoFile.Close();
-}
-
-void CbmL1SttTrackFinder::writedir2current(TObject* obj)
-{
- if (!obj->IsFolder()) obj->Write();
- else {
- TDirectory* cur = gDirectory;
- TDirectory* sub = cur->mkdir(obj->GetName());
- sub->cd();
- TList* listSub = ((TDirectory*) obj)->GetList();
- TIter it(listSub);
- while (TObject* obj1 = it())
- writedir2current(obj1);
- cur->cd();
- }
-}
diff --git a/reco/L1/OffLineInterface/CbmL1SttTrackFinder.h b/reco/L1/OffLineInterface/CbmL1SttTrackFinder.h
deleted file mode 100644
index 117342f301..0000000000
--- a/reco/L1/OffLineInterface/CbmL1SttTrackFinder.h
+++ /dev/null
@@ -1,70 +0,0 @@
-/* Copyright (C) 2008-2009 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
- SPDX-License-Identifier: GPL-3.0-only
- Authors: Sergey Gorbunov [committer] */
-
-#ifndef CBM_L1_STTTRACKFINDER_H
-#define CBM_L1_STTTRACKFINDER_H
-
-#include "CbmStsKFTrackFitter.h"
-
-#include "FairTask.h"
-//#include "TLorentzVector.h"
-#include "TH1.h"
-
-
-class TClonesArray;
-
-class CbmL1SttTrackFinder : public FairTask {
-public:
- /** Constructor **/
- CbmL1SttTrackFinder(const char* name = "CbmL1SttTrackFinder", Int_t iVerbose = 1);
-
- /** Destructor **/
- ~CbmL1SttTrackFinder();
-
- /// * FairTask methods
-
- /** Intialisation at begin of run. To be implemented in the derived class.
- *@value Success If not kSUCCESS, task will be set inactive.
- **/
- InitStatus Init();
-
- /** Reinitialisation.
- *@value Success If not kSUCCESS, task will be set inactive.
- **/
- InitStatus ReInit();
-
- /** Intialise parameter containers.
- **/
- void SetParContainers();
-
- void Exec(Option_t* option);
-
- /** Action after each event. **/
- void Finish();
-
-private:
- TClonesArray* fSttPoints; //! STT MC points
- TClonesArray* fSttHits; //! STT Hits
- TClonesArray* fMuchTracks;
- TClonesArray* fStsTracks;
- TClonesArray* fMCTracks;
- TClonesArray* fSTSTrackMatch;
- TClonesArray* fTrackCollection; //! STT tracks
-
- CbmVertex* fPrimVtx;
- CbmStsKFTrackFitter fStsFitter;
-
- TDirectory* histodir;
-
- void Write();
- void writedir2current(TObject* obj);
-
- TH1F* fhNBranches;
-
-public:
- ClassDef(CbmL1SttTrackFinder, 1);
-};
-
-
-#endif
--
GitLab