MakingProfileHadQE.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 = 63; //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("2");
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("hadronqeMother.pdf");
50  /*
51 
52  //Testing a thing
53  TH1D* py = hist->ProjectionY("py",10,10,"");
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)*25+maxbin;
73  int veryLowBin = maxbin - 25*(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-25*(modex-py->GetBinCenter(yylowbin));
86  double highSide = 25*(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("testThing10.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 > 60) {multiplier = 10;}
135  if (iw > 65) {multiplier = 25;}
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,2.0);
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("hadronqeMotherAndGraph.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("hadronqeMotherLogzAndGraph.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("hadronqeGraph.pdf");
220 
221  //Drawing fit over the 2d plot and graph
222 
223  //Old Tuning Values
224  //double stitch1 = 0.0593168; // GeV
225  //double stitch2 = 0.208919; // GeV
226  //double stitch3 = 1.05000; // GeV
227  //double offset = 0.0534636; // GeV
228  //double slope1 = 0.631518; // unitless
229  //double slope2 = 1.39952; // unitless
230  //double slope3 = 1.73600; // unitless
231  //double slope4 = 2.19266; // unitless
232 
233  //New Tuning Values
234  //double stitch1 = 0.0597385; // GeV
235  //double stitch2 = 0.138985; // GeV
236  //double stitch3 = 0.800025; // GeV
237  //double offset = 0.0514985; // GeV
238  //double slope1 = 0.620593; // unitless
239  //double slope2 = 1.47160; // unitless
240  //double slope3 = 1.80880; // unitless
241  //double slope4 = 2.06219; // unitless
242 
243  //New Tuning Values Rounded for Errors
244  double stitch1 = 0.0597; // GeV
245  double stitch2 = 0.139; // GeV
246  double stitch3 = 0.800; // GeV
247  double offset = 0.0515; // GeV
248  double slope1 = 0.621; // unitless
249  double slope2 = 1.47; // unitless
250  double slope3 = 1.81; // unitless
251  double slope4 = 2.06; // 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,2.0,slope4*2.0+(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  gPad->Update();
288  gPad->SetRightMargin(0.13);
289  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
290  //palette1->SetBBoxCenterX(625);
291  gPad->Modified();
292  gPad->Update();
293  c4->Print("fittingLineOnGraphhadqe.pdf");
294 
295  //2d plot with fitting line
296  TCanvas* c5 = new TCanvas("c5","c5");
297  hist->Draw("colz");
298  fitLine1.Draw();
299  fitLine2.Draw();
300  fitLine3.Draw();
301  fitLine4.Draw();
302  stitchLine1.Draw();
303  stitchLine2.Draw();
304  stitchLine3.Draw();
305  gPad->Update();
306  gPad->SetRightMargin(0.13);
307  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
308  //palette1->SetBBoxCenterX(625);
309  gPad->Modified();
310  gPad->Update();
311  c5->Print("fittingLineOnMotherhadqe.pdf");
312 
313  //2d plot with fitting line on logz
314  TCanvas* c8 = new TCanvas("c8","c8");
315  c8->SetLogz();
316  hist->Draw("colz");
317  fitLine1.Draw();
318  fitLine2.Draw();
319  fitLine3.Draw();
320  fitLine4.Draw();
321  stitchLine1.Draw();
322  stitchLine2.Draw();
323  stitchLine3.Draw();
324  gPad->Update();
325  gPad->SetRightMargin(0.13);
326  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
327  //palette1->SetBBoxCenterX(625);
328  gPad->Modified();
329  gPad->Update();
330  c8->Print("fittingLineOnMotherLogzhadqe.pdf");
331 
332  //Plot of chi squared per ndf
333  TCanvas* c6 = new TCanvas("c6","c6");
334  chiValues->Draw();
335  c6->Print("graphChiSquaredperNDFhadqe.pdf");
336 
337  TFile f("hadronqeProfile.root","RECREATE");
338  susieProfile->Write("1");
339  chiValues->Write();
340 
341 
342 
343 }
344 
345 
346 
void MakingProfileHadQE()
c2
Definition: demo5.py:33
Float_t f1
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24