MakingProfileActCatcherND.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",100, 0.0, 10.0);
36 
37  TH2D* hist = file->Get("3");
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_actCatcher.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 < 154){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,3);
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,400);
180  susieProfile->GetHistogram()->SetMaximum(3);
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_actCatcher.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_actCatcher.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_actCatcher.pdf");
225 
226  //Drawing fit over the 2d plot and graph
227 
228  //Old Tuning Values
229  double offset = 0.171226; // GeV
230  double slope1 = 0.00535381; // GeV/cm
231 
232  //New Tuning Values with Error Rounding
233  //double offset = 0.152; // GeV
234  //double slope1 = 0.00536; // GeV/cm
235 
236  fitLine1 = TLine(0,offset,400,slope1*400+offset);
237  fitLine1.SetLineColor(2);
238 
239 
240  //Graph with fitting line
241  TCanvas* c4 = new TCanvas("c4","c4");
242  susieProfile->Draw("AP");
243  fitLine1.Draw();
244  gPad->Update();
245  gPad->SetRightMargin(0.13);
246  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
247  //palette1->SetBBoxCenterX(625);
248  gPad->Modified();
249  gPad->Update();
250  c4->Print("fittingLineOnGraph_actCatcher.pdf");
251 
252  //2d plot with fitting line
253  TCanvas* c5 = new TCanvas("c5","c5");
254  hist->Draw("colz");
255  fitLine1.Draw();
256  gPad->Update();
257  gPad->SetRightMargin(0.13);
258  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
259  //palette1->SetBBoxCenterX(625);
260  gPad->Modified();
261  gPad->Update();
262  c5->Print("fittingLineOnMother_actCatcher.pdf");
263 
264  //2d logz plot with fitting line
265  TCanvas* c8 = new TCanvas("c8","c8");
266  c8->SetLogz();
267  hist->Draw("colz");
268  fitLine1.Draw();
269  gPad->Update();
270  gPad->SetRightMargin(0.13);
271  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
272  //palette1->SetBBoxCenterX(625);
273  gPad->Modified();
274  gPad->Update();
275  c8->Print("fittingLineOnMotherLogz_actCatcher.pdf");
276 
277  //Plot of chi squared per ndf
278  TCanvas* c6 = new TCanvas("c6","c6");
279  chiValues->Draw();
280  c6->Print("graphChiSquaredperNDF_actCatcher.pdf");
281 
282  TFile f("muonProfile.root","RECREATE");
283  susieProfile->Write("1");
284  chiValues->Write();
285 
286 
287 
288 }
289 
290 
291 
void MakingProfileActCatcherND()
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