MakeProfileMuEFD.C
Go to the documentation of this file.
1 //Run using: root -b -q -l MakeProfileMuEFD.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 = 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("./2DPlotsForFittingFD.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("MuonE_hist");
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  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
52 
53  // make projection in Y view
54  TH1D* py = hist->ProjectionY("py",iw,iw,"");
55  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
56  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
57 
58  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
59 
60  double totalEntries = py->GetEntries();
61  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
62 
63  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
64 
65  //Find the location of width at half max
66  int yyhighbin = maxbin+1;
67  int yylowbin = maxbin-1;
68  // find upper bin of half max
69  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
70  yyhighbin = ihigh;
71  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
72  }
73  // find lower bin of half max
74  for(int ilow = maxbin-1; ilow > 0; --ilow){
75  yylowbin = ilow;
76  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
77  }
78 
79  int multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
80  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
81  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
82 
83  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
84  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
85 
86  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
87 
88  TF1 *f1 = new TF1("customGaus","gaus",0,5);
89  py->Fit("customGaus","OQ","",lowSide,highSide); //O and Q mean do not plot result, quiet mode
90  double mean = f1->GetParameter(1);
91  double error = f1->GetParError(1);
92  double chi = f1->GetChisquare();
93  double ndf = f1->GetNDF();
94 
95  if (ndf > 0) {chiValues->Fill(chi/ndf);}
96  else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
97 
98  if ((chi/ndf) > 10) {std::cout<<" We have chi / ndf of : "<<chi/ndf<<" in the bin: "<<iw<<std::endl;}
99 
100  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
101 
102  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
103  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
104  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
105  //Setting the error on the mean reported by the fit as error in y
106  eyxl[iw-1] = error;
107  eyxh[iw-1] = error;
108 
109  }//End of loop over bins in 2D histogram
110 
111  //Making graph from points
112  TGraphAsymmErrors* customProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
113 
114  customProfile->GetXaxis()->SetNoExponent(kTRUE);
115  customProfile->GetXaxis()->SetTitle("Reco muon track length (m)");
116  customProfile->GetYaxis()->SetTitle("True muon energy (GeV)");
117  customProfile->GetXaxis()->CenterTitle();
118  customProfile->GetYaxis()->CenterTitle();
119  customProfile->SetMarkerStyle(6);
120  customProfile->GetXaxis()->SetLimits(0,15);
121  customProfile->GetHistogram()->SetMaximum(5);
122  customProfile->GetHistogram()->SetMinimum(0);
123  customProfile->SetTitle("");
124 
125  //Overlays 2d plot and graph
126  TCanvas* c1 = new TCanvas("c1","c1");
127  customProfile->SetMarkerColor(18);
128  customProfile->SetLineColor(18);
129  hist->Draw("colz");
130  customProfile->Draw("P");
131  gPad->Update();
132  gPad->SetRightMargin(0.13);
133  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
134  //palette1->SetBBoxCenterX(625);
135  gPad->Modified();
136  gPad->Update();
137  c1->Print("muonMotherAndGraph.pdf");
138 
139  //Overlays 2d loqz plot and graph
140  TCanvas* c2 = new TCanvas("c2","c2");
141  c2->SetLogz();
142  customProfile->SetMarkerColor(13);
143  customProfile->SetLineColor(13);
144  hist->Draw("colz");
145  customProfile->Draw("P");
146  gPad->Update();
147  gPad->SetRightMargin(0.13);
148  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
149  //palette1->SetBBoxCenterX(625);
150  gPad->Modified();
151  gPad->Update();
152  c2->Print("muonMotherLogzAndGraph.pdf");
153 
154  //Just the graph
155  TCanvas* c3 = new TCanvas("c3","c3");
156  customProfile->SetMarkerColor(1);
157  customProfile->SetLineColor(1);
158  customProfile->Draw("AP");
159  gPad->Update();
160  gPad->SetRightMargin(0.13);
161  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
162  //palette1->SetBBoxCenterX(625);
163  gPad->Modified();
164  gPad->Update();
165  c3->Print("muonGraph.pdf");
166 
167  //Drawing fit over the 2d plot and graph
168 
169 
170  //New tuning values rounded for errors
171  double stitch1 = 3.34; // m
172  double stitch2 = 5.39; // m
173  double stitch3 = 10.64; // m
174  double offset = 0.1503; // GeV
175  double slope1 = 0.1910; // GeV/m
176  double slope2 = 0.1987; // GeV/m
177  double slope3 = 0.2039; // GeV/m
178  double slope4 = 0.2159; // GeV/m
179 
180  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
181  fitLine1.SetLineColor(2);
182 
183  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,stitch2,slope2*stitch2+(slope1-slope2)*stitch1+offset);
184  fitLine2.SetLineColor(2);
185 
186  fitLine3 = TLine(stitch2,slope3*stitch2+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset,stitch3,slope3*stitch3+(slope1-slope2)*stitch1+(slope2-slope3)*stitch2+offset);
187  fitLine3.SetLineColor(2);
188 
189  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);
190  fitLine4.SetLineColor(2);
191 
192  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
193  stitchLine2 = TLine(stitch2, 0, stitch2,5.0);
194  stitchLine3 = TLine(stitch3, 0, stitch3,5.0);
195  stitchLine1.SetLineColor(2);
196  stitchLine1.SetLineStyle(7);
197  stitchLine2.SetLineColor(2);
198  stitchLine2.SetLineStyle(7);
199  stitchLine3.SetLineColor(2);
200  stitchLine3.SetLineStyle(7);
201 
202  //Graph with fitting line
203  TCanvas* c4 = new TCanvas("c4","c4");
204  customProfile->Draw("AP");
205  fitLine1.Draw();
206  fitLine2.Draw();
207  fitLine3.Draw();
208  fitLine4.Draw();
209  stitchLine1.Draw();
210  stitchLine2.Draw();
211  stitchLine3.Draw();
212  gPad->Update();
213  gPad->SetRightMargin(0.13);
214  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
215  //palette1->SetBBoxCenterX(625);
216  gPad->Modified();
217  gPad->Update();
218  c4->Print("fittingLineOnGraph.pdf");
219 
220  //2d plot with fitting line
221  TCanvas* c5 = new TCanvas("c5","c5");
222  hist->Draw("colz");
223  fitLine1.Draw();
224  fitLine2.Draw();
225  fitLine3.Draw();
226  fitLine4.Draw();
227  stitchLine1.Draw();
228  stitchLine2.Draw();
229  stitchLine3.Draw();
230  gPad->Update();
231  gPad->SetRightMargin(0.13);
232  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
233  //palette1->SetBBoxCenterX(625);
234  gPad->Modified();
235  gPad->Update();
236  c5->Print("fittingLineOnMother.pdf");
237 
238  //2d logz plot with fitting line
239  TCanvas* c8 = new TCanvas("c8","c8");
240  c8->SetLogz();
241  hist->Draw("colz");
242  fitLine1.Draw();
243  fitLine2.Draw();
244  fitLine3.Draw();
245  fitLine4.Draw();
246  stitchLine1.Draw();
247  stitchLine2.Draw();
248  stitchLine3.Draw();
249  gPad->Update();
250  gPad->SetRightMargin(0.13);
251  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
252  //palette1->SetBBoxCenterX(625);
253  gPad->Modified();
254  gPad->Update();
255  c8->Print("fittingLineOnMotherLogz.pdf");
256 
257  //Plot of chi squared per ndf
258  TCanvas* c6 = new TCanvas("c6","c6");
259  chiValues->Draw();
260  c6->Print("graphChiSquaredperNDF.pdf");
261 
262  TFile f("muonProfile.root","RECREATE");
263  customProfile->Write("1");
264  chiValues->Write();
265 
266 }
c2
Definition: demo5.py:33
Float_t f1
OStream cout
Definition: OStream.cxx:6
void MakeProfileMuEFD()
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24