PlotStyle.h
Go to the documentation of this file.
1 #pragma once
2 
3 // CAFAna includes
4 #include "CAFAna/Core/Ratio.h"
5 #include "CAFAna/Core/Spectrum.h"
10 #include "Utilities/rootlogon.C"
11 
12 // ROOT includes
13 #include "TLine.h"
14 #include "TCanvas.h"
15 #include "TPad.h"
16 #include "TGraph.h"
17 #include "TGaxis.h"
18 #include "TString.h"
19 #include "TLatex.h"
20 #include "TLegend.h"
21 #include "TH1.h"
22 #include "TPaveStats.h"
23 #include "TF1.h"
24 #include "TMath.h"
25 #include "TStatistic.h"
26 #include "TStyle.h"
27 #include "TH2.h"
28 #include "TProfile.h"
29 #include "TPaletteAxis.h"
30 #include "TColor.h"
31 #include "TList.h"
32 #include "TFrame.h"
33 #include "TObjArray.h"
34 #include "TObject.h"
35 #include "TArrow.h"
36 #include "THStack.h"
37 #include "TROOT.h"
38 #include "TRint.h"
39 #include "TLegendEntry.h"
40 
41 // C++ includes
42 #include <iostream>
43 #include <vector>
44 #include <iostream>
45 #include <fstream>
46 #include <iomanip>
47 #include <sstream>
48 #include <string>
49 #include <cstring>
50 #include "TFile.h"
51 #include "TH2D.h"
52 #include "TObject.h"
53 #include "TArray.h"
54 #include "TArrayF.h"
55 #include "TGraphAsymmErrors.h"
56 
57 void MuonFit(ana::Spectrum &plot1, TString zaxisname,
58  TString outname1, TString outname2, TString outname3,
59  Bool_t act=false, Bool_t actandcat=false, Bool_t cat=false,
60  Bool_t printme=false){
61 
62  double pot = 8.09e20;
63  TH2* hist = plot1.ToTH2(pot);
64  if (act){
65  hist->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
66  hist->GetYaxis()->SetTitle("True Muon Energy (GeV)");
67  }
68  if (actandcat){
69  hist->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
70  hist->GetYaxis()->SetTitle("True Muon E - True Muon E in Catcher (GeV)");
71  hist->GetYaxis()->SetTitleSize(0.048);
72  }
73  if (cat){
74  hist->GetXaxis()->SetTitle("Reco Muon Track Length in Catcher (m)");
75  hist->GetYaxis()->SetTitle("True Muon Energy in Catcher (GeV)");
76  }
77 
78  hist->GetZaxis()->SetTitle(zaxisname);
79  static const int nn = 150; //Number of points on graph - should be equal to or smaller than number of bins in x
80  float xxx [nn];
81  float exxl[nn];
82  float exxh[nn];
83  float yyx [nn];
84  float eyxl[nn];
85  float eyxh[nn];
86 
87  // TH1D* chiValues = new TH1D("chiValue",";Chi squared per NDF;Fits",150, 0.0, 5.0);
88 
89  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
90  // make projection in Y view
91  TH1D* py = hist->ProjectionY("py",iw,iw,"");
92  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
93  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
94 
95  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
96 
97  double totalEntries = py->GetEntries();
98  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
99 
100  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
101 
102  //Find the location of width at half max
103  int yyhighbin = maxbin+1;
104  int yylowbin = maxbin-1;
105 
106  // find upper bin of half max
107  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
108  yyhighbin = ihigh;
109  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
110  }
111 
112  // find lower bin of half max
113  for(int ilow = maxbin-1; ilow > 0; --ilow){
114  yylowbin = ilow;
115  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
116  }
117 
118  float multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
119  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
120  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
121 
122  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
123  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
124 
125  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
126 
127  TF1 *f1 = new TF1("customGaus","gaus",0,5);
128  py->Fit("customGaus","OQ","",lowSide,highSide); //O and Q mean do not plot result, quiet mode
129  double mean = f1->GetParameter(1);
130  double error = f1->GetParError(1);
131  double chi = f1->GetChisquare();
132  double ndf = f1->GetNDF();
133 
134  // if (ndf > 0) {chiValues->Fill(chi/ndf);}
135  //else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
136  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
137  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
138  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
139  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
140  //Setting the error on the mean reported by the fit as error in y
141  eyxl[iw-1] = error;
142  eyxh[iw-1] = error;
143 
144  }//End of loop over bins in 2D histogram
145  CenterTitles(hist);
146 
147  //->GetZaxis()->CenterTitle();
148  //Making graph from points
149  TGraphAsymmErrors* customProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
150  customProfile->GetXaxis()->SetNoExponent(kTRUE);
151  customProfile->GetXaxis()->CenterTitle();
152  customProfile->GetYaxis()->CenterTitle();
153  if (act){
154  customProfile->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
155  customProfile->GetYaxis()->SetTitle("True Muon Energy (GeV)");
156  customProfile->GetXaxis()->SetLimits(0,15);
157  customProfile->GetHistogram()->SetMaximum(5);
158  customProfile->GetHistogram()->SetMinimum(0);
159  }
160  if (actandcat){
161  customProfile->GetXaxis()->SetTitle("Reco Muon Track Length (m)");
162  customProfile->GetYaxis()->SetTitle("True Muon E - True Muon E in Catcher (GeV)");
163  customProfile->GetXaxis()->SetLimits(0,15);
164  customProfile->GetHistogram()->SetMaximum(5);
165  customProfile->GetHistogram()->SetMinimum(0);
166  customProfile->GetYaxis()->SetTitleSize(0.048);
167  }
168  if (cat){
169  customProfile->GetXaxis()->SetTitle("Reco Muon Track Length in Catcher (m)");
170  customProfile->GetYaxis()->SetTitle("True Muon Energy in Catcher (GeV)");
171  customProfile->GetXaxis()->SetLimits(0,3);
172  customProfile->GetHistogram()->SetMaximum(3);
173  customProfile->GetHistogram()->SetMinimum(0);
174  }
175 
176  //CenterTitles(customProfile);
177  customProfile->SetTitle("");
178 
179  //Overlays 2d plot and graph
180  TCanvas* c1 = new TCanvas("c1","c1");
181  customProfile->SetMarkerColor(18);
182  customProfile->SetLineColor(18);
183  hist->Draw("colz");
184  customProfile->Draw("P");
185  Simulation();
186  c1->SetRightMargin(0.15);
187  c1->SetLeftMargin(0.15);
188 
189  if (printme){
190  TString outDir = "plots";
191  c1->SaveAs(outDir+"/"+outname1+".png");
192  c1->SaveAs(outDir+"/"+outname1+".pdf");
193  }
194 
195  //Overlays 2d loqz plot and graph
196  TCanvas* c2 = new TCanvas("c2","c2");
197  c2->SetLogz();
198  customProfile->SetMarkerColor(13);
199  customProfile->SetLineColor(13);
200  hist->Draw("colz");
201  customProfile->Draw("P");
202  Simulation();
203  c2->SetRightMargin(0.15);
204  c2->SetLeftMargin(0.15);
205  if (printme){
206  TString outDir = "plots";
207  c2->SaveAs(outDir+"/"+outname2+".png");
208  c2->SaveAs(outDir+"/"+outname2+".pdf");
209  }
210 
211  //Just the graph
212  TCanvas* c3 = new TCanvas("c3","c3");
213  customProfile->SetMarkerColor(1);
214  customProfile->SetLineColor(1);
215  customProfile->Draw("AP");
216  c3->SetRightMargin(0.15);
217  if (act){
218  customProfile->Fit("pol7", "","",0.2, 12.0);
219  }
220  if (actandcat){
221  customProfile->Fit("pol7", "","",3.75, 12.2);
222  }
223 
224  if (cat){
225  customProfile->Fit("pol1", "","",0.02, 2.3);
226  }
227 
228  c3->SetLeftMargin(0.15);
229 
230  if (printme){
231  TString outDir = "plots";
232  c3->SaveAs(outDir+"/"+outname3+".png");
233  c3->SaveAs(outDir+"/"+outname3+".pdf");
234  }
235 }
236 
237 //----------------------------------------------
238 // Hadronic Energy
239 //----------------------------------------------
240 void HadEFit(ana::Spectrum &plot1, TString zaxisname, TString outname1,
241  TString outname2,TString outname3, Bool_t printme=false){
242  double pot = 8.09e20;
243  TH2* hist = plot1.ToTH2(pot);
244  hist->GetXaxis()->SetTitle("Visible Hadronic Energy (GeV)");
245  hist->GetYaxis()->SetTitle("True Neutrino E - Reco Muon E (GeV)");
246  hist->GetZaxis()->SetTitle(zaxisname);
247  hist->GetXaxis()->CenterTitle();
248  hist->GetYaxis()->CenterTitle();
249  hist->GetZaxis()->CenterTitle();
250  static const int nn = 150; //Number of points on graph - should be equal to or smaller than number of bins in x
251  float xxx [nn];
252  float exxl[nn];
253  float exxh[nn];
254  float yyx [nn];
255  float eyxl[nn];
256  float eyxh[nn];
257 
258  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
259 
260  // make projection in Y view
261  TH1D* py = hist->ProjectionY("py",iw,iw,"");
262  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
263  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
264 
265  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
266 
267  double totalEntries = py->GetEntries();
268  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
269 
270  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
271  //Find the location of width at half max
272  int yyhighbin = maxbin+1;
273  int yylowbin = maxbin-1;
274  // find upper bin of half max
275  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
276  yyhighbin = ihigh;
277  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
278  }
279  // find lower bin of half max
280  for(int ilow = maxbin-1; ilow > 0; --ilow){
281  yylowbin = ilow;
282  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
283  }
284  double multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
285  if (iw>100){multiplier = 10.0;}
286  if (iw<3){multiplier = 0.5;}
287  if ((iw<20)&&(iw>2)){multiplier = 1.0;}
288  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
289  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
290 
291  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
292  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
293 
294  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
295  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
296  py->Fit("susieGaus","OQ","",lowSide,highSide); //O and Q mean: do not plot result, quiet mode
297  double mean = f1->GetParameter(1);
298  double error = f1->GetParError(1);
299  double chi = f1->GetChisquare();
300  double ndf = f1->GetNDF();
301 
302  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
303 
304  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
305  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
306  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
307  //Setting the error on the mean reported by the fit as error in y
308  eyxl[iw-1] = error;
309  eyxh[iw-1] = error;
310 
311  }//End of loop over bins in 2D histogram
312 
313  //Making graph from points
314  TGraphAsymmErrors* Profile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
315 
316  Profile->GetXaxis()->SetNoExponent(kTRUE);
317  Profile->GetXaxis()->SetTitle("Visible Hadronic E (GeV)");
318  Profile->GetYaxis()->SetTitle("True Neutrino E - Reco Muon E (GeV)");
319  Profile->GetXaxis()->CenterTitle();
320  Profile->GetYaxis()->CenterTitle();
321  Profile->SetMarkerStyle(6);
322  Profile->GetXaxis()->SetLimits(0,2.0);
323  Profile->GetHistogram()->SetMaximum(3);
324  Profile->GetHistogram()->SetMinimum(0);
325  Profile->SetTitle("");
326  // Profile->Write();
327  // f1.Close();
328 
329  //Overlays 2d plot and graph
330  TCanvas* c1 = new TCanvas("c1","c1");
331  Profile->SetMarkerColor(18);
332  Profile->SetLineColor(18);
333  hist->Draw("colz");
334  Profile->Draw("P");
335  //gPad->Update();
336  Simulation();
337  c1->SetLeftMargin(0.15);
338  c1->SetRightMargin(0.15);
339 
340  if (printme){
341  TString outDir = "plots";
342  c1->SaveAs(outDir+"/"+outname1+".png");
343  c1->SaveAs(outDir+"/"+outname1+".pdf");
344  //c1->SaveAs(outDir+"/"+outname+".root");
345  }
346  TCanvas* c2 = new TCanvas("c2","c2");
347  c2->SetLogz();
348  Profile->SetMarkerColor(13);
349  Profile->SetLineColor(13);
350  hist->Draw("colz");
351  hist->GetYaxis()->SetRangeUser(0,3);
352  Profile->Draw("P");
353  Simulation();
354  //gPad->Update();
355  c2->SetRightMargin(0.15);
356  c2->SetLeftMargin(0.15);
357 
358  if (printme){
359  TString outDir = "plots";
360  c2->SaveAs(outDir+"/"+outname2+".png");
361  c2->SaveAs(outDir+"/"+outname2+".pdf");
362  //c1->SaveAs(outDir+"/"+outname+".root");
363  }
364 
365  TCanvas* c3 = new TCanvas("c3","c3");
366  Profile->SetMarkerColor(1);
367  Profile->SetLineColor(1);
368  Profile->Draw("AP");
369  c3->SetLeftMargin(0.15);
370  Profile->Fit("pol7", "","",0.001, 2.0);
371 
372  if (printme){
373  TString outDir = "plots";
374  c3->SaveAs(outDir+"/"+outname3+".png");
375  c3->SaveAs(outDir+"/"+outname3+".pdf");
376  }
377  } // end of void
378 
void Simulation()
Definition: tools.h:16
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
std::string outDir
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
c2
Definition: demo5.py:33
void CenterTitles(TH1 *h)
#define pot
Float_t f1
void HadEFit(ana::Spectrum &plot1, TString zaxisname, TString outname1, TString outname2, TString outname3, Bool_t printme=false)
Definition: PlotStyle.h:240
void MuonFit(ana::Spectrum &plot1, TString zaxisname, TString outname1, TString outname2, TString outname3, Bool_t act=false, Bool_t actandcat=false, Bool_t cat=false, Bool_t printme=false)
Definition: PlotStyle.h:57
c1
Definition: demo5.py:24