Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
getMeanRMS.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file getMeanRMS.C
1 #include <cmath>
2 #include <vector>
3 #include <TFile.h>
4 #include <TString.h>
5 #include <TLine.h>
6 #include <TGraph.h>
7 #include <TLatex.h>
8 #include <TAxis.h>
9 #include <TGraphErrors.h>
10 #include <TMultiGraph.h>
11 #include <TLegend.h>
12 #include <TCanvas.h>
13 #include <TF1.h>
14 #include <TH1.h>
15 #include <TMath.h>
16 #include <cassert>
17 #include <iostream>
18 #include <fstream>
19 
20 #include "SaveCanvas.C"
21 
22 Double_t
23 langaufun(Double_t *x, Double_t *par)
24 {
25 
26  //Fit parameters:
27  //par[0]=Width (scale) parameter of Landau density
28  //par[1]=Most Probable (MP, location) parameter of Landau density
29  //par[2]=Total area (integral -inf to inf, normalization constant)
30  //par[3]=Width (sigma) of convoluted Gaussian function
31  //
32  //In the Landau distribution (represented by the CERNLIB approximation),
33  //the maximum is located at x=-0.22278298 with the location parameter=0.
34  //This shift is corrected within this function, so that the actual
35  //maximum is identical to the MP parameter.
36 
37  // Numeric constants
38  Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
39  Double_t mpshift = -0.22278298; // Landau maximum location
40 
41  // Control constants
42  Double_t np = 100.0; // number of convolution steps
43  Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
44 
45  // Variables
46  Double_t xx;
47  Double_t mpc;
48  Double_t fland;
49  Double_t sum = 0.0;
50  Double_t xlow, xupp;
51  Double_t step;
52  Double_t i;
53 
54  // MP shift correction
55  mpc = par[1] - mpshift * par[0];
56 
57  // Range of convolution integral
58  xlow = x[0] - sc * par[3];
59  xupp = x[0] + sc * par[3];
60 
61  step = (xupp - xlow) / np;
62 
63  // Convolution integral of Landau and Gaussian by sum
64  for (i = 1.0; i <= np / 2; i++)
65  {
66  xx = xlow + (i - .5) * step;
67 
68  fland = TMath::Landau(xx, mpc, par[0]) / par[0];
69  sum += fland * TMath::Gaus(x[0], xx, par[3]);
70 
71  xx = xupp - (i - .5) * step;
72  fland = TMath::Landau(xx, mpc, par[0]) / par[0];
73  sum += fland * TMath::Gaus(x[0], xx, par[3]);
74  }
75 
76  return (par[2] * step * sum * invsq2pi / par[3]);
77 }
78 
79 void
80 getMeanRMS(void)
81 {
82  ofstream Output;
83  Output.open("meanRMS_v1.txt");
84 
85  char name[20000];
86 
87  int runList[100] =
88  { 2181, 2182, 2184, 2185, 2186, 2187, 2188, 2189, 2190, 2192 };
89  double vert[10] =
90  { 22, 42, 62, 82, 102, 122, 142, 152, 162, 172 };
91  double mu[8][100];
92  double rms[8][100];
93  double err[8][100];
94  double vt[8][100];
95  double up[8][100];
96  double mip[80] =
97  { 43.4997, 45.6349, 40.8449, 47.7902, 48.4663, 39.2068, 48.3755, 83.1102 };
98 
99  for (int j = 0; j < 8; j++)
100  {
101  sprintf(name,
102  "/phenix/u/jinhuang/links/sPHENIX_work/Prototype_2016/Production_1101_EMCalSet2_HCalPR12/beam_0000%d-0000_DSTReader.root_DrawPrototype2MIPAnalysis_EMCDistribution_RAW_Valid_HODO_MIP_Col2_PedestalOther_TimingGood.root",
103  runList[j]);
104  TFile *_file0 = TFile::Open(name);
105 
106  TCanvas *c1 = new TCanvas(Form("getMeanRMS_run%d", runList[j]),
107  Form("getMeanRMS_run%d", runList[j]), 1800, 950);
108  c1->Divide(4, 2);
109 
110  for (int i = 0; i < 8; i++)
111  {
112  double min = 0;
113  double max = 150;
114  //double max=400;
115  //if(runList[j]==2418 || runList[j]==2417)max=400;
116  //if(runList[j]==2551)max=200;
117  //if(runList[j]==2560)max=400;
118  /*TF1 *fitFunc = new TF1("fitFunc","[0]*TMath::Landau(x,[1],[2])+[3]*TMath::Gaus(x,[4],[5])",min,max);
119  //TF1 *fitFunc = new TF1("fitFunc","[0]*TMath::Landau(x,[1],[2])",min,max);
120  //TF1 *fitFunc = new TF1("fitFunc","[0]*TMath::Gaus(x,[1],[2])",min,max);
121  fitFunc->SetParameter(0,100);
122  fitFunc->SetParameter(1,10);
123  fitFunc->SetParameter(2,10);
124  fitFunc->SetParameter(3,10);
125  fitFunc->SetParameter(4,10);
126  fitFunc->SetParameter(5,10);*/
127 
128  c1->cd(i + 1);
129  c1->cd(i + 1)->SetLogy();
130  sprintf(name, "hEnergy_ieta2_iphi%d", i);
131  TH1F *h = (TH1F*) gDirectory->Get(name)->Clone();
132  assert(h);
133  h->Draw();
134 
135  TF1 *ffit = new TF1(TString(h->GetName()) + "_langaufun", langaufun,
136  0, 500, 4);
137  ffit->SetParameters(6, 50, 3000, 15);
138  ffit->SetParNames("Width", "MP", "Area", "GSigma");
139  ffit->SetParLimits(0,0,20);
140  ffit->SetParLimits(1,0,200);
141  ffit->SetParLimits(2,0,20000);
142 // ffit->SetParLimits(3,10,20);
143  ffit->SetParLimits(3,1,30);
144 
145  h->Fit(ffit, "RBM");
146 
147  c1->Update();
148 
149 // Output << col[j] << "\t" << i << "\t" << fitFunc->GetMaximumX(20, 500)
150 // << endl;
151  Output << runList[j] << "\t" << vert[j] << "\t" << i << "\t"
152  << ffit->GetParameter(1) << "\t" << ffit->GetParError(1) << "\t"
153  << h->GetEntries() << endl;
154 // mu[i][j]=h->GetMean();
155  mu[i][j] = ffit->GetParameter(1);
156 // h->Fit("landau");
157 // mu[i][j] = landau->GetParameter(1);
158 // rms[i][j] = h->GetRMS();
159  err[i][j] = ffit->GetParError(1);
160  vt[i][j] = vert[j]; //+i*.1;
161 // up[i][j] =
162 // (h->GetMean() + 1.281 * h->GetRMS() / sqrt(h->GetEntries()))
163 // / mip[i];
164 // //cout<<up[i][j]<<endl;
165 // Output << runList[j] << "\t" << vert[j] << "\t" << i << "\t"
166 // << h->GetMean() << "\t" << h->GetRMS() << "\t" << h->GetEntries()
167 // << "\t" << up[i][j] << endl;
168 
169  //if(vert[j]<160 || vert[j]>170)vt[i][j]=vert[j]+i*.4;
170 
171  }
172 // sprintf(name, "meanRMS_%d_v2.root", runList[j]);
173 // c1->Print(name);
174  SaveCanvas(c1, TString(name) + TString(c1->GetName()), kTRUE);
175  }
176 
177 // double avg[7];
178 // for (int j = 0; j < 7; j++) //veritcal
179 // {
180 // double sum = 0.;
181 // for (int i = 0; i < 8; i++) //tower
182 // {
183 // sum = sum + mu[i][j];
184 // }
185 // avg[j] = sum / 8.;
186 // }
187 
188 // TCanvas c4("c4");
189 // for (int i = 0; i < 8; i++)
190 // {
191 // TGraphErrors *gr = new TGraphErrors(7, vert, avg, 0, 0);
192 // gr->SetMarkerStyle(20);
193 // //gr->SetMaximum(27);
194 // //gr->SetMinimum(23);
195 // gr->SetTitle("");
196 // gr->GetXaxis()->SetTitle("vertical position (mm)");
197 // gr->GetYaxis()->SetTitle("<MPV ADC>");
198 // gr->Draw("ap");
199 // }
200 
201  TLegend * leg = new TLegend(0.91, 0.51, 0.99, 0.89, NULL, "brNDC");
202 
203  TCanvas * c1 = new TCanvas(Form("getMeanRMS_Summary"),
204  Form("getMeanRMS_Summary"), 1800, 950);
205  gPad->DrawFrame(0, 0, 200, 200);
206  for (int i = 0; i < 8; i++)
207  {
208  TGraphErrors *gr = new TGraphErrors(7, vt[i], mu[i], 0, err[i]);
209  //TGraphErrors *gr = new TGraphErrors(11,vt[i],mu[i],0,err[i]);
210  gr->SetMarkerStyle(20);
211  gr->SetMarkerColor(i + 1);
212  //gr->SetMaximum(32);
213  //gr->SetMinimum(-14);
214  gr->SetMaximum(60);
215  gr->SetMinimum(10);
216  gr->SetTitle("");
217  gr->GetXaxis()->SetTitle("vertical position (mm)");
218  gr->GetYaxis()->SetTitle("MPV ADC");
219 // if (i == 0)
220 // gr->Draw("ap");
221 // if (i > 0)
222  gr->Draw("p");
223  //gr->Fit("pol1");
224  //gr->GetXaxis()->TAxis::SetNdivisions(505);
225  //gr->GetXaxis()->SetLimits(160.6,166);
226  //gr->GetXaxis()->SetLimits(130,190);
227  sprintf(name, "row %d", i);
228  leg->AddEntry(gr, name, "p");
229 
230 // gr->Draw("p");
231  }
232  leg->Draw();
233  SaveCanvas(c1, TString(c1->GetName()), kTRUE);
234 
235 // leg2 = new TLegend(0.54, 0.35, 0.62, 0.73, NULL, "brNDC");
236 // TCanvas c3("c3");
237 // for (int i = 0; i < 8; i++)
238 // {
239 // TGraphErrors *gr2 = new TGraphErrors(11, vert, up[i], 0, 0);
240 // //TGraphErrors *gr2 = new TGraphErrors(11,vert,mu[i],0,err[i]);
241 // gr2->SetMarkerStyle(20);
242 // gr2->SetMarkerColor(i + 1);
243 // //gr2->SetMaximum(32);
244 // //gr2->SetMinimum(-14);
245 // gr2->SetMaximum(.4);
246 // gr2->SetMinimum(-.04);
247 // gr2->SetTitle("");
248 // gr2->GetXaxis()->SetTitle("vertical position (mm)");
249 // gr2->GetYaxis()->SetTitle("Cherenkov/MIP");
250 // if (i == 0)
251 // gr2->Draw("ap");
252 // if (i > 0)
253 // gr2->Draw("p");
254 // //gr2->GetXaxis()->TAxis::SetNdivisions(505);
255 // //gr2->GetXaxis()->SetLimits(160.6,166);
256 // gr2->GetXaxis()->SetLimits(130, 190);
257 // sprintf(name, "row %d", i);
258 // leg2->AddEntry(gr2, name, "p");
259 // }
260 // leg2->Draw();
261 
262 }
263