Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DrawTpcPrototypeGenFitTrkFitter_Summary.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DrawTpcPrototypeGenFitTrkFitter_Summary.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 string basePath;
134 const static double x0 = 6 - (2 + 7. / 16); // inch
135 
136 pair<TGraph *, TGraph *> getResolutions(vector<pair<double, string>> datafiles, Color_t color = kRed + 1)
137 {
138  vector<double> x_cm;
139  vector<double> x;
140  vector<double> res;
141  vector<double> res_err;
142  vector<double> res2;
143  vector<double> res2_err;
144 
145  for (const auto &data : datafiles)
146  {
147  cout << __PRETTY_FUNCTION__ << " "
148  << " processing x = " << data.first << " from " << data.second << endl;
149 
150  TFile *f = TFile::Open((TString(basePath + data.second.c_str()) + "_TpcPrototypeGenFitTrkFitter.root_DrawJet_Resolution.root"));
151  assert(f);
152  assert(f->IsOpen());
153  TH1 *hresidual_cor = (TH1 *) f->GetObjectChecked("hresidual_cor", "TH1");
154  assert(hresidual_cor);
155  TF1 *hresidual_dgaus = (TF1 *) hresidual_cor->GetListOfFunctions()->At(0);
156  assert(hresidual_dgaus);
157 
158  if (hresidual_dgaus->GetParameter(2) < hresidual_dgaus->GetParameter(4))
159  {
160  res.push_back(hresidual_dgaus->GetParameter(2) * 1e4);
161  res_err.push_back(hresidual_dgaus->GetParError(2) * 1e4);
162  }
163  else
164  {
165  res.push_back(hresidual_dgaus->GetParameter(4) * 1e4);
166  res_err.push_back(hresidual_dgaus->GetParError(4) * 1e4);
167  }
168  x.push_back(data.first);
169  x_cm.push_back((data.first - x0) * 2.54);
170  res2.push_back(res[res.size() - 1] * res[res.size() - 1]);
171  res2_err.push_back(2 * res[res.size() - 1] * res_err[res.size() - 1]);
172 
173  cout << "x = " << x[res.size() - 1] << " x_cm = " << x_cm[res.size() - 1] << " res = " << res[res.size() - 1] << " res_err = " << res_err[res.size() - 1] << endl;
174  }
175 
176  TGraphErrors *g = new TGraphErrors(x.size(), x_cm.data(), res.data(), 0, res_err.data());
177  g->SetMarkerColor(color);
178  g->SetLineColor(color);
179  g->Print();
180 
181  TGraphErrors *g2 = new TGraphErrors(x.size(), x.data(), res2.data(), 0, res2_err.data());
182  g2->SetMarkerColor(color);
183  g2->SetLineColor(color);
184  g2->Print();
185 
186  return pair<TGraph *, TGraph *>(g, g2);
187 }
188 
189 void DrawTpcPrototypeGenFitTrkFitter_Summary( //
190  string basepath = "/phenix/u/jinhuang/links/sPHENIX_work/TPC/fnal_June2019/SimPadPlaneIter5/")
191 {
192  // gSystem->Load("libtpc2019.so");
193  //
194  basePath = basepath;
195  const double max_res = 300;
196  const double field_off_trans_diffusion = 150; // um/sqrt(cm), From: Prakhar Garg <prakhar.garg@stonybrook.edu> Date: Tue, Jun 18, 2019 at 1:04 PM
197  const double field_on_trans_diffusion = 40; // um/sqrt(cm), From: Prakhar Garg <prakhar.garg@stonybrook.edu> Date: Tue, Jun 18, 2019 at 1:04 PM
198 
199  SetsPhenixStyle();
200  TVirtualFitter::SetDefaultFitter("Minuit2");
201  gStyle->SetLegendTextSize(0);
202 
203  TCanvas *c1 = new TCanvas("DrawTpcPrototypeGenFitTrkFitter_Summary", "DrawTpcPrototypeGenFitTrkFitter_Summary", 1800, 700);
204  c1->Divide(2, 1);
205  int idx = 1;
206  TPad *p = nullptr;
207 
208  p = (TPad *) c1->cd(idx++);
209  c1->Update();
210 
211  TH1 *frame = p->DrawFrame(2, 0, 20, max_res * max_res);
212  frame->SetTitle(";Horizontal Position [in];Resolution^{2}, Cluster w/ 2+ pad [#mum^{2}]");
213  frame->GetYaxis()->SetTitleOffset(1.75);
214  TLine *l = new TLine(x0, 0, x0, max_res * max_res);
215  l->Draw();
216 
217  TLegend *leg1 = new TLegend(.35, .7, .9, .9, "2019 TPC Beam Test Preview");
218 
219  p = (TPad *) c1->cd(idx++);
220  c1->Update();
221  p->DrawFrame(0, 50, 125, max_res)->SetTitle(";Drift Length, L [cm];Resolution, Cluster w/ 2+ pad [#mum]");
222  TLegend *leg2 = new TLegend(.42, .45, .92, .93, "2019 TPC Beam Test Preview");
223 
224  {
225  Color_t color = kRed + 2;
226  const TString name = "Scan2: Gap 300V, GEM 370V";
227 
228  pair<TGraph *, TGraph *> scan2 = getResolutions(
229  {{6, "tpc_beam_00000288-0000.evt"},
230  {6, "tpc_beam_00000292-0000.evt"},
231  {10, "tpc_beam_00000293-0000.evt"},
232  {14, "tpc_beam_00000294-0000.evt"},
233  {18, "tpc_beam_00000295-0000.evt"}},
234  color);
235 
236  c1->cd(1);
237  scan2.second->Draw("p");
238  TF1 *fline_scan2 = new TF1("fline_scan2", "pol1", x0, 20);
239  scan2.second->Fit(fline_scan2, "MR0");
240  fline_scan2->SetLineColor(color);
241  fline_scan2->Draw("same");
242 
243  leg1->AddEntry(scan2.second, name, "p");
244 
245  c1->cd(2);
246  scan2.first->Draw("p");
247  TF1 *fdiff_scan2 = new TF1("fline_scan2", "sqrt([0]*[0] + ([1]*[1]/[2])*x)", 0, 40);
248  fdiff_scan2->SetParameters(70, field_off_trans_diffusion, 20);
249  fdiff_scan2->FixParameter(1, field_off_trans_diffusion);
250  scan2.first->Fit(fdiff_scan2, "MR0");
251  fdiff_scan2->SetLineColor(color);
252  fdiff_scan2->Draw("same");
253 
254  TF1 *fdiff_scan2_FieldON = (TF1 *) fdiff_scan2->Clone("fdiff_scan2_FieldON");
255  fdiff_scan2_FieldON->SetRange(0, 105);
256  fdiff_scan2_FieldON->SetParameter(1, field_on_trans_diffusion);
257  fdiff_scan2_FieldON->SetLineStyle(kDashed);
258  fdiff_scan2_FieldON->Draw("same");
259  // fdiff_scan2_FieldON->Print();
260 
261  leg2->AddEntry(scan2.first, name, "pe");
262  leg2->AddEntry(fdiff_scan2, Form("#sqrt{(%.1f #mum)^{2} + (%.0f #mum/#sqrt{cm} #times#sqrt{L/%.1f})^{2}}", abs(fdiff_scan2->GetParameter(0)), abs(fdiff_scan2->GetParameter(1)), abs(fdiff_scan2->GetParameter(2))), "l");
263  leg2->AddEntry(fdiff_scan2_FieldON, Form("Field = 1.4T, #sigma_{T, Diffusion} = %.0f #mum/#sqrt{cm}", field_on_trans_diffusion), "l");
264 
265  c1->Update();
266  }
267 
268  {
269  Color_t color = kBlue + 2;
270  const TString name = "Scan3: Gap 150V, GEM 380V";
271 
272  pair<TGraph *, TGraph *> scan3 = getResolutions(
273  {{6, "tpc_beam_00000297-0000.evt"},
274  {6, "tpc_beam_00000298-0000.evt"},
275  {10, "tpc_beam_00000299-0000.evt"},
276  {14, "tpc_beam_00000300-0000.evt"},
277  {18, "tpc_beam_00000301-0000.evt"}},
278  color);
279 
280  c1->cd(1);
281  scan3.second->Draw("p");
282  TF1 *fline_scan3 = new TF1("fline_scan3", "pol1", x0, 20);
283  scan3.second->Fit(fline_scan3, "MR0");
284  fline_scan3->SetLineColor(color);
285  fline_scan3->Draw("same");
286  leg1->AddEntry(scan3.second, name, "p");
287 
288  c1->cd(2);
289  scan3.first->Draw("p");
290  TF1 *fdiff_scan3 = new TF1("fline_scan3", "sqrt([0]*[0] + ([1]*[1]/[2])*x)", 0, 40);
291  fdiff_scan3->SetParameters(70, field_off_trans_diffusion, 20);
292  fdiff_scan3->FixParameter(1, field_off_trans_diffusion);
293  scan3.first->Fit(fdiff_scan3, "MR0");
294  fdiff_scan3->SetLineColor(color);
295  fdiff_scan3->Draw("same");
296 
297  TF1 *fdiff_scan3_FieldON = (TF1 *) fdiff_scan3->Clone("fdiff_scan3_FieldON");
298  fdiff_scan3_FieldON->SetRange(0, 105);
299  fdiff_scan3_FieldON->SetParameter(1, field_on_trans_diffusion);
300  fdiff_scan3_FieldON->SetLineStyle(kDashed);
301  fdiff_scan3_FieldON->Draw("same");
302  // fdiff_scan3_FieldON->Print();
303 
304  leg2->AddEntry(scan3.first, name, "pe");
305  leg2->AddEntry(fdiff_scan3, Form("#sqrt{(%.1f #mum)^{2} + (%.0f #mum/#sqrt{cm} #times#sqrt{L/%.1f})^{2}}", abs(fdiff_scan3->GetParameter(0)), abs(fdiff_scan3->GetParameter(1)), abs(fdiff_scan3->GetParameter(2))), "l");
306  leg2->AddEntry(fdiff_scan3_FieldON, Form("Field = 1.4T, #sigma_{T, Diffusion} = %.0f #mum/#sqrt{cm}", field_on_trans_diffusion), "l");
307 
308  c1->Update();
309  }
310 
311  if (0)
312  {
313  Color_t color = kGreen + 2;
314 
315  pair<TGraph *, TGraph *> scan4 = getResolutions(
316  {{6, "tpc_beam_00000357-0000.evt"},
317  {6, "tpc_beam_00000358-0000.evt"},
318  {10, "tpc_beam_00000359-0000.evt"},
319  {10, "tpc_beam_00000360-0000.evt"},
320  {13, "tpc_beam_00000356-0000.evt"},
321  {14, "tpc_beam_00000361-0000.evt"},
322  {14, "tpc_beam_00000362-0000.evt"},
323  {18, "tpc_beam_00000363-0000.evt"},
324  {18, "tpc_beam_00000366-0000.evt"}},
325  color);
326 
327  c1->cd(1);
328  scan4.second->Draw("p");
329  TF1 *fline_scan4 = new TF1("fline_scan4", "pol1", x0, 16);
330  scan4.second->Fit(fline_scan4, "MR0");
331  fline_scan4->SetLineColor(color);
332  fline_scan4->Draw("same");
333 
334  c1->cd(2);
335  scan4.first->Draw("p");
336  TF1 *fdiff_scan4 = new TF1("fline_scan4", "sqrt([0]*[0] + [1]*[1]*x)", 0, 30);
337  scan4.first->Fit(fdiff_scan4, "MR0");
338  fdiff_scan4->SetLineColor(color);
339  fdiff_scan4->Draw("same");
340 
341  c1->Update();
342  }
343 
344  c1->cd(1);
345  leg1->Draw();
346  c1->cd(2);
347  leg2->Draw();
348  SaveCanvas(c1,
349  TString(basePath) + TString(c1->GetName()), true);
350 }