MakingProfileHadCC.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<<"Hiya!"<<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("19");
38 
39  TCanvas* c7 = new TCanvas("c7","c7");
40  hist->Draw("colz");
41  c7->Print("hadronccMother.pdf");
42 
43  int grPoints = 0;
44 
45  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
46 
47  // make projection in Y view
48  TH1D* py = hist->ProjectionY("py",iw,iw,"");
49  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
50  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
51 
52  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
53 
54  double totalEntries = py->GetEntries();
55  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
56 
57  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
58 
59  //Find the location of width at half max
60  int yyhighbin = maxbin+1;
61  int yylowbin = maxbin-1;
62  // find upper bin of half max
63  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
64  yyhighbin = ihigh;
65  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
66  }
67  // find lower bin of half max
68  for(int ilow = maxbin-1; ilow > 0; --ilow){
69  yylowbin = ilow;
70  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
71  }
72 
73  if(abs(yylowbin-maxbin) < 3 || abs(yylowbin-maxbin) < 3){continue; } // Added by me: kills bad points
74 
75  int multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
76  if (iw>124){multiplier = 10.0;}
77  if (iw<3){multiplier = 5.0;}
78  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
79  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
80 
81  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
82  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
83 
84  xxx[grPoints] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
85 
86  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
87  py->Fit("susieGaus","OQ","",lowSide,highSide); //O and Q mean: do not plot result, quiet mode
88  double mean = f1->GetParameter(1);
89  double error = f1->GetParError(1);
90  double chi = f1->GetChisquare();
91  double ndf = f1->GetNDF();
92 
93  if(mean < 0){continue; } // Another security check
94  if(mean > 10){continue; } // And yet another
95 
96  if (ndf > 0) {chiValues->Fill(chi/ndf);}
97  else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
98 
99  yyx[grPoints] = mean; //Set mean of gaussian fit as y location of graph point
100 
101  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
102  exxl[grPoints] = 0; //effVwx->GetBinWidth(iw)/2.0;
103  exxh[grPoints] = 0; //effVwx->GetBinWidth(iw)/2.0;
104  //Setting the error on the mean reported by the fit as error in y
105  eyxl[grPoints] = error;
106  eyxh[grPoints] = error;
107  grPoints++;
108 
109  }//End of loop over bins in 2D histogram
110 
111  //Making graph from points
112  TGraphAsymmErrors* susieProfile = new TGraphAsymmErrors(grPoints,xxx,yyx,exxl,exxh,eyxl,eyxh);
113 
114  susieProfile->GetXaxis()->SetNoExponent(kTRUE);
115  susieProfile->GetXaxis()->SetTitle("Visible Hadronic E (GeV)");
116  susieProfile->GetYaxis()->SetTitle("True Neutrino E - Reco Muon E (GeV)");
117  susieProfile->SetMarkerStyle(6);
118  susieProfile->GetXaxis()->SetLimits(0,2.0);
119  susieProfile->GetHistogram()->SetMaximum(5);
120  susieProfile->GetHistogram()->SetMinimum(0);
121  susieProfile->SetTitle("");
122 
123  //Overlays 2d plot and graph
124  TCanvas* c1 = new TCanvas("c1","c1");
125  susieProfile->SetMarkerColor(18);
126  susieProfile->SetLineColor(18);
127  hist->Draw("colz");
128  susieProfile->Draw("P");
129  c1->Print("hadronccMotherAndGraph.pdf");
130 
131  //Overlays 2d loqz plot and graph
132  TCanvas* c2 = new TCanvas("c2","c2");
133  c2->SetLogz();
134  susieProfile->SetMarkerColor(13);
135  susieProfile->SetLineColor(13);
136  hist->Draw("colz");
137  susieProfile->Draw("P");
138  c2->Print("hadronccMotherLogzAndGraph.pdf");
139 
140  //Just the graph
141  TCanvas* c3 = new TCanvas("c3","c3");
142  susieProfile->SetMarkerColor(1);
143  susieProfile->SetLineColor(1);
144  susieProfile->Draw("AP");
145  c3->Print("hadronccGraph.pdf");
146 
147  //Drawing fit over the 2d plot and graph
148 
149  // This is for a 4-spline, moving to 3:
150 
151  // double stitch1 = 0.203006; // GeV
152  // double stitch2 = 0.275000; // GeV
153  // double stitch3 = 1.15544; // GeV
154  // double offset = 0.232125; // GeV
155  // double slope1 = 0.954274; // unitless
156  // double slope2 = 1.28840; // unitless
157  // double slope3 = 1.66420; // unitless
158  // double slope4 = 1.90552; // unitless
159 
160  // fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
161  // fitLine1.SetLineColor(2);
162 
163  // fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,stitch2,slope2*stitch2+(slope1-slope2)*stitch1+offset);
164  // fitLine2.SetLineColor(2);
165 
166  // fitLine3 = TLine(stitch2,slope3*stitch2+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset,stitch3,slope3*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset);
167  // fitLine3.SetLineColor(2);
168 
169  // 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);
170  // fitLine4.SetLineColor(2);
171 
172  // stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
173  // stitchLine1.SetLineColor(2);
174  // stitchLine1.SetLineStyle(7);
175 
176  // stitchLine2 = TLine(stitch2, 0, stitch2,5.0);
177  // stitchLine2.SetLineColor(2);
178  // stitchLine2.SetLineStyle(7);
179 
180  // stitchLine3 = TLine(stitch3, 0, stitch3,5.0);
181  // stitchLine3.SetLineColor(2);
182  // stitchLine3.SetLineStyle(7);
183 
184  double stitch1 = 3.75525e-02; // GeV
185  double stitch2 = 2.05000e-01; // GeV
186  double stitch3 = 1.22411e+00; // GeV
187  double offset = 2.71346e-01; // GeV
188  double slope1 = 1.26967e+00; // unitless
189  double slope2 = 9.31211e-01; // unitless
190  double slope3 = 1.56865e+00; // unitless
191  double slope4 = 2.30327e+00; // unitless
192 
193  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
194  fitLine1.SetLineColor(2);
195 
196  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,stitch2,slope2*stitch2+(slope1-slope2)*stitch1+offset);
197  fitLine2.SetLineColor(2);
198 
199  fitLine3 = TLine(stitch2,slope3*stitch2+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset,stitch3,slope3*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset);
200  fitLine3.SetLineColor(2);
201 
202  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);
203  fitLine4.SetLineColor(2);
204 
205  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
206  stitchLine1.SetLineColor(2);
207  stitchLine1.SetLineStyle(7);
208 
209  stitchLine2 = TLine(stitch2, 0, stitch2,5.0);
210  stitchLine2.SetLineColor(2);
211  stitchLine2.SetLineStyle(7);
212 
213  stitchLine3 = TLine(stitch3, 0, stitch3,5.0);
214  stitchLine3.SetLineColor(2);
215  stitchLine3.SetLineStyle(7);
216 
217  //Graph with fitting line
218  TCanvas* c4 = new TCanvas("c4","c4");
219  susieProfile->Draw("AP");
220  fitLine1.Draw();
221  fitLine2.Draw();
222  fitLine3.Draw();
223  fitLine4.Draw();
224  stitchLine1.Draw();
225  stitchLine2.Draw();
226  stitchLine3.Draw();
227  c4->Print("fittingLineOnGraphhadcc.pdf");
228 
229  //2d plot with fitting line
230  TCanvas* c5 = new TCanvas("c5","c5");
231  hist->Draw("colz");
232  fitLine1.Draw();
233  fitLine2.Draw();
234  fitLine3.Draw();
235  fitLine4.Draw();
236  stitchLine1.Draw();
237  stitchLine2.Draw();
238  stitchLine3.Draw();
239  c5->Print("fittingLineOnMotherhadcc.pdf");
240 
241  //Plot of chi squared per ndf
242  TCanvas* c6 = new TCanvas("c6","c6");
243  chiValues->Draw();
244  c6->Print("graphChiSquaredperNDFhadcc.pdf");
245 
246  //2d plot with fitting line on logz
247  TCanvas* c8 = new TCanvas("c8","c8");
248  c8->SetLogz();
249  hist->Draw("colz");
250  fitLine1.Draw();
251  fitLine2.Draw();
252  fitLine3.Draw();
253  fitLine4.Draw();
254  stitchLine1.Draw();
255  stitchLine2.Draw();
256  stitchLine3.Draw();
257  c8->Print("fittingLineOnMotherLogzhadcc.pdf");
258 
259  TFile f("hadronccProfile.root","RECREATE");
260  susieProfile->Write("1");
261  chiValues->Write();
262 
263 }
264 
265 
void abs(TH1 *hist)
c2
Definition: demo5.py:33
void MakingProfileHadCC()
Float_t f1
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24