MakingProfileND.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 = 150; //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",150, 0.0, 5.0);
36 
37  TH2D* hist = file->Get("1");
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("muonMother.pdf");
50  /*
51 
52  //Testing a thing
53  TH1D* py = hist->ProjectionY("py",130,130,"");
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  int multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
72  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
73  int veryLowBin = maxbin - multiplier*(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-multiplier*(modex-py->GetBinCenter(yylowbin));
86  double highSide = multiplier*(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  double chi = f1->GetChisquare();
105  double ndf = f1->GetNDF();
106 
107  std::cout<<" This chi per ndf is: "<<chi<<" / "<<ndf<<" = "<<chi/ndf<<std::endl;
108  */
109 
110  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
111 
112  // make projection in Y view
113  TH1D* py = hist->ProjectionY("py",iw,iw,"");
114  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
115  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
116 
117  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
118 
119  double totalEntries = py->GetEntries();
120  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
121 
122  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
123 
124  //Find the location of width at half max
125  int yyhighbin = maxbin+1;
126  int yylowbin = maxbin-1;
127  // find upper bin of half max
128  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
129  yyhighbin = ihigh;
130  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
131  }
132  // find lower bin of half max
133  for(int ilow = maxbin-1; ilow > 0; --ilow){
134  yylowbin = ilow;
135  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
136  }
137 
138  int multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
139  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
140  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
141 
142  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
143  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
144 
145  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
146 
147  TF1 *f1 = new TF1("susieGaus","gaus",0,5);
148  py->Fit("susieGaus","OQ","",lowSide,highSide); //O and Q mean do not plot result, quiet mode
149  double mean = f1->GetParameter(1);
150  double error = f1->GetParError(1);
151  double chi = f1->GetChisquare();
152  double ndf = f1->GetNDF();
153 
154  if (ndf > 0) {chiValues->Fill(chi/ndf);}
155  else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
156 
157  if ((chi/ndf) > 10) {std::cout<<" We have chi / ndf of : "<<chi/ndf<<" in the bin: "<<iw<<std::endl;}
158 
159  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
160 
161  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
162  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
163  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
164  //Setting the error on the mean reported by the fit as error in y
165  eyxl[iw-1] = error;
166  eyxh[iw-1] = error;
167 
168  }//End of loop over bins in 2D histogram
169 
170  //Making graph from points
171  TGraphAsymmErrors* susieProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
172 
173  susieProfile->GetXaxis()->SetNoExponent(kTRUE);
174  susieProfile->GetXaxis()->SetTitle("Reco Muon Track Length (cm)");
175  susieProfile->GetYaxis()->SetTitle("True Muon Energy (GeV)");
176  susieProfile->GetXaxis()->CenterTitle();
177  susieProfile->GetYaxis()->CenterTitle();
178  susieProfile->SetMarkerStyle(6);
179  susieProfile->GetXaxis()->SetLimits(0,1500);
180  susieProfile->GetHistogram()->SetMaximum(5);
181  susieProfile->GetHistogram()->SetMinimum(0);
182  susieProfile->SetTitle("");
183 
184  //Overlays 2d plot and graph
185  TCanvas* c1 = new TCanvas("c1","c1");
186  susieProfile->SetMarkerColor(18);
187  susieProfile->SetLineColor(18);
188  hist->Draw("colz");
189  susieProfile->Draw("P");
190  gPad->Update();
191  gPad->SetRightMargin(0.13);
192  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
193  //palette1->SetBBoxCenterX(625);
194  gPad->Modified();
195  gPad->Update();
196  c1->Print("muonMotherAndGraph.pdf");
197 
198  //Overlays 2d loqz plot and graph
199  TCanvas* c2 = new TCanvas("c2","c2");
200  c2->SetLogz();
201  susieProfile->SetMarkerColor(13);
202  susieProfile->SetLineColor(13);
203  hist->Draw("colz");
204  susieProfile->Draw("P");
205  gPad->Update();
206  gPad->SetRightMargin(0.13);
207  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
208  //palette1->SetBBoxCenterX(625);
209  gPad->Modified();
210  gPad->Update();
211  c2->Print("muonMotherLogzAndGraph.pdf");
212 
213  //Just the graph
214  TCanvas* c3 = new TCanvas("c3","c3");
215  susieProfile->SetMarkerColor(1);
216  susieProfile->SetLineColor(1);
217  susieProfile->Draw("AP");
218  gPad->Update();
219  gPad->SetRightMargin(0.13);
220  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
221  //palette1->SetBBoxCenterX(625);
222  gPad->Modified();
223  gPad->Update();
224  c3->Print("muonGraph.pdf");
225 
226  //Drawing fit over the 2d plot and graph
227 
228 
229  //New tuning values rounded for errors
230  double stitch1 = 334; // cm
231  double stitch2 = 539; // cm
232  double stitch3 = 1064; // cm
233  double offset = 0.1503; // GeV
234  double slope1 = 0.001910; // GeV/cm
235  double slope2 = 0.001987; // GeV/cm
236  double slope3 = 0.002039; // GeV/cm
237  double slope4 = 0.002159; // GeV/cm
238 
239  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
240  fitLine1.SetLineColor(2);
241 
242  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,stitch2,slope2*stitch2+(slope1-slope2)*stitch1+offset);
243  fitLine2.SetLineColor(2);
244 
245  fitLine3 = TLine(stitch2,slope3*stitch2+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset,stitch3,slope3*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset);
246  fitLine3.SetLineColor(2);
247 
248  fitLine4 = TLine(stitch3,slope4*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+(slope3-slope4)*stitch3+offset,1500,slope4*1500+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+(slope3-slope4)*stitch3+offset);
249  fitLine4.SetLineColor(2);
250 
251  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
252  stitchLine2 = TLine(stitch2, 0, stitch2,5.0);
253  stitchLine3 = TLine(stitch3, 0, stitch3,5.0);
254  stitchLine1.SetLineColor(2);
255  stitchLine1.SetLineStyle(7);
256  stitchLine2.SetLineColor(2);
257  stitchLine2.SetLineStyle(7);
258  stitchLine3.SetLineColor(2);
259  stitchLine3.SetLineStyle(7);
260 
261  //Graph with fitting line
262  TCanvas* c4 = new TCanvas("c4","c4");
263  susieProfile->Draw("AP");
264  fitLine1.Draw();
265  fitLine2.Draw();
266  fitLine3.Draw();
267  fitLine4.Draw();
268  stitchLine1.Draw();
269  stitchLine2.Draw();
270  stitchLine3.Draw();
271  gPad->Update();
272  gPad->SetRightMargin(0.13);
273  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
274  //palette1->SetBBoxCenterX(625);
275  gPad->Modified();
276  gPad->Update();
277  c4->Print("fittingLineOnGraph.pdf");
278 
279  //2d plot with fitting line
280  TCanvas* c5 = new TCanvas("c5","c5");
281  hist->Draw("colz");
282  fitLine1.Draw();
283  fitLine2.Draw();
284  fitLine3.Draw();
285  fitLine4.Draw();
286  stitchLine1.Draw();
287  stitchLine2.Draw();
288  stitchLine3.Draw();
289  gPad->Update();
290  gPad->SetRightMargin(0.13);
291  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
292  //palette1->SetBBoxCenterX(625);
293  gPad->Modified();
294  gPad->Update();
295  c5->Print("fittingLineOnMother.pdf");
296 
297  //2d logz plot with fitting line
298  TCanvas* c8 = new TCanvas("c8","c8");
299  c8->SetLogz();
300  hist->Draw("colz");
301  fitLine1.Draw();
302  fitLine2.Draw();
303  fitLine3.Draw();
304  fitLine4.Draw();
305  stitchLine1.Draw();
306  stitchLine2.Draw();
307  stitchLine3.Draw();
308  gPad->Update();
309  gPad->SetRightMargin(0.13);
310  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
311  //palette1->SetBBoxCenterX(625);
312  gPad->Modified();
313  gPad->Update();
314  c8->Print("fittingLineOnMotherLogz.pdf");
315 
316  //Plot of chi squared per ndf
317  TCanvas* c6 = new TCanvas("c6","c6");
318  chiValues->Draw();
319  c6->Print("graphChiSquaredperNDF.pdf");
320 
321  TFile f("muonProfile.root","RECREATE");
322  susieProfile->Write("1");
323  chiValues->Write();
324 
325 
326 
327 }
328 
329 
330 
c2
Definition: demo5.py:33
Float_t f1
void MakingProfileND()
OStream cout
Definition: OStream.cxx:6
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24