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