MakingProfile.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, 10.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,2000);
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  //Old tuning values
229  //double stitch1 = 317.721; // cm
230  //double stitch2 = 539.019; // cm
231  //double stitch3 = 1012.83; // cm
232  //double offset = 0.150468; // GeV
233  //double slope1 = 0.00190560; // GeV/cm
234  //double slope2 = 0.00198726; // GeV/cm
235  //double slope3 = 0.00203830; // GeV/cm
236  //double slope4 = 0.00212812; // GeV/cm
237 
238  //New tuning values
239  //double stitch1 = 334.212; // cm
240  //double stitch2 = 539.481; // cm
241  //double stitch3 = 1064.40; // cm
242  //double offset = 0.150268; // GeV
243  //double slope1 = 0.00191025; // GeV/cm
244  //double slope2 = 0.00198716; // GeV/cm
245  //double slope3 = 0.00203910; // GeV/cm
246  //double slope4 = 0.00215887; // GeV/cm
247 
248  //New tuning values rounded for errors
249  double stitch1 = 334; // cm
250  double stitch2 = 539; // cm
251  double stitch3 = 1064; // cm
252  double offset = 0.1503; // GeV
253  double slope1 = 0.001910; // GeV/cm
254  double slope2 = 0.001987; // GeV/cm
255  double slope3 = 0.002039; // GeV/cm
256  double slope4 = 0.002159; // GeV/cm
257 
258  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
259  fitLine1.SetLineColor(2);
260 
261  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,stitch2,slope2*stitch2+(slope1-slope2)*stitch1+offset);
262  fitLine2.SetLineColor(2);
263 
264  fitLine3 = TLine(stitch2,slope3*stitch2+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset,stitch3,slope3*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset);
265  fitLine3.SetLineColor(2);
266 
267  fitLine4 = TLine(stitch3,slope4*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+(slope3-slope4)*stitch3+offset,2000,slope4*2000+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+(slope3-slope4)*stitch3+offset);
268  fitLine4.SetLineColor(2);
269 
270  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
271  stitchLine2 = TLine(stitch2, 0, stitch2,5.0);
272  stitchLine3 = TLine(stitch3, 0, stitch3,5.0);
273  stitchLine1.SetLineColor(2);
274  stitchLine1.SetLineStyle(7);
275  stitchLine2.SetLineColor(2);
276  stitchLine2.SetLineStyle(7);
277  stitchLine3.SetLineColor(2);
278  stitchLine3.SetLineStyle(7);
279 
280  //Graph with fitting line
281  TCanvas* c4 = new TCanvas("c4","c4");
282  susieProfile->Draw("AP");
283  fitLine1.Draw();
284  fitLine2.Draw();
285  fitLine3.Draw();
286  fitLine4.Draw();
287  stitchLine1.Draw();
288  stitchLine2.Draw();
289  stitchLine3.Draw();
290  gPad->Update();
291  gPad->SetRightMargin(0.13);
292  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
293  //palette1->SetBBoxCenterX(625);
294  gPad->Modified();
295  gPad->Update();
296  c4->Print("fittingLineOnGraph.pdf");
297 
298  //2d plot with fitting line
299  TCanvas* c5 = new TCanvas("c5","c5");
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  c5->Print("fittingLineOnMother.pdf");
315 
316  //2d logz plot with fitting line
317  TCanvas* c8 = new TCanvas("c8","c8");
318  c8->SetLogz();
319  hist->Draw("colz");
320  fitLine1.Draw();
321  fitLine2.Draw();
322  fitLine3.Draw();
323  fitLine4.Draw();
324  stitchLine1.Draw();
325  stitchLine2.Draw();
326  stitchLine3.Draw();
327  gPad->Update();
328  gPad->SetRightMargin(0.13);
329  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
330  //palette1->SetBBoxCenterX(625);
331  gPad->Modified();
332  gPad->Update();
333  c8->Print("fittingLineOnMotherLogz.pdf");
334 
335  //Plot of chi squared per ndf
336  TCanvas* c6 = new TCanvas("c6","c6");
337  chiValues->Draw();
338  c6->Print("graphChiSquaredperNDF.pdf");
339 
340  TFile f("muonProfile.root","RECREATE");
341  susieProfile->Write("1");
342  chiValues->Write();
343 
344 
345 
346 }
347 
348 
349 
void MakingProfile()
Definition: MakingProfile.C:20
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