CbmRoot
L1AlgoPulls.h
Go to the documentation of this file.
1 #ifndef L1AlgoPulls_h
2 #define L1AlgoPulls_h
3 
4 // #define BUILD_HISTO_FOR_EACH_STANTION
5 
6 #ifdef BUILD_HISTO_FOR_EACH_STANTION
7 const int NStations = 8;
8 #else
9 const int NStations = 0;
10 #endif // BUILD_HISTO_FOR_EACH_STANTION
11 
12 
13 #include "TCanvas.h"
14 #include "TF1.h"
15 #include "TStyle.h"
16 
17 #include "CbmL1.h"
18 #include "CbmL1Def.h"
19 #include "L1Algo/L1Algo.h"
20 #include "L1Algo/L1StsHit.h"
21 #include "L1Algo/L1TrackPar.h"
22 #include <iostream>
23 #include <vector>
24 
25 using std::cout;
26 using std::endl;
27 using std::vector;
28 
29 struct TL1TrackParameters {
30  double x, y, tx, ty, qp;
31 
32  static const int NParameters = 5;
33 
34  TL1TrackParameters() {};
35  TL1TrackParameters(L1TrackPar& T, int i)
36  : x(T.x[i]), y(T.y[i]), tx(T.tx[i]), ty(T.ty[i]), qp(T.qp[i]) {};
37  TL1TrackParameters(CbmL1MCPoint& T)
38  : x(T.x), y(T.y), tx(T.px / T.pz), ty(T.py / T.pz), qp(T.q / T.p) {};
39 
40  double operator[](int i) {
41  switch (i) {
42  case 0: return x;
43  case 1: return y;
44  case 2: return tx;
45  case 3: return ty;
46  case 4: return qp;
47  };
48  }
49 
50  TL1TrackParameters operator-(TL1TrackParameters& b) {
51  TL1TrackParameters c;
52  c.x = x - b.x;
53  c.y = y - b.y;
54  c.tx = tx - b.tx;
55  c.ty = ty - b.ty;
56  c.qp = qp - b.qp;
57  return c;
58  }
59 
60  TL1TrackParameters operator/(TL1TrackParameters& b) {
61  TL1TrackParameters c;
62  c.x = x / b.x;
63  c.y = y / b.y;
64  c.tx = tx / b.tx;
65  c.ty = ty / b.ty;
66  c.qp = qp / b.qp;
67  return c;
68  }
69 
70  void Print() {
71  cout << x << " " << y << " " << tx << " " << ty << " " << qp << endl;
72  }
73 };
74 
75 const TString L1TrackParametersNames[TL1TrackParameters::NParameters] = {"x",
76  "y",
77  "tx",
78  "ty",
79  "qp"};
80 
81 class L1AlgoPulls {
82 public:
83  L1AlgoPulls() : fNAllPulls(0) {
84  fGPulls.clear();
85  for (int i = 0; i < NStations; i++)
86  fStaPulls[i].clear();
87  };
88  void Init();
89 
90  // void AddVec(L1TrackPar& T, THitI ih);
91  void AddOne(L1TrackPar& T, int i, THitI ih);
92  void Print(); // fast method to see pulls :)
93  void Build(bool draw = 1);
94 
95  int fNAllPulls; // number of AddOne calls
96 private:
97  void makeUpHisto(TH1* hist, TString title, float& sigma);
98 
99  CbmL1* fL1;
100  vector<TL1TrackParameters> fGPulls; // pulls for all stations
101  vector<TL1TrackParameters> fStaPulls[NStations]; // pulls for each station
102  TH1F* histoPull[(1 + NStations) * TL1TrackParameters::NParameters];
103 
104  vector<TL1TrackParameters> fGRes; // residuals for all stations
105  TH1F* histoRes[(1 + NStations) * TL1TrackParameters::NParameters];
106 
107  static const float TailCut = 5000.; //
108  static const float csCut = 5.; // chi-square cut
109  static const int textFont = 22; // TNewRoman
110  TStyle* histoStyle;
111 };
112 
113 // ===================================================================================
114 
115 void L1AlgoPulls::Init() {
116  fL1 = CbmL1::Instance();
117 
118 
119  static bool first_call = 1;
120  if (first_call) {
121  // TDirectory *curdir = gDirectory;
122  // fL1->histodir->cd();
123  // gDirectory->mkdir("L1AlgoPulls");
124  // gDirectory->cd("L1AlgoPulls");
125 
126  // add global pulls
127  for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
128  TString name = "pull_";
129  name += L1TrackParametersNames[i];
130  histoPull[i] = new TH1F(name, name, 50, -10, 10);
131  }
132 
133 #ifdef BUILD_HISTO_FOR_EACH_STANTION
134  // add station pulls
135  for (int i = TL1TrackParameters::NParameters;
136  i < (NStations + 1) * TL1TrackParameters::NParameters;
137  i++) {
138  int ista = i / TL1TrackParameters::NParameters - 1;
139  TString name = "pull_sta";
140  name += ista;
141  name += "_";
142  name += L1TrackParametersNames[i % TL1TrackParameters::NParameters];
143  histoPull[i] = new TH1F(name, name, 50, -10, 10);
144  }
145 #endif // BUILD_HISTO_FOR_EACH_STANTION \
146  // add global residuals
147  for (int i = 0; i < TL1TrackParameters::NParameters; i++) {
148  TString name = "residual_";
149  name += L1TrackParametersNames[i];
150  float size = 10;
151  switch (i) {
152  case 0: size = .01; break;
153  case 1: size = .01; break;
154  case 2: size = 0.003; break;
155  case 3: size = 0.003; break;
156  case 4: size = 0.1; break;
157  };
158  histoRes[i] = new TH1F(name, name, 50, -size, size);
159  }
160 
161  // define style
162  histoStyle = new TStyle("histoStyle", "Plain Style(no colors/fill areas)");
163 
164  histoStyle->SetTextFont(textFont);
165  histoStyle->SetPadColor(0);
166  histoStyle->SetCanvasColor(0);
167  histoStyle->SetTitleColor(0);
168  histoStyle->SetStatColor(0);
169 
170  histoStyle->SetOptTitle(0); // without main up title
171  histoStyle->SetOptStat(1000001010);
172  // The parameter mode can be = ksiourmen (default = 000001111)
173  // k = 1; kurtosis printed
174  // k = 2; kurtosis and kurtosis error printed
175  // s = 1; skewness printed
176  // s = 2; skewness and skewness error printed
177  // i = 1; integral of bins printed
178  // o = 1; number of overflows printed
179  // u = 1; number of underflows printed
180  // r = 1; rms printed
181  // r = 2; rms and rms error printed
182  // m = 1; mean value printed
183  // m = 2; mean and mean error values printed
184  // e = 1; number of entries printed
185  // n = 1; name of histogram is printed
186 
187  histoStyle->SetOptFit(10001);
188  // The parameter mode can be = pcev (default = 0111)
189  // p = 1; print Probability
190  // c = 1; print Chisquare/Number of degress of freedom
191  // e = 1; print errors (if e=1, v must be 1)
192  // v = 1; print name/values of parameters
193 
194  histoStyle->SetStatW(0.175);
195  histoStyle->SetStatH(0.02);
196  histoStyle->SetStatX(0.95);
197  histoStyle->SetStatY(0.97);
198  histoStyle->SetStatFontSize(0.05);
199 
200  histoStyle->SetStatFont(textFont);
201 
202 
203  // gDirectory->cd("..");
204  // curdir->cd();
205  first_call = 0;
206  }
207 };
208 
209 // inline void L1AlgoPulls::AddVec(L1TrackPar& T_, THitI ih)
210 // {
211 // for (int i = 0; i < fvecLen; i++)
212 // AddOne(T_,i,ih);
213 // }
214 
215 inline void L1AlgoPulls::AddOne(L1TrackPar& T_, int i, THitI ih) {
216  fNAllPulls++;
217  TL1TrackParameters T(T_, i);
218 
219  if (T_.chi2[i] > csCut * T_.NDF[i]) return;
220  // get error
221  TL1TrackParameters err;
222  if (!(finite(T_.C00[i]) && T_.C00[i] > 0)) return;
223  if (!(finite(T_.C11[i]) && T_.C11[i] > 0)) return;
224  if (!(finite(T_.C22[i]) && T_.C22[i] > 0)) return;
225  if (!(finite(T_.C33[i]) && T_.C33[i] > 0)) return;
226  if (!(finite(T_.C44[i]) && T_.C44[i] > 0)) return;
227  err.x = sqrt(T_.C00[i]);
228  err.y = sqrt(T_.C11[i]);
229  err.tx = sqrt(T_.C22[i]);
230  err.ty = sqrt(T_.C33[i]);
231  err.qp = sqrt(T_.C44[i]);
232 
233  // mc data
234  int iMCP = fL1->vHitMCRef[ih];
235  if (iMCP < 0) return;
236  CbmL1MCPoint& mcP = fL1->vMCPoints[iMCP];
237  TL1TrackParameters mc(mcP);
238 
239  // fill residuals
240  TL1TrackParameters res = (mc - T);
241  fGRes.push_back(res);
242 
243  // fill pulls
244  TL1TrackParameters P = res / err;
245  fGPulls.push_back(P);
246 
247 #ifdef BUILD_HISTO_FOR_EACH_STANTION
248  int ista = mcP.iStation - 2; // last hit
249  // int ista = fL1->algo->vSFlag[ fL1->algo->vStsHits[ih].f ]/4 - 2; // last hit
250  fStaPulls[ista].push_back(P);
251 #endif // BUILD_HISTO_FOR_EACH_STANTION
252 };
253 
254 inline void L1AlgoPulls::Print() { // TODO: renew
255  cout << "All pulls: " << fNAllPulls << endl;
256  cout << "Correct pulls: " << fGPulls.size() << endl;
257  cout << "x y tx ty qp" << endl;
258  for (int i = 0; i < fGPulls.size(); i++) {
259  TL1TrackParameters& pull = fGPulls[i];
260  pull.Print();
261  }
262 };
263 
264 inline void L1AlgoPulls::Build(bool draw) {
265  // --- fill histograms ---
266  // global pulls
267  for (int i = 0; i < fGPulls.size(); i++) {
268  TL1TrackParameters& pull = fGPulls[i];
269  for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
270  if (TailCut > fabs(pull[ih])) histoPull[ih]->Fill(pull[ih]);
271  }
272  }
273 #ifdef BUILD_HISTO_FOR_EACH_STANTION
274  // station pulls
275  for (int iSta = 0; iSta < NStations; iSta++) {
276  vector<TL1TrackParameters>& Pulls = fStaPulls[iSta];
277  for (int i = 0; i < Pulls.size(); i++) {
278  TL1TrackParameters& pull = Pulls[i];
279  for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
280  if (TailCut > fabs(pull[ih]))
281  histoPull[(iSta + 1) * TL1TrackParameters::NParameters + ih]->Fill(
282  pull[ih]);
283  }
284  }
285  }
286 #endif // BUILD_HISTO_FOR_EACH_STANTION
287 
288  // global residuals
289  for (int i = 0; i < fGRes.size(); i++) {
290  TL1TrackParameters& res = fGRes[i];
291  for (int ih = 0; ih < TL1TrackParameters::NParameters; ih++) {
292  if (TailCut > fabs(res[ih])) histoRes[ih]->Fill(res[ih]);
293  }
294  }
295 
296  // --- draw histograms --- and save info
297  float pulls[(NStations + 1) * TL1TrackParameters::NParameters]
298  [2], // 0 - sigma, 1 - RMS
299  residuals[(NStations + 1) * TL1TrackParameters::NParameters][2];
300 
301  system("mkdir L1_Pulls -p");
302  chdir("L1_Pulls");
303  TCanvas* c2 = new TCanvas("c2", "c2", 0, 0, 600, 400);
304  c2->cd();
305  histoStyle->cd();
306  for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters;
307  ih++) {
308  makeUpHisto(histoPull[ih], histoPull[ih]->GetName(), pulls[ih][0]);
309  pulls[ih][1] = histoPull[ih]->GetRMS();
310  if (draw) {
311  histoPull[ih]->Draw();
312  TString name = histoPull[ih]->GetName();
313  name += ".pdf";
314  c2->SaveAs(name);
315  }
316  }
317  for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
318  makeUpHisto(histoRes[ih], histoRes[ih]->GetName(), residuals[ih][0]);
319  residuals[ih][1] = histoRes[ih]->GetRMS();
320  if (draw) {
321  histoRes[ih]->Draw();
322  TString name = histoRes[ih]->GetName();
323  name += ".pdf";
324  c2->SaveAs(name);
325  }
326  }
327  chdir("..");
328 
329  // --- print information ---
330  cout << "All entries: " << fNAllPulls << endl;
331  cout << "Correct entries: " << fGPulls.size() << endl;
332  cout << "Pulls sigma & RMS: " << endl;
333  for (int ih = 0; ih < (NStations + 1) * TL1TrackParameters::NParameters;
334  ih++) {
335  int ipar = ih % TL1TrackParameters::NParameters;
336  int ista = ih / TL1TrackParameters::NParameters;
337  if ((ista > 0) && (ipar == 0)) cout << "Station " << ista - 1 << endl;
338  cout << L1TrackParametersNames[ipar] << "\t" << pulls[ih][0] << "\t"
339  << pulls[ih][1] << endl;
340  }
341  cout << "Residuals sigma & RMS: " << endl;
342  for (int ih = 0; ih < (0 + 1) * TL1TrackParameters::NParameters; ih++) {
343  int ipar = ih % TL1TrackParameters::NParameters;
344  cout << L1TrackParametersNames[ipar] << "\t" << residuals[ih][0] << "\t"
345  << residuals[ih][1] << endl;
346  }
347 };
348 
349 inline void L1AlgoPulls::makeUpHisto(TH1* hist, TString title, float& sigma) {
350  if (hist && (hist->GetEntries() != 0)) {
351  TF1* fit = new TF1("fit", "gaus");
352  fit->SetLineColor(2);
353  fit->SetLineWidth(3);
354  hist->Fit(
355  "fit", "", "", hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
356  sigma = fit->GetParameter(2);
357 
358  hist->GetXaxis()->SetLabelFont(textFont);
359  hist->GetXaxis()->SetTitleFont(textFont);
360  hist->GetYaxis()->SetLabelFont(textFont);
361  hist->GetYaxis()->SetTitleFont(textFont);
362 
363  hist->GetXaxis()->SetTitle(title);
364  hist->GetXaxis()->SetTitleOffset(1);
365  hist->GetYaxis()->SetTitle("Entries");
366  hist->GetYaxis()->SetTitleOffset(
367  1.05); // good then entries per bit <= 9999
368  } else {
369  std::cout << " E: Read hists error! " << std::endl;
370  }
371 }
372 
373 #endif
L1AlgoPulls::fNAllPulls
int fNAllPulls
Definition: L1AlgoPulls.h:95
L1AlgoPulls::textFont
static const int textFont
Definition: L1AlgoPulls.h:109
TL1TrackParameters::x
double x
Definition: L1AlgoPulls.h:30
L1Algo.h
TL1TrackParameters::qp
double qp
Definition: L1AlgoPulls.h:30
CbmL1MCPoint
Definition: CbmL1MCPoint.h:23
L1AlgoPulls::histoPull
TH1F * histoPull[(1+NStations) *TL1TrackParameters::NParameters]
Definition: L1AlgoPulls.h:102
TL1TrackParameters::TL1TrackParameters
TL1TrackParameters()
Definition: L1AlgoPulls.h:34
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
L1AlgoPulls
Definition: L1AlgoPulls.h:81
L1AlgoPulls::histoStyle
TStyle * histoStyle
Definition: L1AlgoPulls.h:110
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
L1AlgoPulls::fStaPulls
vector< TL1TrackParameters > fStaPulls[NStations]
Definition: L1AlgoPulls.h:101
L1AlgoPulls::makeUpHisto
void makeUpHisto(TH1 *hist, TString title, float &sigma)
Definition: L1AlgoPulls.h:349
L1AlgoPulls::csCut
static const float csCut
Definition: L1AlgoPulls.h:108
TL1TrackParameters::Print
void Print()
Definition: L1AlgoPulls.h:70
CbmL1.h
TL1TrackParameters::y
double y
Definition: L1AlgoPulls.h:30
NStations
const int NStations
Definition: L1AlgoPulls.h:9
TL1TrackParameters::operator[]
double operator[](int i)
Definition: L1AlgoPulls.h:40
Pulls
void Pulls(int i, int j, int k, double *mc, L1TrackPar &T, fvec qp0, L1FieldRegion &fld)
Definition: L1CADebug.h:15
L1AlgoPulls::fL1
CbmL1 * fL1
Definition: L1AlgoPulls.h:99
L1AlgoPulls::fGPulls
vector< TL1TrackParameters > fGPulls
Definition: L1AlgoPulls.h:100
CbmL1Def.h
THitI
unsigned int THitI
Definition: L1StsHit.h:8
finite
T finite(T x)
Definition: CbmL1Def.h:21
L1TrackPar::C44
fvec C44
Definition: L1TrackPar.h:10
L1TrackPar::chi2
fvec chi2
Definition: L1TrackPar.h:10
CbmL1::Instance
static CbmL1 * Instance()
Definition: CbmL1.h:129
L1TrackPar::NDF
fvec NDF
Definition: L1TrackPar.h:10
L1AlgoPulls::fGRes
vector< TL1TrackParameters > fGRes
Definition: L1AlgoPulls.h:104
TL1TrackParameters::TL1TrackParameters
TL1TrackParameters(CbmL1MCPoint &T)
Definition: L1AlgoPulls.h:37
L1AlgoPulls::AddOne
void AddOne(L1TrackPar &T, int i, THitI ih)
Definition: L1AlgoPulls.h:215
CbmL1MCPoint::iStation
int iStation
Definition: CbmL1MCPoint.h:61
L1TrackPar.h
TL1TrackParameters::ty
double ty
Definition: L1AlgoPulls.h:30
L1TrackPar::C00
fvec C00
Definition: L1TrackPar.h:9
TL1TrackParameters::operator-
TL1TrackParameters operator-(TL1TrackParameters &b)
Definition: L1AlgoPulls.h:50
L1AlgoPulls::histoRes
TH1F * histoRes[(1+NStations) *TL1TrackParameters::NParameters]
Definition: L1AlgoPulls.h:105
L1TrackPar::C33
fvec C33
Definition: L1TrackPar.h:9
L1AlgoPulls::Build
void Build(bool draw=1)
Definition: L1AlgoPulls.h:264
L1AlgoPulls::Print
void Print()
Definition: L1AlgoPulls.h:254
TL1TrackParameters::NParameters
static const int NParameters
Definition: L1AlgoPulls.h:32
L1TrackParametersNames
const TString L1TrackParametersNames[TL1TrackParameters::NParameters]
Definition: L1AlgoPulls.h:75
L1TrackPar
Definition: L1TrackPar.h:6
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
L1AlgoPulls::Init
void Init()
Definition: L1AlgoPulls.h:115
L1StsHit.h
L1TrackPar::C11
fvec C11
Definition: L1TrackPar.h:9
CbmL1::vHitMCRef
vector< int > vHitMCRef
Definition: CbmL1.h:339
CbmL1::vMCPoints
vector< CbmL1MCPoint > vMCPoints
Definition: CbmL1.h:197
L1AlgoPulls::L1AlgoPulls
L1AlgoPulls()
Definition: L1AlgoPulls.h:83
L1AlgoPulls::TailCut
static const float TailCut
Definition: L1AlgoPulls.h:107
CbmL1
Definition: CbmL1.h:113
PairAnalysisStyler::Fill
static Int_t Fill[]
Definition: PairAnalysisStyleDefs.h:82
TL1TrackParameters::operator/
TL1TrackParameters operator/(TL1TrackParameters &b)
Definition: L1AlgoPulls.h:60
TL1TrackParameters::tx
double tx
Definition: L1AlgoPulls.h:30
TL1TrackParameters::TL1TrackParameters
TL1TrackParameters(L1TrackPar &T, int i)
Definition: L1AlgoPulls.h:35
L1TrackPar::C22
fvec C22
Definition: L1TrackPar.h:9
TL1TrackParameters
Definition: L1AlgoPulls.h:29