Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
plot_checks.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file plot_checks.C
1 #include "get_runstr.h"
2 #include <iostream>
3 #include <fstream>
4 #include <TGraphErrors.h>
5 #include <TString.h>
6 
7 int do_zdc = 0;
8 int do_emcal = 0;
9 int do_ohcal = 0;
10 int do_ihcal = 0;
11 int calc_gain_corr = 1;
12 
13 const int MAXRUNS = 20000;
14 Double_t bqmeansum[MAXRUNS]; // bq mean in sum of two arms
15 Double_t bqmeansumerr[MAXRUNS];
16 Double_t bqmeantot[2][MAXRUNS]; // bq mean in each arm
17 Double_t bqmeantoterr[2][MAXRUNS];
18 Double_t bqmean[128][MAXRUNS];
19 Double_t bqmeanerr[128][MAXRUNS];
20 Double_t runindex[MAXRUNS];
21 TGraphErrors *g_meansum {nullptr};
22 TGraphErrors *g_meantot[2] {nullptr};
23 TGraphErrors *g_mean[128] {nullptr};
24 
25 TFile *tfile[MAXRUNS];
26 TH1 *h_bbcqsum[MAXRUNS];
27 TH1 *h_bbcqtot[2][MAXRUNS];
28 TH1 *h_bbcq[MAXRUNS][128];
29 TH2 *h2_bbcq[128];
30 
31 TH1 *h_zdce[MAXRUNS];
32 TH1 *h_emcale[MAXRUNS];
33 TH1 *h_ohcale[MAXRUNS];
34 TH1 *h_ihcale[MAXRUNS];
35 
36 
37 // Get mean gain for each channel across runs, and write out the values
38 void fit_across_runs(TGraphErrors *gainvals[], const char *flist)
39 {
40  TString meangainfname = "results/"; meangainfname.Append( flist );
41  meangainfname.ReplaceAll(".list","_qmean.calib");
42  ofstream meangainfile( meangainfname );
43  TF1 *f_meangain[128] {nullptr};
44  Double_t grp_mean[128];
45  Double_t grp_meanerr[128];
46  for (int ipmt=0; ipmt<128; ipmt++)
47  {
48  TString name = "f_meangain"; name += ipmt;
49  f_meangain[ipmt] = new TF1(name,"pol0",0,1e9);
50  gainvals[ipmt]->Fit(f_meangain[ipmt]);
51  gPad->Update();
52  gPad->Modified();
53  grp_mean[ipmt] = f_meangain[ipmt]->GetParameter(0);
54  grp_meanerr[ipmt] = f_meangain[ipmt]->GetParError(0);
55 
56  meangainfile << ipmt << "\t" << grp_mean[ipmt] << "\t" << grp_meanerr[ipmt] << endl;
57  cout << ipmt << "\t" << grp_mean[ipmt] << "\t" << grp_meanerr[ipmt] << endl;
58  }
59  meangainfile.close();
60 }
61 
62 void plot_checks(const char *flistname = "c.list")
63 {
64 
65  ifstream flist(flistname);
66 
67  TString dir = "results/plot_checks_"; dir += flistname;
68  dir.ReplaceAll(".list","/");
69  TString name = "mkdir -p " + dir;
70  gSystem->Exec( name );
71 
72  TString dstfname;
73  int nruns = 0;
74  while ( flist >> dstfname )
75  {
76  tfile[nruns] = new TFile(dstfname,"READ");
77  h_bbcqsum[nruns] = (TH1*)tfile[nruns]->Get("h_bbcqsum");
78  h_bbcqtot[0][nruns] = (TH1*)tfile[nruns]->Get("h_bbcqtot0");
79  h_bbcqtot[1][nruns] = (TH1*)tfile[nruns]->Get("h_bbcqtot1");
80  bqmeansum[nruns] = h_bbcqsum[nruns]->GetMean();
81  bqmeansumerr[nruns] = h_bbcqsum[nruns]->GetMeanError();
82  bqmeantot[0][nruns] = h_bbcqtot[0][nruns]->GetMean();
83  bqmeantoterr[0][nruns] = h_bbcqtot[0][nruns]->GetMeanError();
84  bqmeantot[1][nruns] = h_bbcqtot[0][nruns]->GetMean();
85  bqmeantoterr[1][nruns] = h_bbcqtot[0][nruns]->GetMeanError();
86 
87  h_zdce[nruns] = (TH1*)tfile[nruns]->Get("h_zdce");
88  h_emcale[nruns] = (TH1*)tfile[nruns]->Get("h_emcale");
89  h_ohcale[nruns] = (TH1*)tfile[nruns]->Get("h_ohcale");
90  h_ihcale[nruns] = (TH1*)tfile[nruns]->Get("h_ihcale");
91 
92  for (int ipmt=0; ipmt<128; ipmt++)
93  {
94  name = "h_bbcq"; name += ipmt;
95  h_bbcq[nruns][ipmt] = (TH1*)tfile[nruns]->Get(name);
96  bqmean[ipmt][nruns] = h_bbcq[nruns][ipmt]->GetMean();
97  bqmeanerr[ipmt][nruns] = h_bbcq[nruns][ipmt]->GetMeanError();
98  }
99 
100  runindex[nruns] = nruns;
101  nruns++;
102  }
103  cout << "Processed " << nruns << " runs." << endl;
104 
105  TCanvas *bc = new TCanvas("bc","mean fits",800,600);
106  TString title;
107  TString meangainfname = dir + flistname;
108  meangainfname.ReplaceAll(".list","_qmean.calib");
109  ofstream meanfile( meangainfname );
110  TF1 *meanfit = new TF1("meanfit","pol0",0,nruns);
111  double runmean[128];
112  double runmeanerr[128];
113  double chi2ndf[128];
114  for (int ipmt=0; ipmt<128; ipmt++)
115  {
116  g_mean[ipmt] = new TGraphErrors(nruns,runindex,bqmean[ipmt],0,bqmeanerr[ipmt]);
117  name = "q_mean"; name += ipmt;
118  title = name;
119  g_mean[ipmt]->SetName( name );
120  g_mean[ipmt]->SetTitle( title );
121 
122  g_mean[ipmt]->Draw("ap");
123  meanfit->SetParameter(0,100);
124  g_mean[ipmt]->Fit(meanfit,"R");
125  gPad->Modified();
126  gPad->Update();
127 
128  runmean[ipmt] = meanfit->GetParameter(0);
129  runmeanerr[ipmt] = meanfit->GetParError(0);
130  Double_t chi2 = meanfit->GetChisquare();
131  Double_t ndf = meanfit->GetNDF();
132  chi2ndf[ipmt] = chi2/ndf;
133 
134  meanfile << ipmt << "\t" << runmean[ipmt] << "\t" << runmeanerr[ipmt] << "\t" << chi2ndf[ipmt] << endl;
135  cout << ipmt << "\t" << runmean[ipmt] << "\t" << runmeanerr[ipmt] << "\t" << chi2ndf[ipmt] << endl;
136  }
137  meanfile.close();
138 
139  // Make vs run index plots
140  TH2 *h2_bbcqsum_vs_run = new TH2F("h2_bbcqsum_vs_run","bbc qsum vs run",nruns,0,nruns,3000,0,3000);
141  TH2 *h2_bbcqs_vs_run = new TH2F("h2_bbcqs_vs_run","bbc south qsum vs run",nruns,0,nruns,1400,0,1400);
142  TH2 *h2_bbcqn_vs_run = new TH2F("h2_bbcqn_vs_run","bbc north qsum vs run",nruns,0,nruns,1400,0,1400);
143  for (int ipmt=0; ipmt<128; ipmt++)
144  {
145  name = "h2_bbcq"; name += ipmt;
146  h2_bbcq[ipmt] = new TH2F(name,name,nruns,0,nruns,1200,0,60);
147  }
148  TH2 *h2_zdce_vs_run {nullptr};
149  TH2 *h2_emcale_vs_run {nullptr};
150  TH2 *h2_ohcale_vs_run {nullptr};
151  TH2 *h2_ihcale_vs_run {nullptr};
152  if ( do_zdc )
153  {
154  int nbinsx = h_zdce[0]->GetNbinsX();
155  Double_t xmin = h_zdce[0]->GetBinLowEdge(1);
156  Double_t xmax = h_zdce[0]->GetBinLowEdge(nbinsx+1);
157  h2_zdce_vs_run = new TH2F("h2_zdce_vs_run","ZDC E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
158  }
159  if ( do_emcal )
160  {
161  int nbinsx = h_emcale[0]->GetNbinsX();
162  Double_t xmin = h_emcale[0]->GetBinLowEdge(1);
163  Double_t xmax = h_emcale[0]->GetBinLowEdge(nbinsx+1);
164  h2_emcale_vs_run = new TH2F("h2_emcale_vs_run","EMC E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
165  }
166  if ( do_ohcal )
167  {
168  int nbinsx = h_ohcale[0]->GetNbinsX();
169  Double_t xmin = h_ohcale[0]->GetBinLowEdge(1);
170  Double_t xmax = h_ohcale[0]->GetBinLowEdge(nbinsx+1);
171  h2_ohcale_vs_run = new TH2F("h2_ohcale_vs_run","OHCAL E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
172  }
173  if ( do_ihcal )
174  {
175  int nbinsx = h_ihcale[0]->GetNbinsX();
176  Double_t xmin = h_ihcale[0]->GetBinLowEdge(1);
177  Double_t xmax = h_ihcale[0]->GetBinLowEdge(nbinsx+1);
178  h2_ihcale_vs_run = new TH2F("h2_ihcale_vs_run","iHCAL E vs run",nruns,0,nruns,nbinsx,xmin,xmax);
179  }
180 
181 
182  for (int irun=0; irun<nruns; irun++)
183  {
184  int nbinsx = h_bbcqsum[irun]->GetNbinsX();
185  for (int ibin=1; ibin<=nbinsx; ibin++)
186  {
187  Float_t val = h_bbcqsum[irun]->GetBinContent(ibin);
188  Float_t qsum = h_bbcqsum[irun]->GetBinCenter(ibin);
189  h2_bbcqsum_vs_run->Fill( irun, qsum, val );
190 
191  val = h_bbcqtot[0][irun]->GetBinContent(ibin);
192  qsum = h_bbcqtot[0][irun]->GetBinCenter(ibin);
193  h2_bbcqs_vs_run->Fill( irun, qsum, val );
194 
195  val = h_bbcqtot[1][irun]->GetBinContent(ibin);
196  qsum = h_bbcqtot[1][irun]->GetBinCenter(ibin);
197  h2_bbcqn_vs_run->Fill( irun, qsum, val );
198  }
199 
200  for (int ipmt=0; ipmt<128; ipmt++)
201  {
202  int nbinsx = h_bbcq[irun][ipmt]->GetNbinsX();
203  for (int ibin=1; ibin<=nbinsx; ibin++)
204  {
205  Double_t val = h_bbcq[irun][ipmt]->GetBinContent(ibin);
206  Double_t q = h_bbcq[irun][ipmt]->GetBinCenter(ibin);
207  h2_bbcq[ipmt]->Fill( (Double_t)irun, q, val );
208  }
209  }
210 
211  if ( do_zdc )
212  {
213  int nbinsx = h_zdce[irun]->GetNbinsX();
214  for (int ibin=1; ibin<=nbinsx; ibin++)
215  {
216  Float_t val = h_zdce[irun]->GetBinContent(ibin);
217  Float_t e = h_zdce[irun]->GetBinCenter(ibin);
218  h2_zdce_vs_run->Fill( irun, e, val );
219  }
220  }
221 
222  if ( do_emcal )
223  {
224  int nbinsx = h_emcale[irun]->GetNbinsX();
225  for (int ibin=1; ibin<=nbinsx; ibin++)
226  {
227  Float_t val = h_emcale[irun]->GetBinContent(ibin);
228  Float_t emce = h_emcale[irun]->GetBinCenter(ibin);
229  h2_emcale_vs_run->Fill( irun, emce, val );
230  }
231  }
232 
233  if ( do_ohcal )
234  {
235  int nbinsx = h_ohcale[irun]->GetNbinsX();
236  for (int ibin=1; ibin<=nbinsx; ibin++)
237  {
238  Float_t val = h_ohcale[irun]->GetBinContent(ibin);
239  Float_t e = h_ohcale[irun]->GetBinCenter(ibin);
240  h2_ohcale_vs_run->Fill( irun, e, val );
241  }
242  }
243 
244  if ( do_ihcal )
245  {
246  int nbinsx = h_ihcale[irun]->GetNbinsX();
247  for (int ibin=1; ibin<=nbinsx; ibin++)
248  {
249  Float_t val = h_ihcale[irun]->GetBinContent(ibin);
250  Float_t e = h_ihcale[irun]->GetBinCenter(ibin);
251  h2_ihcale_vs_run->Fill( irun, e, val );
252  }
253  }
254 
255 
256  }
257 
258  TCanvas *ac[100];
259  int icv = 0;
260 
261  ac[icv] = new TCanvas("mbdq_vs_run2d","MBD Q vs Run",800,600);
262  h2_bbcqsum_vs_run->SetMinimum(1e-5);
263  h2_bbcqsum_vs_run->DrawCopy("colz");
264  gPad->SetLogz(1);
265 
266  ac[++icv] = new TCanvas("mbdqs_vs_run2d","MBD Q.S vs Run",800,600);
267  h2_bbcqs_vs_run->SetMinimum(1e-6);
268  h2_bbcqs_vs_run->DrawCopy("colz");
269  gPad->SetLogz(1);
270 
271  ac[++icv] = new TCanvas("mbdqn_vs_run2d","MBD Q.N vs Run",800,600);
272  h2_bbcqn_vs_run->SetMinimum(1e-6);
273  h2_bbcqn_vs_run->DrawCopy("colz");
274  gPad->SetLogz(1);
275 
276  ac[++icv] = new TCanvas("mbdq_vs_run","MBD Qsum vs Run",800,600);
277  g_meansum = new TGraphErrors(nruns,runindex,bqmeansum,0,bqmeansumerr);
278  name = "q_meansum";
279  title = name;
280  g_meansum->SetName( name );
281  g_meansum->SetTitle( title );
282  g_meansum->SetMarkerStyle(20);
283  g_meansum->Draw("ap");
284 
285  ac[++icv] = new TCanvas("mbdqs_vs_run","MBD Q.S vs Run",800,600);
286  g_meantot[0] = new TGraphErrors(nruns,runindex,bqmeantot[0],0,bqmeantoterr[0]);
287  name = "q_meantot0";
288  title = name;
289  g_meantot[0]->SetName( name );
290  g_meantot[0]->SetTitle( title );
291  g_meantot[0]->SetMarkerStyle(20);
292  g_meantot[0]->Draw("ap");
293 
294  ac[++icv] = new TCanvas("mbdqn_vs_run","MBD Q.N vs Run",800,600);
295  g_meantot[1] = new TGraphErrors(nruns,runindex,bqmeantot[1],0,bqmeantoterr[1]);
296  name = "q_meantot1";
297  title = name;
298  g_meantot[1]->SetName( name );
299  g_meantot[1]->SetTitle( title );
300  g_meantot[1]->SetMarkerStyle(20);
301  g_meantot[1]->Draw("ap");
302 
303  ac[++icv] = new TCanvas("mbdq_vs_run_ch","MBD Q vs Run",1200,800);
304  for (int ipmt=0; ipmt<128; ipmt++)
305  {
306  h2_bbcq[ipmt]->Draw("colz");
307  gPad->SetLogz(1);
308 
309  name = dir; name += h2_bbcq[ipmt]->GetName(); name += ".png";
310  //cout << name << endl;
311  ac[icv]->Print( name );
312  }
313 
314  /*
315  ac[++icv] = new TCanvas("mbdq_vs_run_ch_low","MBD Q vs Run, Low",1200,800);
316  for (int ipmt=0; ipmt<128; ipmt++)
317  {
318  h2_bbcq[ipmt]->GetYaxis()->SetRangeUser(0,4);
319  h2_bbcq[ipmt]->Draw("colz");
320 
321  name = dir; name += h2_bbcq[ipmt]->GetName(); name += "_low.png";
322  //cout << name << endl;
323  ac[icv]->Print( name );
324  }
325  */
326 
327  int run1 = 0;
328  int run2 = 1;
329  /*
330  h_bbcqsum[run1]->Rebin(8);
331  h_bbcqsum[run2]->Rebin(8);
332  h_bbcqsum[run1]->Draw("hist");
333  h_bbcqsum[run2]->SetLineColor(2);
334  h_bbcqsum[run2]->Draw("histsame");
335  gPad->SetLogy(1);
336  */
337 
338  if ( do_zdc )
339  {
340  ac[++icv] = new TCanvas("zdce_vs_run","ZDC E vs Run",800,600);
341  h2_zdce_vs_run->SetMinimum(1e-6);
342  h2_zdce_vs_run->Draw("colz");
343  gPad->SetLogz(1);
344  }
345 
346  if ( do_emcal )
347  {
348  ac[++icv] = new TCanvas("emce_vs_run","EMC E vs Run",800,600);
349  h2_emcale_vs_run->SetMinimum(1e-6);
350  h2_emcale_vs_run->Draw("colz");
351  gPad->SetLogz(1);
352 
353  h_emcale[run1]->Rebin(4);
354  h_emcale[run2]->Rebin(4);
355  h_emcale[run2]->Scale(0.6);
356  h_emcale[run1]->Draw("hist");
357  h_emcale[run2]->SetLineColor(2);
358  h_emcale[run2]->Draw("histsame");
359  gPad->SetLogy(1);
360  }
361 
362  if ( do_ohcal )
363  {
364  ac[++icv] = new TCanvas("ohcal_vs_run","oHCAL E vs Run",800,600);
365  h2_ohcale_vs_run->SetMinimum(1e-6);
366  h2_ohcale_vs_run->Draw("colz");
367  gPad->SetLogz(1);
368 
369  h_ohcale[run1]->Rebin(4);
370  h_ohcale[run2]->Rebin(4);
371  h_ohcale[run2]->Scale(0.6);
372  h_ohcale[run1]->Draw("hist");
373  h_ohcale[run2]->SetLineColor(2);
374  h_ohcale[run2]->Draw("histsame");
375  gPad->SetLogy(1);
376  }
377 
378  if ( do_ihcal )
379  {
380  ac[++icv] = new TCanvas("ihcal_vs_run","iHCAL E vs Run",800,600);
381  h2_ihcale_vs_run->SetMinimum(1e-6);
382  h2_ihcale_vs_run->Draw("colz");
383  gPad->SetLogz(1);
384 
385  h_ihcale[run1]->Rebin(4);
386  h_ihcale[run2]->Rebin(4);
387  h_ihcale[run2]->Scale(0.6);
388  h_ihcale[run1]->Draw("hist");
389  h_ihcale[run2]->SetLineColor(2);
390  h_ihcale[run2]->Draw("histsame");
391  gPad->SetLogy(1);
392  }
393 
394 
395  // find gain corrections
396  if ( calc_gain_corr )
397  {
398  //int refrun = 0; // index of reference run
399  flist.clear();
400  flist.seekg(0);
401  int irun = 0;
402  ofstream gaincorr_calfile;
403  while ( flist >> dstfname )
404  {
405  TString dir = "results/";
406  dir += get_runstr(dstfname);
407  dir += "/";
408  name = "mkdir -p " + dir;
409  gSystem->Exec( name );
410  name = dir + "mbd_gaincorrmean.calib";
411  gaincorr_calfile.open( name );
412 
413  cout << name << endl;
414 
415  for (int ipmt=0; ipmt<128; ipmt++)
416  {
417  double corr = bqmean[ipmt][irun] / runmean[ipmt];
418  //double corr = bqmean[ipmt][irun] / bqmean[refun][ipmt];
419  if ( fabs(corr-1.0)>0.01 )
420  {
421  cout << irun << "\t" << ipmt << "\t" << corr << endl;
422  }
423  gaincorr_calfile << ipmt << "\t" << corr << endl;
424  }
425 
426  gaincorr_calfile.close();
427 
428  irun++;
429  }
430  }
431 
432 }