Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawTpcPrototypeGenFitTrkFitter.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawTpcPrototypeGenFitTrkFitter.C
1 #include "SaveCanvas.C"
2 #include "sPhenixStyle.C"
3 
4 #include <TChain.h>
5 #include <TCut.h>
6 #include <TEfficiency.h>
7 #include <TF1.h>
8 #include <TGraphAsymmErrors.h>
9 #include <TGraphErrors.h>
10 #include <TH2.h>
11 #include <TH3.h>
12 #include <TPolyLine.h>
13 #include <TTree.h>
14 
15 #include <TFile.h>
16 
17 #include <TColor.h>
18 #include <TLatex.h>
19 #include <TLegend.h>
20 #include <TLine.h>
21 #include <TStyle.h>
22 
23 #include <TDirectory.h>
24 #include <TMath.h>
25 #include <TPad.h>
26 #include <TString.h>
27 #include <TTree.h>
28 #include <TVirtualFitter.h>
29 
30 #include <cmath>
31 #include <iostream>
32 
33 using namespace std;
34 
35 // ROOT6 disabled assert. Well....
36 #ifdef assert
37 #undef assert
38 #endif
39 #define assert(exp) \
40  { \
41  if (!(exp)) \
42  { \
43  cout << "Assert (" << #exp << ") failed at " << __FILE__ << " line " << __LINE__ << endl; \
44  exit(1); \
45  } \
46  }
47 
49 void useLogBins(TAxis *axis)
50 {
51  assert(axis);
52  assert(axis->GetXmin() > 0);
53  assert(axis->GetXmax() > 0);
54 
55  const int bins = axis->GetNbins();
56 
57  Axis_t from = log10(axis->GetXmin());
58  Axis_t to = log10(axis->GetXmax());
59  Axis_t width = (to - from) / bins;
60  vector<Axis_t> new_bins(bins + 1);
61 
62  for (int i = 0; i <= bins; i++)
63  {
64  new_bins[i] = TMath::Power(10, from + i * width);
65  }
66  axis->Set(bins, new_bins.data());
67 }
68 
69 TGraphErrors *
70 FitProfile(const TH2 *h2)
71 {
72  TProfile *p2 = h2->ProfileX();
73 
74  int n = 0;
75  double x[1000];
76  double ex[1000];
77  double y[1000];
78  double ey[1000];
79 
80  for (int i = 1; i <= h2->GetNbinsX(); i++)
81  {
82  TH1D *h1 = h2->ProjectionY(Form("htmp_%d", rand()), i, i);
83 
84  if (h1->GetSum() < 30)
85  {
86  cout << "FitProfile - ignore bin " << i << endl;
87  continue;
88  }
89  else
90  {
91  // cout << "FitProfile - fit bin " << i << endl;
92  }
93 
94  TF1 fgaus("fgaus", "gaus", -p2->GetBinError(i) * 4,
95  p2->GetBinError(i) * 4);
96 
97  // TF1 f2(Form("dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
98  // -p2->GetBinError(i) * 4, p2->GetBinError(i) * 4);
99 
100  fgaus.SetParameter(1, p2->GetBinContent(i));
101  fgaus.SetParameter(2, p2->GetBinError(i));
102 
103  h1->Fit(&fgaus, "MQ");
104 
105  TF1 f2(Form("dgaus"), "gaus + [3]",
106  -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
107 
108  f2.SetParameters(fgaus.GetParameter(0), fgaus.GetParameter(1),
109  fgaus.GetParameter(2), 0);
110 
111  h1->Fit(&f2, "MQR");
112 
113  // new TCanvas;
114  // h1->Draw();
115  // fgaus.Draw("same");
116  // break;
117 
118  x[n] = p2->GetBinCenter(i);
119  ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
120  y[n] = fgaus.GetParameter(1);
121  ey[n] = fgaus.GetParError(1);
122 
123  // p2->SetBinContent(i, fgaus.GetParameter(1));
124  // p2->SetBinError(i, fgaus.GetParameter(2));
125 
126  n++;
127  delete h1;
128  }
129 
130  return new TGraphErrors(n, x, y, ex, ey);
131 }
132 
133 TFile *_file0 = NULL;
134 TString description;
135 TTree *T(nullptr);
136 
137 void Resolution(const TCut &cut = "Iteration$ >= 5 && Iteration$ <= 10 && TPCTrack.nCluster>=12 && TPCTrack.clusterSizePhi > 3.5",
138  const double phi_start = -2.905, const double phi_end = -2.885
139 
140 )
141 {
142  TH2 *hresidual_phi = new TH2F("hresidual_phi", "hresidual_phi", 60, phi_start, phi_end, 60, -.2, .2);
143  T->Draw("TPCTrack.clusterResidualPhi:TPCTrack.clusterProjectionPhi>>hresidual_phi",
144  cut, "goff");
145  hresidual_phi->SetTitle(";Global Phi [rad];Phi Residual [cm]");
146 
147  TH2 *hresidual_phi_highres = new TH2F("hresidual_phi_highres", "hresidual_phi", 500, phi_start, phi_end, 1000, -.2, .2);
148  T->Draw("TPCTrack.clusterResidualPhi:TPCTrack.clusterProjectionPhi>>hresidual_phi_highres",
149  cut, "goff");
150  hresidual_phi_highres->SetTitle(";Global Phi [rad];Phi Residual [cm]");
151 
152  TCanvas *c1 = new TCanvas("Resolution", "Resolution", 1800, 860);
153  c1->Divide(3, 1);
154  int idx = 1;
155  TPad *p = nullptr;
156 
157  p = (TPad *) c1->cd(idx++);
158  c1->Update();
159 
160  TGraphErrors *gcenter = FitProfile(hresidual_phi);
161 
162  hresidual_phi->Draw("colz");
163  gcenter->Draw("p");
164  // hresidual_phi
165 
166  // TF1 *fPeriod = new TF1("fPeriod",
167  // "[0] * sin(2*pi*x/(2*pi/12/8/16)) + [1] * sin(2*pi*x/(2*pi/12/8/16)) + [2] + [3]*x + [4]*x*x + [5] * sin(4*pi*x/(2*pi/12/8/16)) + [6] * sin(4*pi*x/(2*pi/12/8/16))+ [7] * sin(3*pi*x/(2*pi/12/8/16)) + [8] * sin(3*pi*x/(2*pi/12/8/16))+ [9] * sin(pi*x/(5*pi/12/8/16)) + [10] * sin(pi*x/(5*pi/12/8/16)) + [11] * x*x*x", -3, 0);
168  TF1 *fPeriod = new TF1("fPeriod",
169  "([0]+[1]*x) * sin(2*pi*x/(2*pi/12/8/16)) + ([1]+[2]*x) * cos(2*pi*x/(2*pi/12/8/16)) + [3] + [4]*x + [5]*x*x + ([6]+[7]*x) * sin(4*pi*x/(2*pi/12/8/16)) + ([8]+[9]*x) * cos(4*pi*x/(2*pi/12/8/16))+ ([10] + [11]*x) * sin(3*pi*x/(2*pi/12/8/16)) + ([12]+[13]*x) * cos(3*pi*x/(2*pi/12/8/16))+ ([14]+[15]*x) * sin(pi*x/(5*pi/12/8/16)) + ([16]+[17]*x) * cos(pi*x/(5*pi/12/8/16)) + ([18]+[19]*x) * x*x*x", -3, 0);
170  gcenter->Fit(fPeriod, "M");
171 
172  p = (TPad *) c1->cd(idx++);
173  c1->Update();
174 
175  TH2 *hresidual_phi_cor = (TH2 *) hresidual_phi->Clone("hresidual_phi_cor");
176  hresidual_phi_cor->Reset();
177  hresidual_phi_cor->SetTitle(";Global Phi [rad];Modulation-Corrected Phi Residual [cm]");
178 
179  double sum = 0;
180  for (int binx = 1; binx <= hresidual_phi_highres->GetNbinsX(); ++binx)
181  for (int biny = 1; biny <= hresidual_phi_highres->GetNbinsY(); ++biny)
182  {
183  const double x = hresidual_phi_highres->GetXaxis()->GetBinCenter(binx);
184  const double y = hresidual_phi_highres->GetYaxis()->GetBinCenter(biny);
185  const double value = hresidual_phi_highres->GetBinContent(binx, biny);
186 
187  sum += value;
188  // const double y_cor = y;
189  const double y_cor = y - fPeriod->Eval(x);
190 
191  hresidual_phi_cor->Fill(x, y_cor, value);
192  }
193  cout << __PRETTY_FUNCTION__ << " sum = " << sum << endl;
194 
195  hresidual_phi_cor->Draw("colz");
196 
197  p = (TPad *) c1->cd(idx++);
198  c1->Update();
199  TH1 *hresidual_cor = hresidual_phi_cor->ProjectionY("hresidual_cor", 1, hresidual_phi_cor->GetNbinsX());
200 
201  hresidual_cor->Draw();
202 
203  TH1 *hresidual = hresidual_phi->ProjectionY("hresidual", 1, hresidual_phi->GetNbinsX());
204  // hresidual->SetMarkerColor(kRed - 3);
205  hresidual->SetLineColor(kRed - 3);
206  hresidual->Draw("same");
207 
208  TF1 *hresidual_fgaus = new TF1("hresidual_fgaus", "gaus", -.2, .2);
209  //
210  TF1 *hresidual_dgaus = new TF1(Form("hresidual_dgaus"), "gaus + [3]*exp(-0.5*((x-[1])/[4])**2) + [5]",
211  -.2, .2);
212  //
213  hresidual_fgaus->SetParameter(1, hresidual->GetMean());
214  hresidual_fgaus->SetParameter(2, hresidual->GetRMS());
215  //
216  hresidual_cor->Fit(hresidual_fgaus, "MQ");
217  //
218  // TF1 f2(Form("dgaus"), "gaus + [3]",
219  // -fgaus.GetParameter(2) * 1.5, fgaus.GetParameter(2) * 1.5);
220  //
221  hresidual_dgaus->SetParameters(hresidual_fgaus->GetParameter(0),
222  hresidual_fgaus->GetParameter(1),
223  hresidual_fgaus->GetParameter(2) / 2,
224  hresidual_fgaus->GetParameter(0) / 20,
225  hresidual_fgaus->GetParameter(2) * 10, 0);
226  //
227  hresidual_cor->Fit(hresidual_dgaus, "M");
228  hresidual_dgaus->SetLineColor(kBlue + 2);
229  const double resolution = TMath::Min(hresidual_dgaus->GetParameter(2), hresidual_dgaus->GetParameter(4));
230 
231  cout << __PRETTY_FUNCTION__ << " hresidual_fgaus->GetParameter(2) = " << hresidual_dgaus->GetParameter(2)
232  << " hresidual_fgaus->GetParameter(4) = " << hresidual_dgaus->GetParameter(4) << " resolution = " << resolution << endl;
233  //
234  // // new TCanvas;
235  // // h1->Draw();
236  // // fgaus.Draw("same");
237  // // break;
238  //
239  // x[n] = p2->GetBinCenter(i);
240  // ex[n] = (p2->GetBinCenter(2) - p2->GetBinCenter(1)) / 2;
241  // y[n] = fgaus.GetParameter(1);
242  // ey[n] = fgaus.GetParError(1);
243 
244  gPad->SetTopMargin(.3);
245  TLegend *leg = new TLegend(.03, .7, .98, .99, description + ": 1-removed residual");
246  leg->AddEntry(hresidual, "Raw residual", "l");
247  leg->AddEntry(hresidual_cor, "Position dependence corected residual", "p");
248  leg->AddEntry(hresidual_dgaus, Form("resolution = %.0f #mum", resolution * 1e4), "l");
249  leg->Draw();
250 
251  SaveCanvas(c1,
252  TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
253 }
254 
255 void TrackQA()
256 {
257  TCanvas *c1 = new TCanvas("TrackQA", "TrackQA", 1800, 860);
258  c1->Divide(3, 2);
259  int idx = 1;
260  TPad *p = nullptr;
261 
262  p = (TPad *) c1->cd(idx++);
263  c1->Update();
264 
265  TH1 *hnTrack = new TH1F("hnTrack", ";# track per event", 17, -.5, 16.5);
266  T->Draw("nTrack>>hnTrack");
267 
268  p = (TPad *) c1->cd(idx++);
269  c1->Update();
270 
271  TH1 *hCluster = new TH1F("hCluster", ";# cluster per track", 16, .5, 16.5);
272  T->Draw("TPCTrack.nCluster>>hCluster");
273 
274  p = (TPad *) c1->cd(idx++);
275  c1->Update();
276 
277  TH1 *hclusterSizePhi = new TH1F("hclusterSizePhi", ";Cluster phi size", 16, .5, 16.5);
278  T->Draw("TPCTrack.clusterSizePhi>>hclusterSizePhi");
279 
280  p = (TPad *) c1->cd(idx++);
281  c1->Update();
282 
283  TH1 *hclusterProjectionPhi = new TH1F("hclusterProjectionPhi", ";Cluster phi [rad]", 100, -3.5, -2.5);
284  T->Draw("TPCTrack.clusterProjectionPhi>>hclusterProjectionPhi");
285 
286  p = (TPad *) c1->cd(idx++);
287  c1->Update();
288 
289  TH1 *hAngle = new TH1F("hAngle", ";Horizontal angle [degree]", 100, -30, 30);
290  T->Draw("atan(TPCTrack.pz/ TPCTrack.px)/pi*180>>hAngle");
291 
292  p = (TPad *) c1->cd(idx++);
293  c1->Update();
294 
295  TH1 *hAngleV = new TH1F("hAngleV", ";Vertical angle [degree]", 100, -30, 30);
296  T->Draw("(atan(TPCTrack.py/ TPCTrack.px) - 2*pi/12/2 )/pi*180>>hAngleV");
297  // TH1 *hresidualRough = new TH1F("hresidualRough", ";Rought phi residual [cm]", 1000, -1, 1);
298  // T->Draw("TPCTrack.clusterResidualPhi>>hresidualRough", "Iteration$ >= 5 && Iteration$ <= 10 && TPCTrack.nCluster>=10");
299 
300  SaveCanvas(c1,
301  TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
302 }
303 
304 void Track3D()
305 {
306  TCanvas *c1 = new TCanvas("Track3D", "Track3D", 1800, 900);
307  c1->Divide(2, 1);
308  int idx = 1;
309  TPad *p = nullptr;
310 
311  p = (TPad *) c1->cd(idx++);
312  c1->Update();
313  // p->SetLogx();
314  // p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
315  //
316  TH3F *h3ClusterOverlay = new TH3F("h3ClusterOverlay", "h3ClusterOverlay", 128, -40, 10, 64, 35, 65, 128, -15, 15);
317  //
318  T->SetAlias("PhiCenter", "pi/12 + pi"); // center line azimuthal angle for TPC sector 0
319  T->Draw("TPCTrack.clusterX*cos(PhiCenter + pi/2) + TPCTrack.clusterY*sin(PhiCenter + pi/2):TPCTrack.clusterX*cos(PhiCenter) + TPCTrack.clusterY*sin(PhiCenter):TPCTrack.clusterZ>>h3ClusterOverlay", "", "BOX2");
320  // "TPCTrack.clusterE", "BOX2");
321  h3ClusterOverlay->SetTitle(";Drift Direction [cm];Pad Row Direction [cm];Azimuth Direction [cm]");
322  h3ClusterOverlay->SetLineWidth(0);
323  //
324  p = (TPad *) c1->cd(idx++);
325  c1->Update();
326  // p->SetLogx();
327  // p->DrawFrame(0, 0, 10, 2, ";Transverse Momentum, p_{T} [GeV/c];Nuclear Modification Factor, R_{AA}");
328 
329  TH1 *h2ClusterOverlay = h3ClusterOverlay->Project3D("zx");
330  h2ClusterOverlay->Draw("colz");
331  // eventT->Draw("Clusters.avg_pady:Clusters.min_sample+Clusters.peak_sample>>h2ClusterOverlay",
332  // "Clusters.peak", "colz");
333  // h2ClusterOverlay->SetTitle(";Time [0-127*50ns];Pads [0-127]");
334  // h2ClusterOverlay->SetLineWidth(0);
335 
336  p->SetTopMargin(.9);
337  TLegend *leg = new TLegend(.15, .9, .95, .99, description + ": accumulated clusters on tracks");
338  leg->Draw();
339 
340  SaveCanvas(c1,
341  TString(_file0->GetName()) + TString(c1->GetName()), false);
342 }
343 
344 void TrackDistortion(const TCut &cut = "TPCTrack.nCluster>=10")
345 {
346  TCanvas *c1 = new TCanvas("TrackDistortion", "TrackDistortion", 1800, 860);
347  c1->Divide(2, 1);
348  int idx = 1;
349  TPad *p = nullptr;
350 
351  p = (TPad *) c1->cd(idx++);
352  c1->Update();
353 
354  TH2 *hPhiDistortion = new TH2F("hPhiDistortion", ";Pad Layers;Phi Residual [cm]", 16, -.5, 15.5, 200, -3, 3);
355  T->Draw("TPCTrack.clusterResidualPhi:Iteration$>>hPhiDistortion", cut, "goff");
356  TGraph *gPhiDistortion = FitProfile(hPhiDistortion);
357  hPhiDistortion->Draw("colz");
358  gPhiDistortion->Draw("p");
359 
360  p = (TPad *) c1->cd(idx++);
361  c1->Update();
362 
363  TH2 *hZDistortion = new TH2F("hZDistortion", ";Pad Layers;Z Residual [cm]", 16, -.5, 15.5, 200, -3, 3);
364  T->Draw("TPCTrack.clusterResidualZ:Iteration$>>hZDistortion", cut, "goff");
365  TGraph *gZDistortion = FitProfile(hZDistortion);
366  hZDistortion->Draw("colz");
367  gZDistortion->Draw("p");
368 
369  SaveCanvas(c1,
370  TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
371 }
372 
373 Double_t langaufun(Double_t *x, Double_t *par)
374 {
375  //Fit parameters:
376  //par[0]=Width (scale) parameter of Landau density
377  //par[1]=Most Probable (MP, location) parameter of Landau density
378  //par[2]=Total area (integral -inf to inf, normalization constant)
379  //par[3]=Width (sigma) of convoluted Gaussian function
380  //
381  //In the Landau distribution (represented by the CERNLIB approximation),
382  //the maximum is located at x=-0.22278298 with the location parameter=0.
383  //This shift is corrected within this function, so that the actual
384  //maximum is identical to the MP parameter.
385 
386  // Numeric constants
387  Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
388  Double_t mpshift = -0.22278298; // Landau maximum location
389 
390  // Control constants
391  Double_t np = 100.0; // number of convolution steps
392  Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
393 
394  // Variables
395  Double_t xx;
396  Double_t mpc;
397  Double_t fland;
398  Double_t sum = 0.0;
399  Double_t xlow, xupp;
400  Double_t step;
401  Double_t i;
402 
403  // MP shift correction
404  mpc = par[1] - mpshift * par[0];
405 
406  // Range of convolution integral
407  xlow = x[0] - sc * par[3];
408  xupp = x[0] + sc * par[3];
409 
410  step = (xupp - xlow) / np;
411 
412  // Convolution integral of Landau and Gaussian by sum
413  for (i = 1.0; i <= np / 2; i++)
414  {
415  xx = xlow + (i - .5) * step;
416  fland = TMath::Landau(xx, mpc, par[0]) / par[0];
417  sum += fland * TMath::Gaus(x[0], xx, par[3]);
418 
419  xx = xupp - (i - .5) * step;
420  fland = TMath::Landau(xx, mpc, par[0]) / par[0];
421  sum += fland * TMath::Gaus(x[0], xx, par[3]);
422  }
423 
424  return (par[2] * step * sum * invsq2pi / par[3]);
425 }
426 
427 TF1 *langaufit(TH1 *his,
428  const vector<double> &fitrange,
429  const vector<double> &startvalues,
430  const vector<double> &parlimitslo, const vector<double> &parlimitshi
431  // ,
432  // const vector<double> &fitparams, const vector<double> &fiterrors,
433  // const vector<double> &ChiSqr, Int_t *NDF
434 )
435 {
436  // Once again, here are the Landau * Gaussian parameters:
437  // par[0]=Width (scale) parameter of Landau density
438  // par[1]=Most Probable (MP, location) parameter of Landau density
439  // par[2]=Total area (integral -inf to inf, normalization constant)
440  // par[3]=Width (sigma) of convoluted Gaussian function
441  //
442  // Variables for langaufit call:
443  // his histogram to fit
444  // fitrange[2] lo and hi boundaries of fit range
445  // startvalues[4] reasonable start values for the fit
446  // parlimitslo[4] lower parameter limits
447  // parlimitshi[4] upper parameter limits
448  // fitparams[4] returns the final fit parameters
449  // fiterrors[4] returns the final fit errors
450  // ChiSqr returns the chi square
451  // NDF returns ndf
452 
453  Int_t i;
454  Char_t FunName[100];
455 
456  sprintf(FunName, "Fitfcn_%s", his->GetName());
457 
458  TF1 *ffitold = (TF1 *) gROOT->GetListOfFunctions()->FindObject(FunName);
459  if (ffitold) delete ffitold;
460 
461  TF1 *ffit = new TF1(FunName, langaufun, fitrange[0], fitrange[1], 4);
462  ffit->SetParameters(startvalues.data());
463  ffit->SetParNames("Width", "MP", "Area", "GSigma");
464 
465  for (i = 0; i < 4; i++)
466  {
467  ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
468  }
469 
470  his->Fit(FunName, "RB0"); // fit within specified range, use ParLimits, do not plot
471 
472  // ffit->GetParameters(fitparams); // obtain fit parameters
473  // for (i=0; i<4; i++) {
474  // fiterrors[i] = ffit->GetParError(i); // obtain fit parameter errors
475  // }
476  // ChiSqr[0] = ffit->GetChisquare(); // obtain chi^2
477  // NDF[0] = ffit->GetNDF(); // obtain ndf
478 
479  return (ffit); // return fit function
480 }
481 
482 void TrackClusterEnergy(const TCut &cut = "TPCTrack.nCluster>=10")
483 {
484  TCanvas *c1 = new TCanvas("TrackClusterEnergy", "TrackClusterEnergy", 1200, 860);
485  // c1->Divide(2, 1);
486  int idx = 1;
487  TPad *p = nullptr;
488 
489  p = (TPad *) c1->cd(idx++);
490  c1->Update();
491  p->SetLogy();
492 
493  TH1 *hClusterEnergy = new TH1F("hClusterEnergy", ";Cluster Energy on good track [ADU];Count / bin", 2500, 0, 5000);
494  T->Draw("TPCTrack.clusterE>>hClusterEnergy", cut, "goff");
495  // TGraph *gPhiDistortion = FitProfile(hPhiDistortion);
496 
497  TF1 *fit = langaufit(hClusterEnergy,
498  {0, 5000},
499  {200, 500, (double) T->GetEntries(cut) * 10, 200},
500  {0, 0, 0, 0},
501  {5000, 5000, 1e10, 5000});
502  fit->SetLineColor(kBlue + 2);
503 
504  hClusterEnergy->Draw();
505  fit->Draw("same");
506  // gPhiDistortion->Draw("p");
507 
508  TLegend *leg = new TLegend(.3, .7, .95, .9, description + ": Cluster Energy on good track");
509  leg->AddEntry(hClusterEnergy, TString("Data: ") + cut.GetTitle(), "l");
510  leg->AddEntry(fit,
511  Form("Langdau * Gauss Fit, #mu= %.0f ADU", fit->GetParameter(1)), "l");
512  leg->Draw();
513 
514  SaveCanvas(c1,
515  TString(_file0->GetName()) + TString(c1->GetName()), kFALSE);
516 }
517 
518 void DrawTpcPrototypeGenFitTrkFitter(
519  const TString infile = "data/tpc_beam/tpc_beam_00000292-0000.evt_TpcPrototypeGenFitTrkFitter.root", const TString desc = "Run 292" //
520 )
521 {
522  gSystem->Load("libtpc2019.so");
523  //
524  SetsPhenixStyle();
525  TVirtualFitter::SetDefaultFitter("Minuit2");
526  gStyle->SetLegendTextSize(0);
527  //
528  if (!_file0)
529  {
530  TString chian_str = infile;
531  chian_str.ReplaceAll("ALL", "*");
532 
533  TChain *t = new TChain("T");
534  const int n = t->Add(chian_str);
535 
536  cout << "Loaded " << n << " root files with eventT in " << chian_str << endl;
537  assert(n > 0);
538 
539  T = t;
540 
541  _file0 = new TFile;
542  _file0->SetName(infile);
543  }
544  //
545  description = desc;
546 
547  TrackQA();
548  TrackDistortion();
549  Track3D();
550  TrackClusterEnergy();
551  // Resolution();
552  Resolution("Iteration$ >= 5 && Iteration$ <= 10 && TPCTrack.nCluster>=12 && TPCTrack.clusterSizePhi > 7");
553 }