From 039cd35278824836543438eb582bbd826497e3be Mon Sep 17 00:00:00 2001
From: sgorbuno <se.gorbunov@gsi.de>
Date: Wed, 14 Jul 2021 22:27:36 +0000
Subject: [PATCH] add L1Algo/utils/* debug utilities to the compilation
---
reco/L1/CMakeLists.txt | 2 +
reco/L1/L1Algo/L1Algo.h | 2 +
.../utils/L1AlgoEfficiencyPerformance.cxx | 8 +
.../utils/L1AlgoEfficiencyPerformance.h | 50 +--
reco/L1/L1Algo/utils/L1AlgoPulls.cxx | 260 +++++++++++++++
reco/L1/L1Algo/utils/L1AlgoPulls.h | 297 ++----------------
6 files changed, 320 insertions(+), 299 deletions(-)
create mode 100644 reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
create mode 100644 reco/L1/L1Algo/utils/L1AlgoPulls.cxx
diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 5e5ababa52..488e129c6c 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -134,6 +134,8 @@ L1Algo/L1Fit.cxx
CbmL1MCTrack.cxx
L1Algo/utils/L1AlgoDraw.cxx
+L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
+L1Algo/utils/L1AlgoPulls.cxx
ParticleFinder/CbmL1PFFitter.cxx
ParticleFinder/CbmL1PFMCParticle.cxx
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index 16833a2976..b6798a39fd 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -149,6 +149,8 @@ public:
void L1KFTrackFitterMuch();
+ float GetMaxInvMom() const { return MaxInvMom[0]; }
+
/// ----- Input data -----
// filled in CbmL1::ReadEvent();
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.cxx b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
new file mode 100644
index 0000000000..e090a46b72
--- /dev/null
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.cxx
@@ -0,0 +1,8 @@
+/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer] */
+
+#include "L1AlgoEfficiencyPerformance.h"
+
+template class L1AlgoEfficiencyPerformance<2>;
+template class L1AlgoEfficiencyPerformance<3>;
diff --git a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
index 683da7e350..dc62da303a 100644
--- a/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
+++ b/reco/L1/L1Algo/utils/L1AlgoEfficiencyPerformance.h
@@ -14,23 +14,22 @@
#include "CbmL1.h"
#include "CbmL1Counters.h"
-#include "CbmL1Def.h"
#include <iostream>
#include <vector>
#include "L1Algo.h"
+#include "L1Def.h"
using std::cout;
using std::endl;
using std::vector;
-const int MaxNStations = 10;
template<int NHits>
class L1Tracklet {
public:
- L1Tracklet() : mcTrackId(-1), iStation(-1) {};
+ L1Tracklet() = default;
static bool compare(const L1Tracklet<NHits>& a, const L1Tracklet<NHits>& b)
{
@@ -41,15 +40,15 @@ public:
bool operator!=(const L1Tracklet<NHits>& a) { return !operator==(a); }
bool operator<(const L1Tracklet<NHits>& a) { return compare(*this, a); }
- int mcTrackId;
- int iStation;
+ int mcTrackId {-1};
+ int iStation {-1};
};
// reconstructed Tracklet
template<int NHits>
class L1RTracklet : public L1Tracklet<NHits> {
public:
- L1RTracklet() {}
+ L1RTracklet() = default;
L1RTracklet(THitI* ih_) : L1Tracklet<NHits>()
{
for (int iih = 0; iih < NHits; iih++) {
@@ -72,14 +71,14 @@ public:
}
bool operator!=(const L1RTracklet<NHits>& a) { return !operator==(a); }
- THitI i[NHits]; // indices of left, middle and right hits
+ THitI i[NHits] {0}; // indices of left, middle and right hits
};
// MC Tracklet
template<int NHits>
class L1MTracklet : public L1Tracklet<NHits> {
public:
- L1MTracklet() : isReconstructable(false) {}
+ L1MTracklet() = default;
void AddReconstructed(int recoId) { recoIds.push_back(recoId); };
void SetAsReconstructable() { isReconstructable = true; };
@@ -89,13 +88,13 @@ public:
bool IsReconstructable() { return isReconstructable; }; // is it possible to reconstruct
bool IsPrimary() { return mcMotherTrackId < 0; }; // is it primary or secondary
- int mcMotherTrackId;
- float p; // momentum
- vector<int> recoIds; // indices of reco tracklets in sorted recoArray
+ int mcMotherTrackId {-1};
+ float p {0.f}; // momentum
+ vector<int> recoIds {}; // indices of reco tracklets in sorted recoArray
- vector<int> hitIds[NHits]; // indices of hits in L1::vStsHits. Separated by stations
+ vector<int> hitIds[NHits] {}; // indices of hits in L1::vStsHits. Separated by stations
private:
- bool isReconstructable; // is it reconstructed
+ bool isReconstructable {false}; // is it reconstructed
};
// tracklets categories
@@ -118,7 +117,7 @@ class L1AlgoEfficiencyPerformance {
public:
typedef L1RTracklet<NHits> L1RecoTracklet;
typedef L1MTracklet<NHits> L1MCTracklet;
- L1AlgoEfficiencyPerformance() {};
+ L1AlgoEfficiencyPerformance() = default;
void Init();
bool AddOne(THitI* ih_); // add one recoTracklets. Return is it reconstructed or not
@@ -127,20 +126,25 @@ public:
void Print(TString title = "Triplets performance statistic",
bool station = 0); // Print efficiencies
+ static constexpr int MaxNStations = 10;
+
private:
+ L1AlgoEfficiencyPerformance<NHits>(const L1AlgoEfficiencyPerformance<NHits>&);
+ L1AlgoEfficiencyPerformance<NHits>& operator=(const L1AlgoEfficiencyPerformance<NHits>&);
+
void FillMC(); // collect mcTracklets from mcTracks
void MatchTracks();
- CbmL1* fL1;
+ CbmL1* fL1 {nullptr};
- vector<L1RecoTracklet> recoTracklets;
- vector<L1MCTracklet> mcTracklets;
+ vector<L1RecoTracklet> recoTracklets {};
+ vector<L1MCTracklet> mcTracklets {};
- TL1AlgoEfficiencies ntra; // current event efficiencies
- TL1AlgoEfficiencies NTRA; // efficiencies statistic
+ TL1AlgoEfficiencies ntra {}; // current event efficiencies
+ TL1AlgoEfficiencies NTRA {}; // efficiencies statistic
- TL1AlgoEfficiencies ntra_sta[MaxNStations]; // for each station separately
- TL1AlgoEfficiencies NTRA_STA[MaxNStations];
+ TL1AlgoEfficiencies ntra_sta[MaxNStations] {}; // for each station separately
+ TL1AlgoEfficiencies NTRA_STA[MaxNStations] {};
};
// ===================================================================================
@@ -202,7 +206,7 @@ void L1AlgoEfficiencyPerformance<NHits>::FillMC()
trlet.mcTrackId = mtra.ID;
trlet.mcMotherTrackId = mtra.mother_ID;
trlet.p = mtra.p;
- if (mtra.p > fL1->MinRecoMom) trlet.SetAsReconstructable();
+ if (mtra.p > 1. / fL1->algo->GetMaxInvMom()) trlet.SetAsReconstructable();
mcTracklets.push_back(trlet);
} // for Iter = stations
@@ -324,7 +328,7 @@ void L1AlgoEfficiencyPerformance<NHits>::CalculateEff()
// count tracklets
ntra_sta[iSta].Inc(reco, "total");
- if (mtra.p > fL1->MinRefMom) { // reference tracks
+ if (mtra.p > 1.) { // reference tracks
ntra_sta[iSta].Inc(reco, "fast");
if (mtra.IsPrimary()) { // reference primary
ntra_sta[iSta].Inc(reco, "fast_prim");
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.cxx b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
new file mode 100644
index 0000000000..258c820b91
--- /dev/null
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.cxx
@@ -0,0 +1,260 @@
+/* Copyright (C) 2010-2021 Frankfurt Institute for Advanced Studies, Goethe-Universität Frankfurt, Frankfurt
+ SPDX-License-Identifier: GPL-3.0-only
+ Authors: Sergey Gorbunov [committer] */
+
+#include "L1AlgoPulls.h"
+
+const TString L1TrackParametersNames[TL1TrackParameters::NParameters] = {"x", "y", "tx", "ty", "qp"};
+
+void L1AlgoPulls::Init()
+{
+ fL1 = CbmL1::Instance();
+
+ static bool first_call = 1;
+ if (first_call) {
+ // TDirectory *curdir = gDirectory;
+ // fL1->histodir->cd();
+ // gDirectory->mkdir("L1AlgoPulls");
+ // gDirectory->cd("L1AlgoPulls");
+
+ // add global pulls
+ for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
+ TString name = "pull_";
+ name += L1TrackParametersNames[i];
+ histoPull[i] = new TH1F(name, name, 50, -10, 10);
+ }
+
+#ifdef BUILD_HISTO_FOR_EACH_STANTION
+ // add station pulls
+ for (int i = TL1TrackParameters::NParameters; i < (NStations + 1) * TL1TrackParameters::NParameters; i++) {
+ int ista = i / TL1TrackParameters::NParameters - 1;
+ TString name = "pull_sta";
+ name += ista;
+ name += "_";
+ name += L1TrackParametersNames[i % TL1TrackParameters::NParameters];
+ histoPull[i] = new TH1F(name, name, 50, -10, 10);
+ }
+#endif
+
+ // add global residuals
+ for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
+ TString name = "residual_";
+ name += L1TrackParametersNames[i];
+ float size = 10;
+ switch (i) {
+ case 0: size = .01; break;
+ case 1: size = .01; break;
+ case 2: size = 0.003; break;
+ case 3: size = 0.003; break;
+ case 4: size = 0.1; break;
+ };
+ histoRes[i] = new TH1F(name, name, 50, -size, size);
+ }
+
+ // define style
+ histoStyle = new TStyle("histoStyle", "Plain Style(no colors/fill areas)");
+
+ histoStyle->SetTextFont(textFont);
+ histoStyle->SetPadColor(0);
+ histoStyle->SetCanvasColor(0);
+ histoStyle->SetTitleColor(0);
+ histoStyle->SetStatColor(0);
+
+ histoStyle->SetOptTitle(0); // without main up title
+ histoStyle->SetOptStat(1000001010);
+ // The parameter mode can be = ksiourmen (default = 000001111)
+ // k = 1; kurtosis printed
+ // k = 2; kurtosis and kurtosis error printed
+ // s = 1; skewness printed
+ // s = 2; skewness and skewness error printed
+ // i = 1; integral of bins printed
+ // o = 1; number of overflows printed
+ // u = 1; number of underflows printed
+ // r = 1; rms printed
+ // r = 2; rms and rms error printed
+ // m = 1; mean value printed
+ // m = 2; mean and mean error values printed
+ // e = 1; number of entries printed
+ // n = 1; name of histogram is printed
+
+ histoStyle->SetOptFit(10001);
+ // The parameter mode can be = pcev (default = 0111)
+ // p = 1; print Probability
+ // c = 1; print Chisquare/Number of degress of freedom
+ // e = 1; print errors (if e=1, v must be 1)
+ // v = 1; print name/values of parameters
+
+ histoStyle->SetStatW(0.175);
+ histoStyle->SetStatH(0.02);
+ histoStyle->SetStatX(0.95);
+ histoStyle->SetStatY(0.97);
+ histoStyle->SetStatFontSize(0.05);
+
+ histoStyle->SetStatFont(textFont);
+
+
+ // gDirectory->cd("..");
+ // curdir->cd();
+ first_call = 0;
+ }
+};
+
+// void L1AlgoPulls::AddVec(L1TrackPar& T_, THitI ih)
+// {
+// for (int i = 0; i < fvecLen; i++)
+// AddOne(T_,i,ih);
+// }
+
+void L1AlgoPulls::AddOne(L1TrackPar& T_, int i, THitI ih)
+{
+ fNAllPulls++;
+ TL1TrackParameters T(T_, i);
+
+ if (T_.chi2[i] > csCut * T_.NDF[i]) return;
+ // get error
+ TL1TrackParameters err;
+ if (!(finite(T_.C00[i]) && T_.C00[i] > 0)) return;
+ if (!(finite(T_.C11[i]) && T_.C11[i] > 0)) return;
+ if (!(finite(T_.C22[i]) && T_.C22[i] > 0)) return;
+ if (!(finite(T_.C33[i]) && T_.C33[i] > 0)) return;
+ if (!(finite(T_.C44[i]) && T_.C44[i] > 0)) return;
+ err.x = sqrt(T_.C00[i]);
+ err.y = sqrt(T_.C11[i]);
+ err.tx = sqrt(T_.C22[i]);
+ err.ty = sqrt(T_.C33[i]);
+ err.qp = sqrt(T_.C44[i]);
+
+ // mc data
+ int iMCP = fL1->vHitMCRef[ih];
+ if (iMCP < 0) return;
+ CbmL1MCPoint& mcP = fL1->vMCPoints[iMCP];
+ TL1TrackParameters mc(mcP);
+
+ // fill residuals
+ TL1TrackParameters res = (mc - T);
+ fGRes.push_back(res);
+
+ // fill pulls
+ TL1TrackParameters P = res / err;
+ fGPulls.push_back(P);
+
+#ifdef BUILD_HISTO_FOR_EACH_STANTION
+ int ista = mcP.iStation - 2; // last hit
+ // int ista = fL1->algo->vSFlag[ fL1->algo->vStsHits[ih].f ]/4 - 2; // last hit
+ fStaPulls[ista].push_back(P);
+#endif // BUILD_HISTO_FOR_EACH_STANTION
+};
+
+void L1AlgoPulls::Print()
+{ // TODO: renew
+ cout << "All pulls: " << fNAllPulls << endl;
+ cout << "Correct pulls: " << fGPulls.size() << endl;
+ cout << "x y tx ty qp" << endl;
+ for (int i = 0; i < fGPulls.size(); i++) {
+ TL1TrackParameters& pull = fGPulls[i];
+ pull.Print();
+ }
+};
+
+void L1AlgoPulls::Build(bool draw)
+{
+ // --- fill histograms ---
+ // global pulls
+ for (int i = 0; i < fGPulls.size(); i++) {
+ TL1TrackParameters& pull = fGPulls[i];
+ for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
+ if (TailCut > fabs(pull[ih])) histoPull[ih]->Fill(pull[ih]);
+ }
+ }
+#ifdef BUILD_HISTO_FOR_EACH_STANTION
+ // station pulls
+ for (int iSta = 0; iSta < NStations; iSta++) {
+ vector<TL1TrackParameters>& Pulls = fStaPulls[iSta];
+ for (int i = 0; i < Pulls.size(); i++) {
+ TL1TrackParameters& pull = Pulls[i];
+ for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
+ if (TailCut > fabs(pull[ih])) histoPull[(iSta + 1) * TL1TrackParameters::NParameters + ih]->Fill(pull[ih]);
+ }
+ }
+ }
+#endif // BUILD_HISTO_FOR_EACH_STANTION
+
+ // global residuals
+ for (int i = 0; i < fGRes.size(); i++) {
+ TL1TrackParameters& res = fGRes[i];
+ for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
+ if (TailCut > fabs(res[ih])) histoRes[ih]->Fill(res[ih]);
+ }
+ }
+
+ // --- draw histograms --- and save info
+ float pulls[(NStations + 1) * TL1TrackParameters::NParameters][2], // 0 - sigma, 1 - RMS
+ residuals[(NStations + 1) * TL1TrackParameters::NParameters][2];
+
+ system("mkdir L1_Pulls -p");
+ chdir("L1_Pulls");
+ TCanvas* c2 = new TCanvas("c2", "c2", 0, 0, 600, 400);
+ c2->cd();
+ histoStyle->cd();
+ for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters; ih++) {
+ makeUpHisto(histoPull[ih], histoPull[ih]->GetName(), pulls[ih][0]);
+ pulls[ih][1] = histoPull[ih]->GetRMS();
+ if (draw) {
+ histoPull[ih]->Draw();
+ TString name = histoPull[ih]->GetName();
+ name += ".pdf";
+ c2->SaveAs(name);
+ }
+ }
+ for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
+ makeUpHisto(histoRes[ih], histoRes[ih]->GetName(), residuals[ih][0]);
+ residuals[ih][1] = histoRes[ih]->GetRMS();
+ if (draw) {
+ histoRes[ih]->Draw();
+ TString name = histoRes[ih]->GetName();
+ name += ".pdf";
+ c2->SaveAs(name);
+ }
+ }
+ chdir("..");
+
+ // --- print information ---
+ cout << "All entries: " << fNAllPulls << endl;
+ cout << "Correct entries: " << fGPulls.size() << endl;
+ cout << "Pulls sigma & RMS: " << endl;
+ for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters; ih++) {
+ int ipar = ih % TL1TrackParameters::NParameters;
+ int ista = ih / TL1TrackParameters::NParameters;
+ if ((ista > 0) && (ipar == 0)) cout << "Station " << ista - 1 << endl;
+ cout << L1TrackParametersNames[ipar] << "\t" << pulls[ih][0] << "\t" << pulls[ih][1] << endl;
+ }
+ cout << "Residuals sigma & RMS: " << endl;
+ for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
+ int ipar = ih % TL1TrackParameters::NParameters;
+ cout << L1TrackParametersNames[ipar] << "\t" << residuals[ih][0] << "\t" << residuals[ih][1] << endl;
+ }
+}
+
+void L1AlgoPulls::makeUpHisto(TH1* hist, TString title, float& sigma)
+{
+ if (hist && (hist->GetEntries() != 0)) {
+ TF1* fit = new TF1("fit", "gaus");
+ fit->SetLineColor(2);
+ fit->SetLineWidth(3);
+ hist->Fit("fit", "", "", hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
+ sigma = fit->GetParameter(2);
+
+ hist->GetXaxis()->SetLabelFont(textFont);
+ hist->GetXaxis()->SetTitleFont(textFont);
+ hist->GetYaxis()->SetLabelFont(textFont);
+ hist->GetYaxis()->SetTitleFont(textFont);
+
+ hist->GetXaxis()->SetTitle(title);
+ hist->GetXaxis()->SetTitleOffset(1);
+ hist->GetYaxis()->SetTitle("Entries");
+ hist->GetYaxis()->SetTitleOffset(1.05); // good then entries per bit <= 9999
+ }
+ else {
+ std::cout << " E: Read hists error! " << std::endl;
+ }
+}
diff --git a/reco/L1/L1Algo/utils/L1AlgoPulls.h b/reco/L1/L1Algo/utils/L1AlgoPulls.h
index ee34a33631..46065b8b28 100644
--- a/reco/L1/L1Algo/utils/L1AlgoPulls.h
+++ b/reco/L1/L1Algo/utils/L1AlgoPulls.h
@@ -15,7 +15,6 @@ const int NStations = 0;
#include "CbmL1.h"
-#include "CbmL1Def.h"
#include "TCanvas.h"
#include "TF1.h"
@@ -25,6 +24,7 @@ const int NStations = 0;
#include <vector>
#include "L1Algo.h"
+#include "L1Def.h"
#include "L1TrackPar.h"
using std::cout;
@@ -32,9 +32,9 @@ using std::endl;
using std::vector;
struct TL1TrackParameters {
- double x, y, tx, ty, qp;
+ double x {0.}, y {0.}, tx {0.}, ty {0.}, qp {0.};
- static const int NParameters = 5;
+ static constexpr int NParameters = 5;
TL1TrackParameters() {};
TL1TrackParameters(L1TrackPar& T, int i) : x(T.x[i]), y(T.y[i]), tx(T.tx[i]), ty(T.ty[i]), qp(T.qp[i]) {};
@@ -48,6 +48,7 @@ struct TL1TrackParameters {
case 2: return tx;
case 3: return ty;
case 4: return qp;
+ default: return 0.f;
};
}
@@ -76,16 +77,10 @@ struct TL1TrackParameters {
void Print() { cout << x << " " << y << " " << tx << " " << ty << " " << qp << endl; }
};
-const TString L1TrackParametersNames[TL1TrackParameters::NParameters] = {"x", "y", "tx", "ty", "qp"};
-
class L1AlgoPulls {
public:
- L1AlgoPulls() : fNAllPulls(0)
- {
- fGPulls.clear();
- for (int i = 0; i < NStations; i++)
- fStaPulls[i].clear();
- };
+ L1AlgoPulls() = default;
+
void Init();
// void AddVec(L1TrackPar& T, THitI ih);
@@ -93,277 +88,27 @@ public:
void Print(); // fast method to see pulls :)
void Build(bool draw = 1);
- int fNAllPulls; // number of AddOne calls
-private:
- void makeUpHisto(TH1* hist, TString title, float& sigma);
-
- CbmL1* fL1;
- vector<TL1TrackParameters> fGPulls; // pulls for all stations
- vector<TL1TrackParameters> fStaPulls[NStations]; // pulls for each station
- TH1F* histoPull[(1 + NStations) * TL1TrackParameters::NParameters];
-
- vector<TL1TrackParameters> fGRes; // residuals for all stations
- TH1F* histoRes[(1 + NStations) * TL1TrackParameters::NParameters];
-
- static const float TailCut = 5000.; //
- static const float csCut = 5.; // chi-square cut
- static const int textFont = 22; // TNewRoman
- TStyle* histoStyle;
-};
-
-// ===================================================================================
-
-void L1AlgoPulls::Init()
-{
- fL1 = CbmL1::Instance();
-
-
- static bool first_call = 1;
- if (first_call) {
- // TDirectory *curdir = gDirectory;
- // fL1->histodir->cd();
- // gDirectory->mkdir("L1AlgoPulls");
- // gDirectory->cd("L1AlgoPulls");
-
- // add global pulls
- for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
- TString name = "pull_";
- name += L1TrackParametersNames[i];
- histoPull[i] = new TH1F(name, name, 50, -10, 10);
- }
-
-#ifdef BUILD_HISTO_FOR_EACH_STANTION
- // add station pulls
- for (int i = TL1TrackParameters::NParameters; i < (NStations + 1) * TL1TrackParameters::NParameters; i++) {
- int ista = i / TL1TrackParameters::NParameters - 1;
- TString name = "pull_sta";
- name += ista;
- name += "_";
- name += L1TrackParametersNames[i % TL1TrackParameters::NParameters];
- histoPull[i] = new TH1F(name, name, 50, -10, 10);
- }
-#endif // BUILD_HISTO_FOR_EACH_STANTION \
- // add global residuals
- for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
- TString name = "residual_";
- name += L1TrackParametersNames[i];
- float size = 10;
- switch (i) {
- case 0: size = .01; break;
- case 1: size = .01; break;
- case 2: size = 0.003; break;
- case 3: size = 0.003; break;
- case 4: size = 0.1; break;
- };
- histoRes[i] = new TH1F(name, name, 50, -size, size);
- }
-
- // define style
- histoStyle = new TStyle("histoStyle", "Plain Style(no colors/fill areas)");
-
- histoStyle->SetTextFont(textFont);
- histoStyle->SetPadColor(0);
- histoStyle->SetCanvasColor(0);
- histoStyle->SetTitleColor(0);
- histoStyle->SetStatColor(0);
-
- histoStyle->SetOptTitle(0); // without main up title
- histoStyle->SetOptStat(1000001010);
- // The parameter mode can be = ksiourmen (default = 000001111)
- // k = 1; kurtosis printed
- // k = 2; kurtosis and kurtosis error printed
- // s = 1; skewness printed
- // s = 2; skewness and skewness error printed
- // i = 1; integral of bins printed
- // o = 1; number of overflows printed
- // u = 1; number of underflows printed
- // r = 1; rms printed
- // r = 2; rms and rms error printed
- // m = 1; mean value printed
- // m = 2; mean and mean error values printed
- // e = 1; number of entries printed
- // n = 1; name of histogram is printed
-
- histoStyle->SetOptFit(10001);
- // The parameter mode can be = pcev (default = 0111)
- // p = 1; print Probability
- // c = 1; print Chisquare/Number of degress of freedom
- // e = 1; print errors (if e=1, v must be 1)
- // v = 1; print name/values of parameters
-
- histoStyle->SetStatW(0.175);
- histoStyle->SetStatH(0.02);
- histoStyle->SetStatX(0.95);
- histoStyle->SetStatY(0.97);
- histoStyle->SetStatFontSize(0.05);
-
- histoStyle->SetStatFont(textFont);
-
-
- // gDirectory->cd("..");
- // curdir->cd();
- first_call = 0;
- }
-};
-
-// inline void L1AlgoPulls::AddVec(L1TrackPar& T_, THitI ih)
-// {
-// for (int i = 0; i < fvecLen; i++)
-// AddOne(T_,i,ih);
-// }
+ int fNAllPulls {0}; // number of AddOne calls
-inline void L1AlgoPulls::AddOne(L1TrackPar& T_, int i, THitI ih)
-{
- fNAllPulls++;
- TL1TrackParameters T(T_, i);
-
- if (T_.chi2[i] > csCut * T_.NDF[i]) return;
- // get error
- TL1TrackParameters err;
- if (!(finite(T_.C00[i]) && T_.C00[i] > 0)) return;
- if (!(finite(T_.C11[i]) && T_.C11[i] > 0)) return;
- if (!(finite(T_.C22[i]) && T_.C22[i] > 0)) return;
- if (!(finite(T_.C33[i]) && T_.C33[i] > 0)) return;
- if (!(finite(T_.C44[i]) && T_.C44[i] > 0)) return;
- err.x = sqrt(T_.C00[i]);
- err.y = sqrt(T_.C11[i]);
- err.tx = sqrt(T_.C22[i]);
- err.ty = sqrt(T_.C33[i]);
- err.qp = sqrt(T_.C44[i]);
-
- // mc data
- int iMCP = fL1->vHitMCRef[ih];
- if (iMCP < 0) return;
- CbmL1MCPoint& mcP = fL1->vMCPoints[iMCP];
- TL1TrackParameters mc(mcP);
-
- // fill residuals
- TL1TrackParameters res = (mc - T);
- fGRes.push_back(res);
-
- // fill pulls
- TL1TrackParameters P = res / err;
- fGPulls.push_back(P);
-
-#ifdef BUILD_HISTO_FOR_EACH_STANTION
- int ista = mcP.iStation - 2; // last hit
- // int ista = fL1->algo->vSFlag[ fL1->algo->vStsHits[ih].f ]/4 - 2; // last hit
- fStaPulls[ista].push_back(P);
-#endif // BUILD_HISTO_FOR_EACH_STANTION
-};
-
-inline void L1AlgoPulls::Print()
-{ // TODO: renew
- cout << "All pulls: " << fNAllPulls << endl;
- cout << "Correct pulls: " << fGPulls.size() << endl;
- cout << "x y tx ty qp" << endl;
- for (int i = 0; i < fGPulls.size(); i++) {
- TL1TrackParameters& pull = fGPulls[i];
- pull.Print();
- }
-};
-
-inline void L1AlgoPulls::Build(bool draw)
-{
- // --- fill histograms ---
- // global pulls
- for (int i = 0; i < fGPulls.size(); i++) {
- TL1TrackParameters& pull = fGPulls[i];
- for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
- if (TailCut > fabs(pull[ih])) histoPull[ih]->Fill(pull[ih]);
- }
- }
-#ifdef BUILD_HISTO_FOR_EACH_STANTION
- // station pulls
- for (int iSta = 0; iSta < NStations; iSta++) {
- vector<TL1TrackParameters>& Pulls = fStaPulls[iSta];
- for (int i = 0; i < Pulls.size(); i++) {
- TL1TrackParameters& pull = Pulls[i];
- for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
- if (TailCut > fabs(pull[ih])) histoPull[(iSta + 1) * TL1TrackParameters::NParameters + ih]->Fill(pull[ih]);
- }
- }
- }
-#endif // BUILD_HISTO_FOR_EACH_STANTION
+private:
+ L1AlgoPulls(const L1AlgoPulls&);
+ L1AlgoPulls& operator=(const L1AlgoPulls&);
- // global residuals
- for (int i = 0; i < fGRes.size(); i++) {
- TL1TrackParameters& res = fGRes[i];
- for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
- if (TailCut > fabs(res[ih])) histoRes[ih]->Fill(res[ih]);
- }
- }
+ void makeUpHisto(TH1* hist, TString title, float& sigma);
- // --- draw histograms --- and save info
- float pulls[(NStations + 1) * TL1TrackParameters::NParameters][2], // 0 - sigma, 1 - RMS
- residuals[(NStations + 1) * TL1TrackParameters::NParameters][2];
+ CbmL1* fL1 {nullptr};
+ vector<TL1TrackParameters> fGPulls {}; // pulls for all stations
+ vector<TL1TrackParameters> fStaPulls[NStations] {}; // pulls for each station
+ TH1F* histoPull[(1 + NStations) * TL1TrackParameters::NParameters] {nullptr};
- system("mkdir L1_Pulls -p");
- chdir("L1_Pulls");
- TCanvas* c2 = new TCanvas("c2", "c2", 0, 0, 600, 400);
- c2->cd();
- histoStyle->cd();
- for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters; ih++) {
- makeUpHisto(histoPull[ih], histoPull[ih]->GetName(), pulls[ih][0]);
- pulls[ih][1] = histoPull[ih]->GetRMS();
- if (draw) {
- histoPull[ih]->Draw();
- TString name = histoPull[ih]->GetName();
- name += ".pdf";
- c2->SaveAs(name);
- }
- }
- for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
- makeUpHisto(histoRes[ih], histoRes[ih]->GetName(), residuals[ih][0]);
- residuals[ih][1] = histoRes[ih]->GetRMS();
- if (draw) {
- histoRes[ih]->Draw();
- TString name = histoRes[ih]->GetName();
- name += ".pdf";
- c2->SaveAs(name);
- }
- }
- chdir("..");
+ vector<TL1TrackParameters> fGRes {}; // residuals for all stations
+ TH1F* histoRes[(1 + NStations) * TL1TrackParameters::NParameters] {nullptr};
- // --- print information ---
- cout << "All entries: " << fNAllPulls << endl;
- cout << "Correct entries: " << fGPulls.size() << endl;
- cout << "Pulls sigma & RMS: " << endl;
- for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters; ih++) {
- int ipar = ih % TL1TrackParameters::NParameters;
- int ista = ih / TL1TrackParameters::NParameters;
- if ((ista > 0) && (ipar == 0)) cout << "Station " << ista - 1 << endl;
- cout << L1TrackParametersNames[ipar] << "\t" << pulls[ih][0] << "\t" << pulls[ih][1] << endl;
- }
- cout << "Residuals sigma & RMS: " << endl;
- for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
- int ipar = ih % TL1TrackParameters::NParameters;
- cout << L1TrackParametersNames[ipar] << "\t" << residuals[ih][0] << "\t" << residuals[ih][1] << endl;
- }
+ static constexpr float TailCut = 5000.; //
+ static constexpr float csCut = 5.; // chi-square cut
+ static constexpr int textFont = 22; // TNewRoman
+ TStyle* histoStyle {nullptr};
};
-inline void L1AlgoPulls::makeUpHisto(TH1* hist, TString title, float& sigma)
-{
- if (hist && (hist->GetEntries() != 0)) {
- TF1* fit = new TF1("fit", "gaus");
- fit->SetLineColor(2);
- fit->SetLineWidth(3);
- hist->Fit("fit", "", "", hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
- sigma = fit->GetParameter(2);
-
- hist->GetXaxis()->SetLabelFont(textFont);
- hist->GetXaxis()->SetTitleFont(textFont);
- hist->GetYaxis()->SetLabelFont(textFont);
- hist->GetYaxis()->SetTitleFont(textFont);
-
- hist->GetXaxis()->SetTitle(title);
- hist->GetXaxis()->SetTitleOffset(1);
- hist->GetYaxis()->SetTitle("Entries");
- hist->GetYaxis()->SetTitleOffset(1.05); // good then entries per bit <= 9999
- }
- else {
- std::cout << " E: Read hists error! " << std::endl;
- }
-}
#endif
--
GitLab