Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FitVerticalScan.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FitVerticalScan.C
1 // $Id: $
2 
11 #include <cmath>
12 #include <vector>
13 #include <TFile.h>
14 #include <TString.h>
15 #include <TLine.h>
16 #include <TGraphErrors.h>
17 #include <TLatex.h>
18 #include <TGraphErrors.h>
19 #include <TMultiGraph.h>
20 #include <TLegend.h>
21 #include <TF1.h>
22 #include <cassert>
23 #include <iostream>
24 #include <fstream>
25 
26 #include "SaveCanvas.C"
27 #include "SetOKStyle.C"
28 using namespace std;
29 
30 vector<TGraphErrors *> vg_src;
31 
32 void
33 OnlineAnalysis(TString filename =
34  "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/column3_214x.txt")
35 {
36  std::map<int, float> runtoz;
37 
38  runtoz[2138] = 29.7;
39  runtoz[2140] = 55.2;
40  runtoz[2141] = 80.6;
41  runtoz[2142] = 106.3;
42  runtoz[2143] = 131.5;
43  runtoz[2144] = 156.6;
44  runtoz[2146] = 182.4;
45 
46  runtoz[2171] = 142.1;
47  runtoz[2172] = 122.0;
48  runtoz[2174] = 102.1;
49  runtoz[2176] = 81.9;
50  runtoz[2177] = 61.9;
51  runtoz[2178] = 42.1;
52  runtoz[2179] = 21.9;
53 
54  runtoz[2181] = 21.85;
55  runtoz[2182] = 41.99;
56  runtoz[2183] = 61.87;
57  runtoz[2184] = 61.87;
58  runtoz[2185] = 82.12;
59  runtoz[2186] = 101.9;
60  runtoz[2187] = 121.9;
61  runtoz[2188] = 142.0;
62  runtoz[2189] = 152.0;
63  runtoz[2190] = 162.1;
64  runtoz[2191] = 172.0;
65  runtoz[2192] = 172.0;
66 
67  TCanvas *c1 = new TCanvas("c1", "Vertical scan", 640, 480);
68 
69  TLegend *leg = new TLegend(0.9122257, 0.6483516, 0.9968652, 0.9010989, NULL,
70  "brNDC");
71  leg->SetFillColor(0);
72 
73  int rn;
74  int tower;
75  float mpv;
76  float dmpv;
77 
78  TGraphErrors *gr[8];
79 
80  vector<float> v_x[8];
81  vector<float> v_y[8];
82  vector<float> v_dx[8];
83  vector<float> v_dy[8];
84 
85  ifstream sumfile;
86  sumfile.open(filename);
87 
88  std::string line;
89  int linenum = 0;
90  int index;
91  while (sumfile >> rn >> tower >> mpv >> dmpv)
92  {
93  index = (tower / 8);
94  cout << index << " " << rn << " " << runtoz[rn] << " " << mpv << " "
95  << dmpv << endl;
96  v_x[index].push_back(runtoz[rn]);
97  v_dx[index].push_back(0.2);
98  v_y[index].push_back(mpv);
99  v_dy[index].push_back(dmpv);
100  linenum++;
101  }
102 
103  TMultiGraph *mg = new TMultiGraph();
104 
105  TString towername;
106  int col = 1;
107  int marker = 21;
108  for (int i = 0; i < 8; i++)
109  {
110  gr[i] = new TGraphErrors(v_x[i].size(), &v_x[i][0], &v_y[i][0],
111  &v_dx[i][0], &v_dy[i][0]);
112  gr[i]->SetMarkerStyle(marker);
113  gr[i]->SetMarkerSize(0.4);
114  gr[i]->SetLineColor(col);
115  gr[i]->SetMarkerColor(col);
116  gr[i]->Fit("expo");
117  gr[i]->GetFunction("expo")->SetLineColor(col);
118 
119  vg_src.push_back((TGraphErrors *) gr[i]->Clone());
120 
121  towername = "Tower ";
122  towername += i;
123  leg->AddEntry(gr[i], towername, "LP");
124  mg->Add(gr[i]);
125  col++;
126  marker++;
127  }
128 
129  mg->Draw("AP");
130  TString title = "Longitudinal EMCAL scan (";
131  title += filename;
132  title += ")";
133  mg->SetTitle(title);
134  mg->GetXaxis()->SetTitle("Vertical table position (mm)");
135  mg->GetYaxis()->SetTitle("MPV of 200 GeV p (ADC counts)");
136  leg->Draw();
137 
138  c1->Update();
139  c1->Modified();
140 
141 }
142 
143 void
144 FitVerticalScan(const TString input =
145 // "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/mpvPerTower.root"
146  "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/getMeanRMS_Summary.root")
147 {
148 // gROOT->SetStyle("Modern");
149  SetOKStyle();
150 
151  TFile * _file0 = new TFile(input);
152  assert(_file0->IsOpen());
153 
154 // TCanvas * c2 = (TCanvas *) _file0->Get("c2");
155  TCanvas * c2 = (TCanvas *) _file0->Get("getMeanRMS_Summary");
156  assert(c2);
157 
158  c2->Draw();
159 
160  int i = 2;
161  vg_src.push_back(
162  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
163  vg_src.push_back(
164  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
165  vg_src.push_back(
166  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
167  vg_src.push_back(
168  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
169  vg_src.push_back(
170  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
171  vg_src.push_back(
172  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
173  vg_src.push_back(
174  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
175  vg_src.push_back(
176  (TGraphErrors *) c2->GetListOfPrimitives()->At(i++)->Clone());
177 
178 // OnlineAnalysis();
179 
180  TGraphErrors * g_final = new TGraphErrors((vg_src[0])->GetN());
181  g_final->SetMarkerSize(3);
182  g_final->SetMarkerStyle(kFullCircle);
183  g_final->SetMarkerColor(kBlue + 3);
184  g_final->SetFillColor(kWhite);
185  g_final->SetLineWidth(3);
186  g_final->SetFillStyle(0);
187 
188 // TF1 * f_att = new TF1("f_att","(exp(-((14-x)/[0])) + 0.3 * exp(-((x-10 + 13)/[0])) )*[1]",0,20);
189 // f_att->SetParameters(105,1);
190 // f_att->SetLineColor(kBlue+3);
191 // f_att->SetLineWidth(4);
192 // TF1 * f_att = new TF1("f_att","(exp(-((14-x)/[0])) + 0.3 * exp(-((x + 13)/[0])) )*[1]",0,20);
193 // f_att->SetParameters(105,1);
194  TF1 * f_att = new TF1("f_att", "(exp(-((13.9-x)/[0])))*[1]", 0, 20);
195  f_att->SetParameters(105, 1);
196 
197  f_att->SetLineColor(kBlue + 3);
198  f_att->SetLineWidth(4);
199 
200  const int N_rows = 8;
201 
202  TGraphErrors * g_src = vg_src[0];
203 // g_src->Print();
204  for (int j = 0; j < (g_src->GetN()); ++j)
205  {
206  assert(g_src);
207 
208  (g_final->GetY())[j] = 0;
209  (g_final->GetEY())[j] = 0;
210  (g_final->GetEX())[j] = 0;
211  (g_final->GetX())[j] = (g_src->GetX())[j] / 10 - 0.8;
212  }
213 
214  for (int i = 0; i < N_rows; ++i)
215  {
216  new TCanvas();
217 
218  TGraphErrors * g_src = vg_src[i];
219 
220  assert(g_src);
221 
222 // g_src->Set
223  g_src->DrawClone("ap*");
224 
225  double sum = 0;
226  for (int j = 0; j < (g_src->GetN()); ++j)
227  {
228  sum += (g_src->GetY())[j];
229  }
230  double average = sum / g_src->GetN();
231 
232  for (int j = 0; j < (g_src->GetN()); ++j)
233  {
234  (g_src->GetY())[j] /= average;
235  (g_src->GetEY())[j] /= average;
236 
237  (g_final->GetY())[j] += (g_src->GetY())[j];
238  (g_final->GetEY())[j] += (g_src->GetEY())[j] * (g_src->GetEY())[j];
239  }
240  }
241 
242  for (int j = 0; j < (g_src->GetN()); ++j)
243  {
244  (g_final->GetY())[j] /= N_rows;
245  (g_final->GetEY())[j] = sqrt((g_final->GetEY())[j]) / N_rows;
246  }
247 
248  TText * t;
249  TCanvas *c1 = new TCanvas("FitVerticalScan", "FitVerticalScan", 1800, 650);
250 
251  gPad->SetGridx(0);
252  gPad->SetGridy(0);
253  gPad->SetRightMargin(0.99);
254  gPad->SetLeftMargin(0.15);
255  gPad->SetBottomMargin(0.2);
256 
257  TH1 * frame = gPad->DrawFrame(0, .9, 13.9, 1.1);
258  frame->SetTitle(
259  ";Beam position along the length of module (cm);MIP amplitude (A. U.)");
260  frame->GetXaxis()->SetTitleOffset(1);
261  frame->GetXaxis()->SetTitleSize(0.08);
262  frame->GetXaxis()->SetLabelSize(0.08);
263  frame->GetYaxis()->SetTitleOffset(.7);
264  frame->GetYaxis()->SetTitleSize(0.08);
265  frame->GetYaxis()->SetLabelSize(0.08);
266  frame->GetXaxis()->CenterTitle();
267 
268 // f_att->DrawClone("same");
269  g_final->Draw("p");
270  g_final->Fit(f_att, "MR");
271 // g_final->Fit("pol1");
272  g_final->Print();
273 
274  TLegend *leg = new TLegend(.15, .7, .8, .9);
275  leg->AddEntry(g_final, "Data, Averaged over 8 towers","ep");
276  leg->AddEntry(f_att,
277  Form("Fit, C#timesExp[ -(13.9 - x) / L_{eff} ], L_{eff} = %.0f#pm%.0f cm",
278  f_att->GetParameter(0), f_att->GetParError(0)),"l");
279 
280  leg->Draw();
281 
282  TLatex latex;
283  latex.SetTextSize(0.08);
284  latex.SetTextAlign(12);
285  latex.DrawLatex(0,0.92,"#leftarrow Mirror side");
286  latex.SetTextAlign(32);
287  latex.DrawLatex(13.9,0.92,"Lightguide side #rightarrow");
288 
289  SaveCanvas(c1, TString(_file0->GetName()) + TString(c1->GetName()), kTRUE);
290 }
291