MakingProfileHadNonQEND.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 = 115; //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("7");
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",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  std::cout<<"high bin "<<yyhighbin<<" very low bin "<<yylowbin<<std::endl;
72  double multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
73  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
74  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
75  std::cout<<"very high bin "<<veryHighBin<<" very low bin "<<veryLowBin<<std::endl;
76 
77  //Count number of entries in truncated histogram
78  int numEntries = 0;
79  for(int iBin = veryLowBin; iBin <= veryHighBin; ++iBin){
80  int entries = py->GetBinContent(iBin);
81  std::cout<<" i am bin "<<iBin<<" and I have "<<entries<<" entries"<<std::endl;
82  numEntries = numEntries + entries;
83  }
84  std::cout<<" number of entries total is "<<numEntries<<std::endl;
85 
86  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
87  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
88  std::cout<<"low side "<<lowSide<<" high side "<<highSide<<std::endl;
89 
90  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
91 
92  f1->SetParameters(50,0.3,0.05);
93 
94  TCanvas* c6 = new TCanvas("c6","c6");
95  py->Fit("susieGaus","","",lowSide,highSide);
96  c6->Print("testThing.pdf");
97 
98  double mean = f1->GetParameter(1);
99  double sigma = f1->GetParameter(2);
100  double error = sigma / sqrt(numEntries);
101  double chrisError = f1->GetParError(1);
102 
103  std::cout<<" The mean is "<<mean<<" and the sigma is "<<sigma<<" and error "<<error<<" and chris error is "<<chrisError<<std::endl;
104  */
105 
106  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
107 
108  // make projection in Y view
109  TH1D* py = hist->ProjectionY("py",iw,iw,"");
110  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
111  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
112 
113  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
114 
115  double totalEntries = py->GetEntries();
116  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
117 
118  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
119 
120  //Find the location of width at half max
121  int yyhighbin = maxbin+1;
122  int yylowbin = maxbin-1;
123  // find upper bin of half max
124  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
125  yyhighbin = ihigh;
126  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
127  }
128  // find lower bin of half max
129  for(int ilow = maxbin-1; ilow > 0; --ilow){
130  yylowbin = ilow;
131  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
132  }
133 
134  double multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
135  if (iw>100){multiplier = 10.0;}
136  if (iw<3){multiplier = 0.5;}
137  if ((iw<20)&&(iw>2)){multiplier = 1.0;}
138  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
139  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
140 
141  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
142  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
143 
144  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
145 
146  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
147  py->Fit("susieGaus","OQ","",lowSide,highSide); //O and Q mean: do not plot result, quiet mode
148  double mean = f1->GetParameter(1);
149  double error = f1->GetParError(1);
150  double chi = f1->GetChisquare();
151  double ndf = f1->GetNDF();
152 
153  if (ndf > 0) {chiValues->Fill(chi/ndf);}
154  else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
155 
156  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
157 
158  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
159  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
160  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
161  //Setting the error on the mean reported by the fit as error in y
162  eyxl[iw-1] = error;
163  eyxh[iw-1] = error;
164 
165  }//End of loop over bins in 2D histogram
166 
167  //Making graph from points
168  TGraphAsymmErrors* susieProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
169 
170  susieProfile->GetXaxis()->SetNoExponent(kTRUE);
171  susieProfile->GetXaxis()->SetTitle("Visible Hadronic E (GeV)");
172  susieProfile->GetYaxis()->SetTitle("True Neutrino E - Reco Muon E (GeV)");
173  susieProfile->GetXaxis()->CenterTitle();
174  susieProfile->GetYaxis()->CenterTitle();
175  susieProfile->SetMarkerStyle(6);
176  susieProfile->GetXaxis()->SetLimits(0,2.0);
177  susieProfile->GetHistogram()->SetMaximum(5);
178  susieProfile->GetHistogram()->SetMinimum(0);
179  susieProfile->SetTitle("");
180 
181  //Overlays 2d plot and graph
182  TCanvas* c1 = new TCanvas("c1","c1");
183  susieProfile->SetMarkerColor(18);
184  susieProfile->SetLineColor(18);
185  hist->Draw("colz");
186  susieProfile->Draw("P");
187  gPad->Update();
188  gPad->SetRightMargin(0.13);
189  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
190  //palette1->SetBBoxCenterX(625);
191  gPad->Modified();
192  gPad->Update();
193  c1->Print("hadronnonqeMotherAndGraph.pdf");
194 
195  //Overlays 2d loqz plot and graph
196  TCanvas* c2 = new TCanvas("c2","c2");
197  c2->SetLogz();
198  susieProfile->SetMarkerColor(13);
199  susieProfile->SetLineColor(13);
200  hist->Draw("colz");
201  susieProfile->Draw("P");
202  gPad->Update();
203  gPad->SetRightMargin(0.13);
204  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
205  //palette1->SetBBoxCenterX(625);
206  gPad->Modified();
207  gPad->Update();
208  c2->Print("hadronnonqeMotherLogzAndGraph.pdf");
209 
210  //Just the graph
211  TCanvas* c3 = new TCanvas("c3","c3");
212  susieProfile->SetMarkerColor(1);
213  susieProfile->SetLineColor(1);
214  susieProfile->Draw("AP");
215  gPad->Update();
216  gPad->SetRightMargin(0.13);
217  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
218  //palette1->SetBBoxCenterX(625);
219  gPad->Modified();
220  gPad->Update();
221  c3->Print("hadronnonqeGraph.pdf");
222 
223  //Drawing fit over the 2d plot and graph
224 
225  //Old Tune Values
226  //double stitch1 = 0.225007; // GeV
227  //double offset = 0.202837; // GeV
228  //double slope1 = 1.21606; // unitless
229  //double slope2 = 1.56778; // unitless
230 
231  //New Tune Values
232  //double stitch1 = 0.168785; // GeV
233  //double offset = 0.253713; // GeV
234  //double slope1 = 1.08492; // unitless
235  //double slope2 = 1.88869; // unitless
236 
237  //New Tune Values with Error Rounding
238  double stitch1 = 0.169; // GeV
239  double offset = 0.254; // GeV
240  double slope1 = 1.080; // unitless
241  double slope2 = 1.889; // unitless
242 
243  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
244  fitLine1.SetLineColor(2);
245 
246  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,2.0,slope2*2.0+(slope1-slope2)*stitch1+offset);
247  fitLine2.SetLineColor(2);
248 
249 
250  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
251  stitchLine1.SetLineColor(2);
252  stitchLine1.SetLineStyle(7);
253 
254 
255  //Graph with fitting line
256  TCanvas* c4 = new TCanvas("c4","c4");
257  susieProfile->Draw("AP");
258  fitLine1.Draw();
259  fitLine2.Draw();
260  stitchLine1.Draw();
261  gPad->Update();
262  gPad->SetRightMargin(0.13);
263  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
264  //palette1->SetBBoxCenterX(625);
265  gPad->Modified();
266  gPad->Update();
267  c4->Print("fittingLineOnGraphhadnonqe.pdf");
268 
269  //2d plot with fitting line
270  TCanvas* c5 = new TCanvas("c5","c5");
271  hist->Draw("colz");
272  fitLine1.Draw();
273  fitLine2.Draw();
274  stitchLine1.Draw();
275  gPad->Update();
276  gPad->SetRightMargin(0.13);
277  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
278  //palette1->SetBBoxCenterX(625);
279  gPad->Modified();
280  gPad->Update();
281  c5->Print("fittingLineOnMotherhadnonqe.pdf");
282 
283  //Plot of chi squared per ndf
284  TCanvas* c6 = new TCanvas("c6","c6");
285  chiValues->Draw();
286  c6->Print("graphChiSquaredperNDFhadnonqe.pdf");
287 
288  //2d plot with fitting line on logz
289  TCanvas* c8 = new TCanvas("c8","c8");
290  c8->SetLogz();
291  hist->Draw("colz");
292  fitLine1.Draw();
293  fitLine2.Draw();
294  stitchLine1.Draw();
295  gPad->Update();
296  gPad->SetRightMargin(0.13);
297  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
298  //palette1->SetBBoxCenterX(625);
299  gPad->Modified();
300  gPad->Update();
301  c8->Print("fittingLineOnMotherLogzhadnonqe.pdf");
302 
303  TFile f("hadronnonqeProfile.root","RECREATE");
304  susieProfile->Write("1");
305  chiValues->Write();
306 
307 
308 
309 }
310 
311 
312 
void MakingProfileHadNonQEND()
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