From 26e4d3f51111914b0df9957bce0eec6cea052d3b Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Fri, 12 Aug 2022 15:25:51 +0200
Subject: [PATCH] L1: added station index to hit; added phi distributions to
CbmL1Performance
---
reco/L1/CbmL1.cxx | 2 +-
reco/L1/CbmL1MCTrack.cxx | 2 +-
reco/L1/CbmL1Performance.cxx | 42 +++++++++++++++++++++++++++++-
reco/L1/L1Algo/L1CATrackFinder.cxx | 3 +--
reco/L1/L1Algo/L1ClonesMerger.cxx | 4 +--
reco/L1/L1Algo/L1TrackExtender.cxx | 16 ++++++------
reco/L1/L1Algo/L1TrackFitter.cxx | 23 ++++++++--------
7 files changed, 66 insertions(+), 26 deletions(-)
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 9be1b80d69..3b4896bee5 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -986,7 +986,7 @@ void CbmL1::Reconstruct(CbmEvent* event)
#ifdef USE_EVENT_NUMBER
h.n = mcp.event;
#endif
- const int ista = (*fpAlgo->fStripFlag)[h.f] / 4;
+ const int ista = h.iSt;
const L1Station& sta = fpAlgo->GetParameters()->GetStation(ista);
if (std::find(strips.begin(), strips.end(), h.f) != strips.end()) { // separate strips
diff --git a/reco/L1/CbmL1MCTrack.cxx b/reco/L1/CbmL1MCTrack.cxx
index e45a6bf590..cb065fef29 100644
--- a/reco/L1/CbmL1MCTrack.cxx
+++ b/reco/L1/CbmL1MCTrack.cxx
@@ -125,7 +125,7 @@ void CbmL1MCTrack::CalculateHitCont()
for (int ih = 0; ih < nhits; ih++) {
int jh = Hits[ih];
const L1Hit& h = (*algo->vHits)[jh];
- int ista = (*algo->fStripFlag)[h.f] / 4;
+ int ista = h.iSt;
if (ista - istaold == 1) ncont++;
else if (ista - istaold > 1) {
diff --git a/reco/L1/CbmL1Performance.cxx b/reco/L1/CbmL1Performance.cxx
index 820efc06a5..d5721ae82c 100644
--- a/reco/L1/CbmL1Performance.cxx
+++ b/reco/L1/CbmL1Performance.cxx
@@ -43,6 +43,7 @@
#include "TH2.h"
#include "TProfile.h"
#include <TFile.h>
+#include "TMath.h"
#include <iostream>
#include <list>
@@ -530,6 +531,18 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
static TH2F *h2_rest_nhits_vs_mom_prim, *h2_rest_nhits_vs_mom_sec, *h2_rest_fstation_vs_mom_prim,
*h2_rest_fstation_vs_mom_sec, *h2_rest_lstation_vs_fstation_prim, *h2_rest_lstation_vs_fstation_sec;
+ static TH1F* h_reg_phi_prim;
+ static TH1F* h_reg_phi_sec;
+ static TH1F* h_acc_phi_prim;
+ static TH1F* h_acc_phi_sec;
+ static TH1F* h_reco_phi_prim;
+ static TH1F* h_reco_phi_sec;
+ static TH1F* h_rest_phi_prim;
+ static TH1F* h_rest_phi_sec;
+ static TH1F* h_ghost_phi;
+ static TH1F* h_reco_phi;
+ static TH1F* h_notfound_phi;
+
static bool first_call = 1;
if (first_call) {
@@ -579,6 +592,15 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_rest_mom_prim = new TH1F("h_rest_mom_prim", "Momentum of not found primary tracks", 500, 0.0, 5.0);
h_rest_mom_sec = new TH1F("h_rest_mom_sec", "Momentum of not found secondary tracks", 500, 0.0, 5.0);
+ h_reg_phi_prim = new TH1F("h_reg_phi_prim", "Azimuthal angle of registered primary tracks", 1000, -3.15, 3.15);
+ h_reg_phi_sec = new TH1F("h_reg_phi_sec", "Azimuthal angle of registered secondary tracks", 1000, -3.15, 3.15);
+ h_acc_phi_prim = new TH1F("h_acc_phi_prim", "Azimuthal angle of accepted primary tracks", 1000, -3.15, 3.15);
+ h_acc_phi_sec = new TH1F("h_acc_phi_sec", "Azimuthal angle of accepted secondary tracks", 1000, -3.15, 3.15);
+ h_reco_phi_prim = new TH1F("h_reco_phi_prim", "Azimuthal angle of reconstructed primary tracks", 1000, -3.15, 3.15);
+ h_reco_phi_sec = new TH1F("h_reco_phi_sec", "Azimuthal angle of reconstructed secondary tracks", 1000, -3.15, 3.15);
+ h_rest_phi_prim = new TH1F("h_rest_phi_prim", "Azimuthal angle of not found primary tracks", 1000, -3.15, 3.15);
+ h_rest_phi_sec = new TH1F("h_rest_phi_sec", "Azimuthal angle of not found secondary tracks", 1000, -3.15, 3.15);
+
h_reg_nhits_prim = new TH1F("h_reg_nhits_prim", "Number of hits in registered primary track", 51, -0.1, 10.1);
h_reg_nhits_sec = new TH1F("h_reg_nhits_sec", "Number of hits in registered secondary track", 51, -0.1, 10.1);
h_acc_nhits_prim = new TH1F("h_acc_nhits_prim", "Number of hits in accepted primary track", 51, -0.1, 10.1);
@@ -589,6 +611,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_rest_nhits_sec = new TH1F("h_rest_nhits_sec", "Number of hits in not found secondary track", 51, -0.1, 10.1);
h_ghost_mom = new TH1F("h_ghost_mom", "Momentum of ghost tracks", 500, 0.0, 5.0);
+ h_ghost_phi = new TH1F("h_ghost_phi", "Azimuthal angle of ghost tracks", 1000, -3.15, 3.15);
h_ghost_nhits = new TH1F("h_ghost_nhits", "Number of hits in ghost track", 51, -0.1, 10.1);
h_ghost_fstation = new TH1F("h_ghost_fstation", "First station of ghost track", 50, -0.5, 10.0);
h_ghost_chi2 = new TH1F("h_ghost_chi2", "Chi2/NDF of ghost track", 50, -0.5, 10.0);
@@ -599,6 +622,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_ghost_purity = new TH1F("h_ghost_purity", "Ghost: percentage of correct hits", 100, -0.5, 100.5);
h_reco_mom = new TH1F("h_reco_mom", "Momentum of reco track", 50, 0.0, 5.0);
+ h_reco_phi = new TH1F("h_reco_phi", "Azimuthal angle of reco track", 1000, -3.15, 3.15);
h_reco_nhits = new TH1F("h_reco_nhits", "Number of hits in reco track", 50, 0.0, 10.0);
h_reco_station = new TH1F("h_reco_station", "First station of reco track", 50, -0.5, 10.0);
h_reco_chi2 = new TH1F("h_reco_chi2", "Chi2/NDF of reco track", 50, -0.5, 10.0);
@@ -631,6 +655,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_sec_r = new TH1F("h_sec_r", "R of sec MC track at the first hit", 50, 0.0, 15.0);
h_notfound_mom = new TH1F("h_notfound_mom", "Momentum of not found track", 50, 0.0, 5.0);
+ h_notfound_phi = new TH1F("h_notfound_phi", "Momentum of not found track", 50, 0.0, 5.0);
h_notfound_nhits = new TH1F("h_notfound_nhits", "Num hits of not found track", 50, 0.0, 10.0);
h_notfound_station = new TH1F("h_notfound_station", "First station of not found track", 50, 0.0, 10.0);
h_notfound_r = new TH1F("h_notfound_r", "R at first hit of not found track", 50, 0.0, 15.0);
@@ -772,7 +797,11 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
CbmL1Track* prtra = &(*rtraIt);
if ((prtra->Hits).size() < 1) continue;
{ // fill histos
- if (fabs(prtra->T[4]) > 1.e-10) h_reco_mom->Fill(fabs(1.0 / prtra->T[4]));
+ if (fabs(prtra->T[4]) > 1.e-10) h_reco_mom->Fill(fabs(1.0 / prtra->T[4])); // TODO: Is it a right precision? In FairTrackParam it is 1.e-4 (S.Zharko)
+ // NOTE: p = (TMath::Abs(fQp) > 1.e-4) ? 1. / TMath::Abs(fQp) : 1.e4; // FairTrackParam::Momentum(TVector3)
+ // h_reco_mom->Fill(TMath::Abs(prtra->T[4] > 1.e-4) ? 1. / TMath::Abs(prtra->T[4]) : 1.e+4); // this should be correct
+
+ h_reco_phi->Fill(TMath::ATan2(-prtra->T[3], -prtra->T[2])); // TODO: What is precision?
h_reco_nhits->Fill((prtra->Hits).size());
CbmL1HitStore& mh = fvHitStore[prtra->Hits[0]];
h_reco_station->Fill(mh.iStation);
@@ -803,6 +832,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_ghost_purity->Fill(100 * prtra->GetMaxPurity());
if (fabs(prtra->T[4]) > 1.e-10) {
h_ghost_mom->Fill(fabs(1.0 / prtra->T[4]));
+ h_ghost_phi->Fill(atan(prtra->T[3] / prtra->T[2])); // phi = atan(py / px) = atan(ty / tx)
h_ghost_Rmom->Fill(fabs(1.0 / prtra->T[4]));
}
h_ghost_nhits->Fill((prtra->Hits).size());
@@ -838,6 +868,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
double momentum = mtra.p;
double pt = sqrt(mtra.px * mtra.px + mtra.py * mtra.py);
double theta = acos(mtra.pz / mtra.p) * 180 / 3.1415;
+ double phi = TMath::ATan2(-mtra.py, -mtra.px);
h_mcp->Fill(momentum);
h_nmchits->Fill(nmchits);
@@ -850,6 +881,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_reg_MCmom->Fill(momentum);
if (mtra.IsPrimary()) {
h_reg_mom_prim->Fill(momentum);
+ h_reg_phi_prim->Fill(phi);
h_reg_prim_MCmom->Fill(momentum);
h_reg_nhits_prim->Fill(nSta);
h2_reg_nhits_vs_mom_prim->Fill(momentum, nSta);
@@ -858,6 +890,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
else {
h_reg_mom_sec->Fill(momentum);
+ h_reg_phi_sec->Fill(phi);
h_reg_sec_MCmom->Fill(momentum);
h_reg_nhits_sec->Fill(nSta);
if (momentum > 0.01) {
@@ -875,6 +908,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
h_acc_MCmom->Fill(momentum);
if (mtra.IsPrimary()) {
h_acc_mom_prim->Fill(momentum);
+ h_acc_phi_prim->Fill(phi);
h_acc_prim_MCmom->Fill(momentum);
h_acc_nhits_prim->Fill(nSta);
h2_acc_nhits_vs_mom_prim->Fill(momentum, nSta);
@@ -883,6 +917,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
else {
h_acc_mom_sec->Fill(momentum);
+ h_acc_phi_sec->Fill(phi);
h_acc_sec_MCmom->Fill(momentum);
h_acc_nhits_sec->Fill(nSta);
if (momentum > 0.01) {
@@ -936,6 +971,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
if (mtra.mother_ID < 0) { // primary
h_reco_mom_prim->Fill(momentum);
+ h_reco_phi_prim->Fill(phi);
h_reco_prim_MCmom->Fill(momentum);
h_reco_nhits_prim->Fill(nSta);
h2_reco_nhits_vs_mom_prim->Fill(momentum, nSta);
@@ -944,6 +980,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
else {
h_reco_mom_sec->Fill(momentum);
+ h_reco_phi_sec->Fill(phi);
h_reco_sec_MCmom->Fill(momentum);
h_reco_nhits_sec->Fill(nSta);
if (momentum > 0.01) {
@@ -955,6 +992,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
else {
h_notfound_mom->Fill(momentum);
+ h_notfound_phi->Fill(phi);
p_eff_all_vs_mom->Fill(momentum, 0.0);
p_eff_all_vs_nhits->Fill(nMCHits, 0.0);
p_eff_all_vs_pt->Fill(pt, 0.0);
@@ -970,6 +1008,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
if (mtra.mother_ID < 0) { // primary
h_rest_mom_prim->Fill(momentum);
+ h_rest_phi_prim->Fill(phi);
h_rest_nhits_prim->Fill(nSta);
h2_rest_nhits_vs_mom_prim->Fill(momentum, nSta);
h2_rest_fstation_vs_mom_prim->Fill(momentum, fh.iStation + 1);
@@ -977,6 +1016,7 @@ void CbmL1::HistoPerformance() // TODO: check if works correctly. Change vHitRe
}
else {
h_rest_mom_sec->Fill(momentum);
+ h_rest_phi_sec->Fill(phi);
h_rest_nhits_sec->Fill(nSta);
if (momentum > 0.01) {
h2_rest_nhits_vs_mom_sec->Fill(momentum, nSta);
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index 34ee032900..0747f42398 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -2258,8 +2258,7 @@ void L1Algo::CATrackFinder()
continue; //ghost suppression//find track with 3 hits only if it was created from a chain of triplets, but not from only one triplet
if (kGlobal != fTrackingMode && kMcbm != fTrackingMode) {
- if ((firstTripletLevel == 0)
- && (GetFStation((*fStripFlag)[(*vHitsUnused)[first_trip.GetLHit()].f]) != 0))
+ if ((firstTripletLevel == 0) && ((*vHitsUnused)[first_trip.GetLHit()].iSt != 0))
continue; // ghost supression // collect only MAPS tracks-triplets CHECK!!!
}
}
diff --git a/reco/L1/L1Algo/L1ClonesMerger.cxx b/reco/L1/L1Algo/L1ClonesMerger.cxx
index 26d1b3f94f..8329386dc7 100644
--- a/reco/L1/L1Algo/L1ClonesMerger.cxx
+++ b/reco/L1/L1Algo/L1ClonesMerger.cxx
@@ -64,10 +64,10 @@ void L1ClonesMerger::Exec(L1Vector<L1Track>& extTracks, L1Vector<L1HitIndex_t>&
#endif
for (int iTr = 0; iTr < nTracks; iTr++) {
firstHit[iTr] = start_hit;
- firstStation[iTr] = (*frAlgo.fStripFlag)[(*frAlgo.vHits)[extRecoHits[start_hit]].f] / 4;
+ firstStation[iTr] = (*frAlgo.vHits)[extRecoHits[start_hit]].iSt;
start_hit += extTracks[iTr].NHits - 1;
lastHit[iTr] = start_hit;
- lastStation[iTr] = (*frAlgo.fStripFlag)[(*frAlgo.vHits)[extRecoHits[start_hit]].f] / 4;
+ lastStation[iTr] = (*frAlgo.vHits)[extRecoHits[start_hit]].iSt;
start_hit++;
isStored[iTr] = false;
diff --git a/reco/L1/L1Algo/L1TrackExtender.cxx b/reco/L1/L1Algo/L1TrackExtender.cxx
index 9d38fe5ea5..77be8f855f 100644
--- a/reco/L1/L1Algo/L1TrackExtender.cxx
+++ b/reco/L1/L1Algo/L1TrackExtender.cxx
@@ -44,9 +44,9 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
const L1Hit& hit1 = (*vHits)[hits[iFirstHit + step]];
const L1Hit& hit2 = (*vHits)[hits[iFirstHit + 2 * step]];
- int ista0 = GetFStation((*fStripFlag)[hit0.f]);
- int ista1 = GetFStation((*fStripFlag)[hit1.f]);
- int ista2 = GetFStation((*fStripFlag)[hit2.f]);
+ int ista0 = hit0.iSt;
+ int ista1 = hit1.iSt;
+ int ista2 = hit2.iSt;
const L1Station& sta0 = fParameters.GetStation(ista0);
const L1Station& sta1 = fParameters.GetStation(ista1);
@@ -121,7 +121,7 @@ void L1Algo::BranchFitterFast(const L1Branch& t, L1TrackPar& T, const bool dir,
for (int i = iFirstHit + step; step * i <= step * iLastHit; i += step) {
const L1Hit& hit = (*vHits)[hits[i]];
ista_prev = ista;
- ista = GetFStation((*fStripFlag)[hit.f]);
+ ista = hit.iSt;
const L1Station& sta = fParameters.GetStation(ista);
@@ -199,15 +199,15 @@ void L1Algo::FindMoreHits(L1Branch& t, L1TrackPar& T, const bool dir,
const signed short int step = -2 * static_cast<int>(dir) + 1; // increment for station index
const int iFirstHit = (dir) ? 2 : t.NHits - 3;
- // int ista = GetFStation((*fStripFlag)[(*vHits)[t.Hits[iFirstHit]].f]) + 2*step; // current station. set to the end of track
+ // int ista = (*vHits)[t.Hits[iFirstHit]].iSt + 2*step; // current station. set to the end of track
const L1Hit& hit0 = (*vHits)[t.fHits[iFirstHit]]; // optimize
const L1Hit& hit1 = (*vHits)[t.fHits[iFirstHit + step]];
const L1Hit& hit2 = (*vHits)[t.fHits[iFirstHit + 2 * step]];
- const int ista0 = GetFStation((*fStripFlag)[hit0.f]);
- const int ista1 = GetFStation((*fStripFlag)[hit1.f]);
- const int ista2 = GetFStation((*fStripFlag)[hit2.f]);
+ const int ista0 = hit0.iSt;
+ const int ista1 = hit1.iSt;
+ const int ista2 = hit2.iSt;
const L1Station& sta0 = fParameters.GetStation(ista0);
const L1Station& sta1 = fParameters.GetStation(ista1);
diff --git a/reco/L1/L1Algo/L1TrackFitter.cxx b/reco/L1/L1Algo/L1TrackFitter.cxx
index f65c48511d..cc3f85ee64 100644
--- a/reco/L1/L1Algo/L1TrackFitter.cxx
+++ b/reco/L1/L1Algo/L1TrackFitter.cxx
@@ -56,9 +56,10 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
const L1Hit& hit1 = (*vHits)[hits[nHits - 2]];
const L1Hit& hit2 = (*vHits)[hits[nHits - 3]];
- int ista0 = (*fStripFlag)[hit0.f] / 4;
- int ista1 = (*fStripFlag)[hit1.f] / 4;
- int ista2 = (*fStripFlag)[hit2.f] / 4;
+ int ista0 = hit0.iSt;
+ int ista1 = hit1.iSt;
+ int ista2 = hit2.iSt;
+
//cout<<"back: ista012="<<ista0<<" "<<ista1<<" "<<ista2<<endl;
const L1Station& sta0 = fParameters.GetStation(ista0);
@@ -124,7 +125,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
for (int i = nHits - 2; i >= 0; i--) {
// if( fabs(T.qp[0])>2. ) break; // iklm. Don't know it need for
const L1Hit& hit = (*vHits)[hits[i]];
- ista = (*fStripFlag)[hit.f] / 4;
+ ista = hit.iSt;
const L1Station& sta = fParameters.GetStation(ista);
@@ -192,10 +193,10 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
const L1Hit& hit0 = (*vHits)[hits[0]];
const L1Hit& hit1 = (*vHits)[hits[1]];
const L1Hit& hit2 = (*vHits)[hits[2]];
-
- int ista0 = GetFStation((*fStripFlag)[hit0.f]);
- int ista1 = GetFStation((*fStripFlag)[hit1.f]);
- int ista2 = GetFStation((*fStripFlag)[hit2.f]);
+
+ int ista0 = hit0.iSt;
+ int ista1 = hit1.iSt;
+ int ista2 = hit2.iSt;
const L1Station& sta0 = fParameters.GetStation(ista0);
const L1Station& sta1 = fParameters.GetStation(ista1);
@@ -256,7 +257,7 @@ void L1Algo::KFTrackFitter_simple() // TODO: Add pipe.
for (int i = 1; i < nHits; i++) {
const L1Hit& hit = (*vHits)[hits[i]];
- ista = (*fStripFlag)[hit.f] / 4;
+ ista = hit.iSt;
const L1Station& sta = fParameters.GetStation(ista);
fvec u = hit.u;
fvec v = hit.v;
@@ -423,7 +424,7 @@ void L1Algo::L1KFTrackFitter()
int iSta[L1Constants::size::kMaxNstations];
for (i = 0; i < nHitsTrack; i++) {
const L1Hit& hit = (*vHits)[fRecoHits[start_hit++]];
- const int ista = (*fStripFlag)[hit.f] / 4;
+ const int ista = hit.iSt;
iSta[i] = ista;
w[ista][iVec] = 1.;
if (ista > fNstationsBeforePipe) w_time[ista][iVec] = 1.;
@@ -893,7 +894,7 @@ void L1Algo::L1KFTrackFitterMuch()
int nHitsTrackField = 0;
for (i = 0; i < nHitsTrack; i++) {
const L1Hit& hit = (*vHits)[fRecoHits[start_hit++]];
- const int ista = (*fStripFlag)[hit.f] / 4;
+ const int ista = hit.iSt;
if (ista < fNfieldStations) nHitsTrackField++;
iSta[i] = ista;
w[ista][iVec] = 1.;
--
GitLab