MakingProfileHadNonQE.C
Go to the documentation of this file.
1 //Run using: root -b -q -l MakingProfile.C
2 
3 //This script is where I make the graph from 2d plots that I shall fit. I also then
4 //make nice versions of the graph and 2d plot with fit lines on it.
5 
6 #include <iostream>
7 #include <fstream>
8 #include <iomanip>
9 #include <sstream>
10 #include <string>
11 #include <cstring>
12 #include "TFile.h"
13 #include "TH2D.h"
14 #include "TObject.h"
15 #include "TArray.h"
16 #include "TArrayF.h"
17 #include "TGraphAsymmErrors.h"
18 
19 
21 {
22  std::cout<<"Hello Susie!"<<std::endl;
23 
24  static const int nn = 130; //Number of points on graph - should be equal to or smaller than number of bins in x
25  float xxx[nn];
26  float exxl[nn];
27  float exxh[nn];
28  float yyx[nn];
29  float eyxl[nn];
30  float eyxh[nn];
31 
32  gStyle->SetOptStat(0);
33  TFile* file = new TFile("./2DPlotsForFitting.root","READ");
34 
35  TH1D* chiValues = new TH1D("chiValue",";Chi squared per NDF;Fits",100, 0.0, 10.0);
36 
37  TH2D* hist = file->Get("3");
38 
39  TCanvas* c7 = new TCanvas("c7","c7");
40  hist->GetXaxis()->CenterTitle();
41  hist->GetYaxis()->CenterTitle();
42  hist->Draw("colz");
43  gPad->Update();
44  gPad->SetRightMargin(0.13);
45  TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
46  palette1->SetBBoxCenterX(625);
47  gPad->Modified();
48  gPad->Update();
49  c7->Print("hadronnonqeMother.pdf");
50 
51  /*
52  //Testing a thing
53  TH1D* py = hist->ProjectionY("py",2,2,"");
54  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
55  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
56  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
57  std::cout<<" max bin is: "<<maxbin<<" maxcontent: "<<maxcont<<" modex " <<modex<<std::endl;
58 
59  int yyhighbin = maxbin+1;
60  int yylowbin = maxbin-1;
61  // find upper bin of half max
62  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
63  yyhighbin = ihigh;
64  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
65  }
66  // find lower bin of half max
67  for(int ilow = maxbin-1; ilow > 0; --ilow){
68  yylowbin = ilow;
69  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
70  }
71 
72  int veryHighBin = (yyhighbin - maxbin)*5+maxbin;
73  int veryLowBin = maxbin - 5*(maxbin-yylowbin);
74  std::cout<<"very high bin "<<veryHighBin<<" very low bin "<<veryLowBin<<std::endl;
75 
76  //Count number of entries in truncated histogram
77  int numEntries = 0;
78  for(int iBin = veryLowBin; iBin <= veryHighBin; ++iBin){
79  int entries = py->GetBinContent(iBin);
80  std::cout<<" i am bin "<<iBin<<" and I have "<<entries<<" entries"<<std::endl;
81  numEntries = numEntries + entries;
82  }
83  std::cout<<" number of entries total is "<<numEntries<<std::endl;
84 
85  double lowSide = modex-5*(modex-py->GetBinCenter(yylowbin));
86  double highSide = 5*(py->GetBinCenter(yyhighbin)-modex)+modex;
87  std::cout<<"low side "<<lowSide<<" high side "<<highSide<<std::endl;
88 
89  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
90 
91  f1->SetParameters(50,0.3,0.05);
92 
93  TCanvas* c6 = new TCanvas("c6","c6");
94  py->Fit("susieGaus","","",lowSide,highSide);
95  c6->Print("testThing.pdf");
96 
97  double mean = f1->GetParameter(1);
98  double sigma = f1->GetParameter(2);
99  double error = sigma / sqrt(numEntries);
100  double chrisError = f1->GetParError(1);
101 
102  std::cout<<" The mean is "<<mean<<" and the sigma is "<<sigma<<" and error "<<error<<" and chris error is "<<chrisError<<std::endl;
103  */
104 
105  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
106 
107  // make projection in Y view
108  TH1D* py = hist->ProjectionY("py",iw,iw,"");
109  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
110  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
111 
112  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
113 
114  double totalEntries = py->GetEntries();
115  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
116 
117  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
118 
119  //Find the location of width at half max
120  int yyhighbin = maxbin+1;
121  int yylowbin = maxbin-1;
122  // find upper bin of half max
123  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
124  yyhighbin = ihigh;
125  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
126  }
127  // find lower bin of half max
128  for(int ilow = maxbin-1; ilow > 0; --ilow){
129  yylowbin = ilow;
130  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
131  }
132 
133  int multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
134  if (iw>124){multiplier = 10.0;}
135  if (iw<3){multiplier = 5.0;}
136  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
137  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
138 
139  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
140  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
141 
142  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
143 
144  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
145  py->Fit("susieGaus","OQ","",lowSide,highSide); //O and Q mean: do not plot result, quiet mode
146  double mean = f1->GetParameter(1);
147  double error = f1->GetParError(1);
148  double chi = f1->GetChisquare();
149  double ndf = f1->GetNDF();
150 
151  if (ndf > 0) {chiValues->Fill(chi/ndf);}
152  else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
153 
154  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
155 
156  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
157  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
158  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
159  //Setting the error on the mean reported by the fit as error in y
160  eyxl[iw-1] = error;
161  eyxh[iw-1] = error;
162 
163  }//End of loop over bins in 2D histogram
164 
165  //Making graph from points
166  TGraphAsymmErrors* susieProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
167 
168  susieProfile->GetXaxis()->SetNoExponent(kTRUE);
169  susieProfile->GetXaxis()->SetTitle("Visible Hadronic E (GeV)");
170  susieProfile->GetYaxis()->SetTitle("True Neutrino E - Reco Muon E (GeV)");
171  susieProfile->GetXaxis()->CenterTitle();
172  susieProfile->GetYaxis()->CenterTitle();
173  susieProfile->SetMarkerStyle(6);
174  susieProfile->GetXaxis()->SetLimits(0,1.98);
175  susieProfile->GetHistogram()->SetMaximum(5);
176  susieProfile->GetHistogram()->SetMinimum(0);
177  susieProfile->SetTitle("");
178 
179  //Overlays 2d plot and graph
180  TCanvas* c1 = new TCanvas("c1","c1");
181  susieProfile->SetMarkerColor(18);
182  susieProfile->SetLineColor(18);
183  hist->Draw("colz");
184  susieProfile->Draw("P");
185  gPad->Update();
186  gPad->SetRightMargin(0.13);
187  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
188  //palette1->SetBBoxCenterX(625);
189  gPad->Modified();
190  gPad->Update();
191  c1->Print("hadronnonqeMotherAndGraph.pdf");
192 
193  //Overlays 2d loqz plot and graph
194  TCanvas* c2 = new TCanvas("c2","c2");
195  c2->SetLogz();
196  susieProfile->SetMarkerColor(13);
197  susieProfile->SetLineColor(13);
198  hist->Draw("colz");
199  susieProfile->Draw("P");
200  gPad->Update();
201  gPad->SetRightMargin(0.13);
202  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
203  //palette1->SetBBoxCenterX(625);
204  gPad->Modified();
205  gPad->Update();
206  c2->Print("hadronnonqeMotherLogzAndGraph.pdf");
207 
208  //Just the graph
209  TCanvas* c3 = new TCanvas("c3","c3");
210  susieProfile->SetMarkerColor(1);
211  susieProfile->SetLineColor(1);
212  susieProfile->Draw("AP");
213  gPad->Update();
214  gPad->SetRightMargin(0.13);
215  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
216  //palette1->SetBBoxCenterX(625);
217  gPad->Modified();
218  gPad->Update();
219  c3->Print("hadronnonqeGraph.pdf");
220 
221  //Drawing fit over the 2d plot and graph
222 
223  //Old Tuning Values
224  //double stitch1 = 0.203006; // GeV
225  //double stitch2 = 0.275000; // GeV
226  //double stitch3 = 1.15544; // GeV
227  //double offset = 0.232125; // GeV
228  //double slope1 = 0.954274; // unitless
229  //double slope2 = 1.28840; // unitless
230  //double slope3 = 1.66420; // unitless
231  //double slope4 = 1.90552; // unitless
232 
233  //New Tuning Values
234  //double stitch1 = 0.131748; // GeV
235  //double stitch2 = 0.222564; // GeV
236  //double stitch3 = 0.963558; // GeV
237  //double offset = 0.241399; // GeV
238  //double slope1 = 0.943280; // unitless
239  //double slope2 = 1.20946; // unitless
240  //double slope3 = 1.78345; // unitless
241  //double slope4 = 2.19532; // unitless
242 
243  //New Tuning Values with Error Rounding
244  double stitch1 = 0.132; // GeV
245  double stitch2 = 0.223; // GeV
246  double stitch3 = 0.964; // GeV
247  double offset = 0.241; // GeV
248  double slope1 = 0.94; // unitless
249  double slope2 = 1.21; // unitless
250  double slope3 = 1.78; // unitless
251  double slope4 = 2.20; // unitless
252 
253  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
254  fitLine1.SetLineColor(2);
255 
256  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,stitch2,slope2*stitch2+(slope1-slope2)*stitch1+offset);
257  fitLine2.SetLineColor(2);
258 
259  fitLine3 = TLine(stitch2,slope3*stitch2+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset,stitch3,slope3*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset);
260  fitLine3.SetLineColor(2);
261 
262  fitLine4 = TLine(stitch3,slope4*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+(slope3-slope4)*stitch3+offset,1.98,slope4*1.98+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+(slope3-slope4)*stitch3+offset);
263  fitLine4.SetLineColor(2);
264 
265  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
266  stitchLine1.SetLineColor(2);
267  stitchLine1.SetLineStyle(7);
268 
269  stitchLine2 = TLine(stitch2, 0, stitch2,5.0);
270  stitchLine2.SetLineColor(2);
271  stitchLine2.SetLineStyle(7);
272 
273  stitchLine3 = TLine(stitch3, 0, stitch3,5.0);
274  stitchLine3.SetLineColor(2);
275  stitchLine3.SetLineStyle(7);
276 
277  //Graph with fitting line
278  TCanvas* c4 = new TCanvas("c4","c4");
279  susieProfile->Draw("AP");
280  fitLine1.Draw();
281  fitLine2.Draw();
282  fitLine3.Draw();
283  fitLine4.Draw();
284  stitchLine1.Draw();
285  stitchLine2.Draw();
286  stitchLine3.Draw();
287  c4->Print("fittingLineOnGraphhadnonqe.pdf");
288 
289  //2d plot with fitting line
290  TCanvas* c5 = new TCanvas("c5","c5");
291  hist->Draw("colz");
292  fitLine1.Draw();
293  fitLine2.Draw();
294  fitLine3.Draw();
295  fitLine4.Draw();
296  stitchLine1.Draw();
297  stitchLine2.Draw();
298  stitchLine3.Draw();
299  gPad->Update();
300  gPad->SetRightMargin(0.13);
301  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
302  //palette1->SetBBoxCenterX(625);
303  gPad->Modified();
304  gPad->Update();
305  c5->Print("fittingLineOnMotherhadnonqe.pdf");
306 
307  //Plot of chi squared per ndf
308  TCanvas* c6 = new TCanvas("c6","c6");
309  chiValues->Draw();
310  c6->Print("graphChiSquaredperNDFhadnonqe.pdf");
311 
312  //2d plot with fitting line on logz
313  TCanvas* c8 = new TCanvas("c8","c8");
314  c8->SetLogz();
315  hist->Draw("colz");
316  fitLine1.Draw();
317  fitLine2.Draw();
318  fitLine3.Draw();
319  fitLine4.Draw();
320  stitchLine1.Draw();
321  stitchLine2.Draw();
322  stitchLine3.Draw();
323  gPad->Update();
324  gPad->SetRightMargin(0.13);
325  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
326  //palette1->SetBBoxCenterX(625);
327  gPad->Modified();
328  gPad->Update();
329  c8->Print("fittingLineOnMotherLogzhadnonqe.pdf");
330 
331  TFile f("hadronnonqeProfile.root","RECREATE");
332  susieProfile->Write("1");
333  chiValues->Write();
334 
335 
336 
337 }
338 
339 
340 
c2
Definition: demo5.py:33
void MakingProfileHadNonQE()
Float_t f1
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24