Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GetResolutions.C
Go to the documentation of this file.
1 // Script to extract time-walk constants for the BCAL
2 
3 #include <TFile.h>
4 #include <TDirectory.h>
5 
6 namespace GetResolutionsNS {
7  //Leave this global so the accesors don't need the pointer as an argument
8  TFile *thisFile=NULL;
9 
10  // Accessor functions to grab histograms from our file
11  // (Makes things easy with the HistogramTools fills)
12  TH1I * Get1DHistogram(const char * plugin, const char * directoryName, const char * name, bool print = true){
13  TH1I * histogram;
14  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
15  thisFile->GetObject(fullName, histogram);
16  if (histogram == 0){
17  if (print) cout << "Unable to find histogram " << fullName.Data() << endl;
18  return NULL;
19  }
20  return histogram;
21  }
22 
23  TH2I * Get2DHistogram(const char * plugin, const char * directoryName, const char * name){
24  TH2I * histogram;
25  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
26  thisFile->GetObject(fullName, histogram);
27  if (histogram == 0){
28  cout << "Unable to find histogram " << fullName.Data() << endl;
29  return NULL;
30  }
31  //cout << "Found histogram " << fullName.Data() << endl;
32  return histogram;
33  }
34 
35  TH2D * Get2DWeightedHistogram(const char * plugin, const char * directoryName, const char * name){
36  TH2D * histogram;
37  TString fullName = TString(plugin) + "/" + TString(directoryName) + "/" + TString(name);
38  thisFile->GetObject(fullName, histogram);
39  if (histogram == 0){
40  cout << "Unable to find histogram " << fullName.Data() << endl;
41  return NULL;
42  }
43  //cout << "Found histogram " << fullName.Data() << endl;
44  return histogram;
45  }
46 }
47 
48 void GetXYfromSize(int size, int &x, int &y) {
49  x=1;
50  y=1;
51  while (x*y<size) {
52  if (x==y) x++;
53  else
54  if (x>y) y++;
55  }
56 }
57 
58 
59 void GetResolutionHistograms(TH2 *hist2d, TH1D *rmshist, TH1D *fithist, TH1D *meanhist, TH1D *fitmeanhist, int minhits, TString tag, bool print);
60 int GetColor(int number);
61 void FormatHist(TH1* hist, int color, int markerstyle);
62 
63 void makeplot(std::string plot_title, TString tag, std::string plot_label,
64  std::vector<std::string> histname_plot,
65  std::vector<std::string> histdir_plot, std::string plot_type,
66  std::vector<std::string> legend_plot, bool logx=1) {
67  char canvasname[255];
68  sprintf(canvasname,"canvas_%s_%s",plot_label.c_str(),plot_type.c_str());
69  TCanvas *canvas_showers_sigma = new TCanvas(canvasname,canvasname);
70  gPad->SetLogx(logx);
71  gPad->SetGridx();
72  gPad->SetGridy();
73 
74  char fullhistname[255], plotfilename[255];
75  TLegend *leg_showers_sigma = new TLegend(0.75,0.75,0.99,0.99);
76  for (int i=0; i<histname_plot.size(); i++) {
77  sprintf(fullhistname,"%s_%s_%s",histdir_plot[i].c_str(),histname_plot[i].c_str(),plot_type.c_str());
78  TH1D* histpointer = (TH1D*)gDirectory->Get(fullhistname);
79  if (histpointer==NULL) {
80  printf(" no %s\n",fullhistname);
81  } else {
82  if (i==0) {
83  histpointer->SetTitle(plot_title.c_str());
84  histpointer->Draw();
85  }
86  else histpointer->Draw("same");
87  if (logx) histpointer->SetAxisRange(0.02,5);
88  leg_showers_sigma->AddEntry(histpointer,legend_plot[i].c_str(),"ep");
89  }
90  }
91  leg_showers_sigma->Draw();
92  sprintf(plotfilename,"plots/results/%s_%s_%s.png",plot_label.c_str(),plot_type.c_str(),tag.Data());
93  canvas_showers_sigma->Print(plotfilename);
94 }
95 
96 
97 // Do the extraction of the actual constants
98 void GetResolutions(int run = 2931, TString filename = "hd_root.root", TString tag = "default", bool summary=1){
99  printf("GetResolutions\n");
100  // Open our input and output file
101  printf("Opening for input \"%s\"\n",filename.Data());
102  //TFile *fil = new TFile(filename.Data());
103  GetResolutionsNS::thisFile = TFile::Open(filename.Data());
104  // Check to make sure it is open
105  if ( GetResolutionsNS::thisFile == NULL) {
106  cout << "Unable to open file " << filename.Data() << "...Exiting" << endl;
107  return;
108  } else {
109  cout << "Opened file " << filename.Data() << "...good" << endl;
110  }
111  //GetResolutionsNS::thisFile = fil;
112 
113  char outfilename[255];
114  sprintf(outfilename,"Resolutions_Results_%s.root",tag.Data());
115  printf("Opening for output %s\n",outfilename);
116  TFile *outputFile = TFile::Open(outfilename, "RECREATE");
117  // outputFile->mkdir("Fits");
118  outputFile->mkdir("ResultOverview");
119 
120 
121  gStyle->SetOptStat(2220);
122  gStyle->SetOptFit(1);
123 
124  int pad=0, color=0;
125 
126  // Showers
127 
128  std::vector<std::string> histname;
129  std::vector<std::string> histdir;
130  std::vector<int> minhits;
131 
132  // Make list of histograms to get resolutions from
133  histdir.push_back("Showers");
134  histname.push_back("PionShowers_q-");
135  minhits.push_back(1000);
136  histdir.push_back("Showers");
137  histname.push_back("AllShowers_q0");
138  minhits.push_back(1000);
139  if (!summary) { // only do all the fits not in summary mode
140  histdir.push_back("Points_deltaTVsEnergy");
141  histname.push_back("AllPoints_q-");
142  minhits.push_back(1000);
143  histdir.push_back("Points_deltaTVsEnergy");
144  histname.push_back("Layer1_q-");
145  minhits.push_back(1000);
146  histdir.push_back("Points_deltaTVsEnergy");
147  histname.push_back("Layer2_q-");
148  minhits.push_back(1000);
149  histdir.push_back("Points_deltaTVsEnergy");
150  histname.push_back("Layer3_q-");
151  minhits.push_back(1000);
152  histdir.push_back("Points_deltaTVsEnergy");
153  histname.push_back("Layer4_q-");
154  minhits.push_back(1000);
155  histdir.push_back("Points_deltaTVsEnergy");
156  histname.push_back("AllPoints_q0");
157  minhits.push_back(1000);
158  histdir.push_back("Points_deltaTVsEnergy");
159  histname.push_back("Layer1_q0");
160  minhits.push_back(1000);
161  histdir.push_back("Points_deltaTVsEnergy");
162  histname.push_back("Layer2_q0");
163  minhits.push_back(1000);
164  histdir.push_back("Points_deltaTVsEnergy");
165  histname.push_back("Layer3_q0");
166  minhits.push_back(1000);
167  histdir.push_back("Points_deltaTVsEnergy");
168  histname.push_back("Layer4_q0");
169  minhits.push_back(1000);
170  histdir.push_back("Points_deltaTVsShowerEnergy");
171  histname.push_back("AllPoints_q-");
172  minhits.push_back(1000);
173  histdir.push_back("Points_deltaTVsShowerEnergy");
174  histname.push_back("Layer1_q-");
175  minhits.push_back(1000);
176  histdir.push_back("Points_deltaTVsShowerEnergy");
177  histname.push_back("Layer2_q-");
178  minhits.push_back(1000);
179  histdir.push_back("Points_deltaTVsShowerEnergy");
180  histname.push_back("Layer3_q-");
181  minhits.push_back(1000);
182  histdir.push_back("Points_deltaTVsShowerEnergy");
183  histname.push_back("Layer4_q-");
184  minhits.push_back(1000);
185  histdir.push_back("Showers");
186  histname.push_back("PionShowersVsP_q-");
187  minhits.push_back(1000);
188  histdir.push_back("Showers");
189  histname.push_back("PionShowersVsZ_q-");
190  minhits.push_back(1000);
191  histdir.push_back("Showers");
192  histname.push_back("AllShowersVsZ_q0");
193  minhits.push_back(1000);
194 
195 
196  // histdir.push_back("Points_deltaTVsEnergy");
197  // histname.push_back("AllPoints_TDC_q-");
198  // minhits.push_back(1000);
199  // histdir.push_back("Points_deltaTVsEnergy");
200  // histname.push_back("AllPoints_ADC_q-");
201  // minhits.push_back(1000);
202  // histdir.push_back("Points_deltaTVsEnergy");
203  // histname.push_back("Layer1_TDC_q-");
204  // minhits.push_back(1000);
205  // histdir.push_back("Points_deltaTVsEnergy");
206  // histname.push_back("Layer2_TDC_q-");
207  // minhits.push_back(1000);
208  // histdir.push_back("Points_deltaTVsEnergy");
209  // histname.push_back("Layer3_TDC_q-");
210  // minhits.push_back(1000);
211  // histdir.push_back("Points_deltaTVsEnergy");
212  // histname.push_back("Layer1_ADC_q-");
213  // minhits.push_back(1000);
214  // histdir.push_back("Points_deltaTVsEnergy");
215  // histname.push_back("Layer2_ADC_q-");
216  // minhits.push_back(1000);
217  // histdir.push_back("Points_deltaTVsEnergy");
218  // histname.push_back("Layer3_ADC_q-");
219  // minhits.push_back(1000);
220  // histdir.push_back("Points_deltaTVsEnergy");
221  // histname.push_back("Layer4_ADC_q-");
222  // minhits.push_back(1000);
223  // histdir.push_back("Hits_deltaTVsE");
224  // histname.push_back("AllHits_TDC_q-");
225  // minhits.push_back(1000);
226  // histdir.push_back("Hits_deltaTVsE");
227  // histname.push_back("AllHits_ADC_q-");
228  // minhits.push_back(1000);
229 
230  histdir.push_back("Target Time");
231  histname.push_back("deltaTVsCell_q-");
232  minhits.push_back(10);
233  histdir.push_back("Target Time");
234  histname.push_back("deltaTVsCell_q-_Eweight");
235  minhits.push_back(10);
236  histdir.push_back("Target Time");
237  histname.push_back("deltaTVsCell_q0");
238  minhits.push_back(10);
239  histdir.push_back("Target Time");
240  histname.push_back("deltaTVsCell_q0_Eweight");
241  minhits.push_back(10);
242  }
243 
244  std::vector<TH2*> TimeVsEnergy;
245  std::vector<TH1D*> TimeSigma;
246  std::vector<TH1D*> TimeRMS;
247  std::vector<TH1D*> TimeMean;
248  std::vector<TH1D*> TimeFitMean;
249 
250  // Get the resolutions for the histograms
251  int xpads, ypads;
252  GetXYfromSize(histname.size(),xpads,ypads);
253  TCanvas *canvas_resolutions = new TCanvas("canvas_resolutions","canvas_resolutions",1200,800);
254  canvas_resolutions->Divide(xpads,ypads,0.001,0.001);
255  TCanvas *canvas_means = new TCanvas("canvas_means","canvas_means",1200,800);
256  canvas_means->Divide(xpads,ypads,0.001,0.001);
257  for (int i=0; i<histname.size(); i++) {
258  auto locTimeVsEnergy = GetResolutionsNS::Get2DHistogram("BCAL_Global_Offsets", histdir[i].c_str(), histname[i].c_str());
259  if (locTimeVsEnergy != NULL) {
260  TimeVsEnergy.push_back(locTimeVsEnergy);
261  } else {
262  printf("Looking for a weighted histogram\n");
263  auto locTimeVsEnergy = GetResolutionsNS::Get2DWeightedHistogram("BCAL_Global_Offsets", histdir[i].c_str(), histname[i].c_str());
264  if (locTimeVsEnergy != NULL) {
265  //printf("was able to recover\n");
266  TimeVsEnergy.push_back(locTimeVsEnergy);
267  }
268  }
269  //printf("size: %i\n",TimeVsEnergy.size());
270  if (TimeVsEnergy[i] == NULL) {
271  TimeSigma.push_back(NULL);
272  TimeRMS.push_back(NULL);
273  TimeMean.push_back(NULL);
274  TimeFitMean.push_back(NULL);
275  } else {
276  //printf("in loop %s %s %s\n","BCAL_Global_Offsets", histdir[i].c_str(), histname[i].c_str());
277  //gPad->SetLogx();
278  color = GetColor(i+1);
279  TAxis * xaxis = TimeVsEnergy[i]->GetXaxis();
280  int nBinsX = TimeVsEnergy[i]->GetNbinsX();
281  char fullhistname[255];
282  sprintf(fullhistname,"%s_%s_sigma",histdir[i].c_str(),histname[i].c_str());
283  //printf("%s\n",fullhistname);
284  TH1D* locTimeSigma = TimeVsEnergy[i]->ProjectionX(fullhistname);
285  TH1D* locTimeRMS = TimeVsEnergy[i]->ProjectionX(Form("%s_%s_RMS",histdir[i].c_str(),histname[i].c_str()));
286  TH1D* locTimeMean = TimeVsEnergy[i]->ProjectionX(Form("%s_%s_mean",histdir[i].c_str(),histname[i].c_str()));
287  TH1D* locTimeFitMean = TimeVsEnergy[i]->ProjectionX(Form("%s_%s_fitmean",histdir[i].c_str(),histname[i].c_str()));
288  TimeSigma.push_back(locTimeSigma);
289  TimeRMS.push_back(locTimeRMS);
290  TimeMean.push_back(locTimeMean);
291  TimeFitMean.push_back(locTimeFitMean);
292 
293  // TimeSigma[i] = new TH1D("TimeSigma","Charged shower time resolution;E [GeV];T_{Target}-t_{RF} [ns]",
294  // nBinsX,xaxis->GetXmin(),xaxis->GetXmax());
295  // TimeRMS[i] = new TH1D("TimeRMS","Charged shower time resolution;E [GeV];T_{Target}-t_{RF} [ns]",
296  // nBinsX,xaxis->GetXmin(),xaxis->GetXmax());
297  // TimeMean[i] = new TH1D("TimeMean","Charged shower time resolution;E [GeV];T_{Target}-t_{RF} [ns]",
298  // nBinsX,xaxis->GetXmin(),xaxis->GetXmax());
299  // TimeFitMean[i] = new TH1D("TimeFitMean","Charged shower time resolution;E [GeV];T_{Target}-t_{RF} [ns]",
300  // nBinsX,xaxis->GetXmin(),xaxis->GetXmax());
301 
302  FormatHist(TimeSigma[i],color,20);
303  TimeSigma[i]->GetYaxis()->SetTitle("#sigmat [ns]");
304  FormatHist(TimeRMS[i],color,22);
305  TimeRMS[i]->GetYaxis()->SetTitle("#sigmat [ns]");
306  FormatHist(TimeMean[i],color,20);
307  TimeMean[i]->GetYaxis()->SetTitle("#deltat [ns]");
308  FormatHist(TimeFitMean[i],color,22);
309  TimeFitMean[i]->GetYaxis()->SetTitle("#deltat [ns]");
310  printf("%2i: ",i);
311  TString newtag = tag + "_" + histname[i].c_str();
312  //printf("%s\n",newtag.Data());
313  GetResolutionHistograms(TimeVsEnergy[i],TimeRMS[i],TimeSigma[i],TimeMean[i],TimeFitMean[i],minhits[i],newtag.Data(),0);
314  canvas_resolutions->cd(i+1);
315  TimeRMS[i]->Draw();
316  TimeSigma[i]->Draw("same");
317  canvas_means->cd(i+1);
318  TimeFitMean[i]->Draw();
319  //TimeMean[i]->Draw("same");
320  }
321  }
322  char plotfilename[255];
323  if (!summary) {
324  sprintf(plotfilename,"plots/results/resolution_pad_%s.png",tag.Data());
325  canvas_resolutions->Print(plotfilename);
326  sprintf(plotfilename,"plots/results/means_pad_%s.png",tag.Data());
327  canvas_means->Print(plotfilename);
328  }
329 
330  /*
331  histdir_plot.push_back("Points_deltaTVsShowerEnergy");
332  histname_plot.push_back("AllPoints_q-");
333  legend_plot.push_back("");
334  histdir_plot.push_back("Points_deltaTVsShowerEnergy");
335  histname_plot.push_back("Layer1_q-");
336  legend_plot.push_back("");
337  histdir_plot.push_back("Points_deltaTVsShowerEnergy");
338  histname_plot.push_back("Layer2_q-");
339  legend_plot.push_back("");
340  histdir_plot.push_back("Points_deltaTVsShowerEnergy");
341  histname_plot.push_back("Layer3_q-");
342  legend_plot.push_back("");
343  histdir_plot.push_back("Points_deltaTVsShowerEnergy");
344  histname_plot.push_back("Layer4_q-");
345  legend_plot.push_back("");
346 */
347 
348  gStyle->SetOptStat(0);
349  gStyle->SetOptFit(0);
350  char title[255], fullhistname[255];
351  std::string title_plot;
352 
353  std::vector<std::string> histname_plot;
354  std::vector<std::string> histdir_plot;
355  std::vector<std::string> legend_plot;
356 
357  histdir_plot.push_back("Showers");
358  histname_plot.push_back("PionShowers_q-");
359  legend_plot.push_back("Charged showers");
360  histdir_plot.push_back("Showers");
361  histname_plot.push_back("AllShowers_q0");
362  legend_plot.push_back("Neutral showers");
363  title_plot = "BCAL Shower time resolution";
364  makeplot(title_plot,tag,"shower",histname_plot,histdir_plot,"sigma",legend_plot);
365  title_plot = "BCAL Shower time offset";
366  makeplot(title_plot,tag,"shower",histname_plot,histdir_plot,"fitmean",legend_plot);
367  histdir_plot.clear();
368  histname_plot.clear();
369  legend_plot.clear();
370 
371  if (!summary) {
372  histdir_plot.push_back("Points_deltaTVsEnergy");
373  histname_plot.push_back("AllPoints_q-");
374  legend_plot.push_back("All points");
375  histdir_plot.push_back("Points_deltaTVsEnergy");
376  histname_plot.push_back("Layer1_q-");
377  legend_plot.push_back("layer 1");
378  histdir_plot.push_back("Points_deltaTVsEnergy");
379  histname_plot.push_back("Layer2_q-");
380  legend_plot.push_back("layer 2");
381  histdir_plot.push_back("Points_deltaTVsEnergy");
382  histname_plot.push_back("Layer3_q-");
383  legend_plot.push_back("layer 3");
384  histdir_plot.push_back("Points_deltaTVsEnergy");
385  histname_plot.push_back("Layer4_q-");
386  legend_plot.push_back("layer 4");
387  title_plot = "Charged points time resolution";
388  makeplot(title_plot,tag,"charged_points",histname_plot,histdir_plot,"sigma",legend_plot);
389  title_plot = "Charged points time offset";
390  makeplot(title_plot,tag,"charged_points",histname_plot,histdir_plot,"fitmean",legend_plot);
391  histdir_plot.clear();
392  histname_plot.clear();
393  legend_plot.clear();
394 
395 
396  histdir_plot.push_back("Points_deltaTVsEnergy");
397  histname_plot.push_back("AllPoints_q0");
398  legend_plot.push_back("All points");
399  histdir_plot.push_back("Points_deltaTVsEnergy");
400  histname_plot.push_back("Layer1_q0");
401  legend_plot.push_back("Layer 1");
402  histdir_plot.push_back("Points_deltaTVsEnergy");
403  histname_plot.push_back("Layer2_q0");
404  legend_plot.push_back("Layer 2");
405  histdir_plot.push_back("Points_deltaTVsEnergy");
406  histname_plot.push_back("Layer3_q0");
407  legend_plot.push_back("Layer 3");
408  histdir_plot.push_back("Points_deltaTVsEnergy");
409  histname_plot.push_back("Layer4_q0");
410  legend_plot.push_back("Layer 4");
411  title_plot = "Neutral points time resolution";
412  makeplot(title_plot,tag,"neutral_points",histname_plot,histdir_plot,"sigma",legend_plot);
413  title_plot = "Neutral points time offset";
414  makeplot(title_plot,tag,"neutral_points",histname_plot,histdir_plot,"fitmean",legend_plot);
415  histdir_plot.clear();
416  histname_plot.clear();
417  legend_plot.clear();
418 
419  histdir_plot.push_back("Showers");
420  histname_plot.push_back("AllShowersVsZ_q0");
421  legend_plot.push_back("Neutral showers");
422  histdir_plot.push_back("Showers");
423  histname_plot.push_back("PionShowersVsZ_q-");
424  legend_plot.push_back("Charged showers");
425  title_plot = "Shower time resolution";
426  makeplot(title_plot,tag,"shower_vsZ",histname_plot,histdir_plot,"sigma",legend_plot,0);
427  title_plot = "Shower time offset";
428  makeplot(title_plot,tag,"shower_vsZ",histname_plot,histdir_plot,"fitmean",legend_plot,0);
429  histdir_plot.clear();
430  histname_plot.clear();
431  legend_plot.clear();
432 
433  histdir_plot.push_back("Showers");
434  histname_plot.push_back("PionShowersVsP_q-");
435  legend_plot.push_back("Charged showers");
436  title_plot = "Shower time resolution";
437  makeplot(title_plot,tag,"shower_vsP",histname_plot,histdir_plot,"sigma",legend_plot);
438  title_plot = "Shower time offset";
439  makeplot(title_plot,tag,"shower_vsP",histname_plot,histdir_plot,"fitmean",legend_plot);
440  histdir_plot.clear();
441  histname_plot.clear();
442  legend_plot.clear();
443  }
444 
445  outputFile->Write();
446  GetResolutionsNS::thisFile->Close();
447  return;
448 }
449 
450 
451 
452 
453 
454 
455 void GetResolutionHistograms(TH2 *hist2d, TH1D *rmshist, TH1D *fitsigmahist, TH1D *meanhist, TH1D *fitmeanhist, int minhits,
456  TString tag, bool print) {
457 
458  TAxis * xaxis = hist2d->GetXaxis();
459  TAxis * yaxis = hist2d->GetYaxis();
460  int nBinsX = hist2d->GetNbinsX();
461  int nBinsY = hist2d->GetNbinsY();
462 
463  printf("bins (%i,%i) %s %s\n",nBinsX,nBinsY,hist2d->GetName(),hist2d->GetTitle());
464 
465  TH1D *mulitbin_projY = hist2d->ProjectionY("mulitbin_projY", 1, 1); // structure to accumulate data bin by bin until enough to fit
466  mulitbin_projY->Reset();
467  TH1D *xmean = hist2d->ProjectionX("xmean"); // structure to accumulate average X value of bin
468  xmean->Reset();
469  // reset the histograms that will be returned
470  rmshist->Reset();
471  fitsigmahist->Reset();
472  meanhist->Reset();
473  fitmeanhist->Reset();
474 
475  TCanvas *canvas_fit = new TCanvas("canvas_fit","canvas_fit");
476  gPad->SetLogy();
477 
478  for (int i = 1 ; i <= nBinsX; i++){
479  TH1D *singlebin_projY = hist2d->ProjectionY("temp", i, i);
480  mulitbin_projY->Add(singlebin_projY); // accumulate bin
481  float center = yaxis->GetBinCenter(mulitbin_projY->GetMaximumBin());
482  // To avoid fitting the tails, only fit a region 1 ns around the peak
483  float width=1; // in ns
484  float mintime = center-width;
485  float maxtime = center+width;
486 
487  char title[255];
488  sprintf(title,"E=%.3f GeV (%s)",xaxis->GetBinCenter(i),tag.Data());
489  mulitbin_projY->SetTitle(title);
490 
491  int entries_multi = mulitbin_projY->GetEntries();
492  float int_multi = mulitbin_projY->Integral(yaxis->FindBin(mintime),yaxis->FindBin(maxtime));
493  float int_single = singlebin_projY->Integral(yaxis->FindBin(mintime),yaxis->FindBin(maxtime));
494  xmean->SetBinContent(i,int_single); // used to find weighted center of new composite bin
495  float weightedmean = xmean->GetMean();
496  // int bin = xmean->GetXaxis()->FindBin(xmean);
497  int bin = xmean->FindBin(weightedmean);
498  //if (int_multi<minhits || (i==nBinsX && int_multi<(minhits/4.))) {
499  if (entries_multi<minhits || (i==nBinsX && entries_multi<(minhits/4.))) {
500  //printf("%4i %4i %8.1f %8.1f %6i %6i %7.2f %2i\n",i,nBinsX,int_single,int_multi,entries_multi,minhits,weightedmean,bin);
501  continue;
502  }
503  //if (int_multi != int_single && i!=bin) printf("%4i %i %i %.2f %i\n", i,int_single,int_multi,weightedmean,bin);
504 
505  mulitbin_projY->SetAxisRange(-5,5);
506  double mean = mulitbin_projY->GetMean();
507  double meanerr = mulitbin_projY->GetMeanError();
508  double RMS = mulitbin_projY->GetRMS();
509  double RMSerr = mulitbin_projY->GetRMSError();
510  rmshist->SetBinContent(bin,RMS);
511  rmshist->SetBinError(bin,RMSerr);
512  meanhist->SetBinContent(bin,mean);
513  meanhist->SetBinError(bin,meanerr);
514 
515  mulitbin_projY->SetAxisRange(mintime,maxtime);
516  TF1 *f_gaus = new TF1("f_gaus","gaus",mintime,maxtime);
517  f_gaus->SetLineColor(kRed);
518  f_gaus->SetLineWidth(1);
519 
520  mulitbin_projY->Fit(f_gaus,"RQ");
521  double fitmean = f_gaus->GetParameter(1);
522  double fitmeanerr = f_gaus->GetParError(1);
523  double sigma = f_gaus->GetParameter(2);
524  double sigmaerr = f_gaus->GetParError(2);
525  //printf("mean %6.3f %.3f RMS %6.3f %.3f fit mean %6.3f %.3f sigma %6.3f %.3f\n",
526  // mean,meanerr,sigma,sigmaerr,fitmean,fitmeanerr,sigma,sigmaerr);
527  fitsigmahist->SetBinContent(bin,sigma);
528  fitsigmahist->SetBinError(bin,sigmaerr);
529  fitmeanhist->SetBinContent(bin,fitmean);
530  fitmeanhist->SetBinError(bin,fitmeanerr);
531  if (print) {
532  char plotname[255];
533  gStyle->SetPadRightMargin(0.20);
534  gStyle->SetStatW(0.20);
535  mulitbin_projY->SetAxisRange(-10,10);
536  sprintf(plotname,"plots/resolutions/fit_%s_%i.png",tag.Data(),i);
537  canvas_fit->Print(plotname);
538  }
539  mulitbin_projY->Reset();
540  xmean->Reset();
541  }
542  delete canvas_fit;
543 }
544 
545 
546 
547 int GetColor(int number)
548 {
549  int color = number;
550  if (color > 4) color++;
551  if (color > 9) color = 791-9+number;
552  //printf("%i %i\n",number,color);
553  return color;
554 }
555 
556 void FormatHist(TH1* hist, int color, int markerstyle)
557 {
558  hist->SetMarkerStyle(markerstyle);
559  hist->SetMarkerSize(0.8);
560  hist->SetMarkerColor(color);
561  hist->SetLineColor(color);
562  hist->SetLineWidth(2);
563 }
hist2d
TH1F * hist2d[Idx+1]
Definition: readhist.C:16
GetResolutions
void GetResolutions(int run=2931, TString filename="hd_root.root", TString tag="default", bool summary=1)
Definition: GetResolutions.C:98
GetResolutionHistograms
void GetResolutionHistograms(TH2 *hist2d, TH1D *rmshist, TH1D *fithist, TH1D *meanhist, TH1D *fitmeanhist, int minhits, TString tag, bool print)
Definition: GetResolutions.C:455
pad
TPad * pad
Definition: psc_mon.C:73
x
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
string
char string[256]
Definition: plot_epics_dump.C:24
sprintf
sprintf(text,"Post KinFit Cut")
filename
TString filename
Definition: Read_bcal_hadronic_eff2.C:55
y
#define y
run
int run
Definition: JEventProcessor_PS_E_calib.cc:60
GetResolutionsNS::Get2DWeightedHistogram
TH2D * Get2DWeightedHistogram(const char *plugin, const char *directoryName, const char *name)
Definition: GetResolutions.C:35
size
Definition: DFDCPseudo.h:53
FormatHist
void FormatHist(TH1 *hist, int color, int markerstyle)
Definition: GetResolutions.C:556
GetColor
int GetColor(int number)
Definition: GetResolutions.C:547
GetResolutionsNS::Get2DHistogram
TH2I * Get2DHistogram(const char *plugin, const char *directoryName, const char *name)
Definition: GetResolutions.C:23
sigma
Double_t sigma[NCHANNELS]
Definition: st_tw_resols.C:37
GetXYfromSize
void GetXYfromSize(int size, int &x, int &y)
Definition: GetResolutions.C:48
GetResolutionsNS::Get1DHistogram
TH1I * Get1DHistogram(const char *plugin, const char *directoryName, const char *name, bool print=true)
Definition: GetResolutions.C:12
makeplot
void makeplot(std::string plot_title, TString tag, std::string plot_label, std::vector< std::string > histname_plot, std::vector< std::string > histdir_plot, std::string plot_type, std::vector< std::string > legend_plot, bool logx=1)
Definition: GetResolutions.C:63
printf
printf("string=%s", string)
hist
TH1F * hist[Idx+1]
Definition: readhist.C:10