MakeProfileHadEND.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 = 115; //Number of points on graph - should be equal to or smaller than number of bins in x
25  // static const int nn = 57; //Number of points on graph - should be equal to or smaller than number of bins in x
26  float xxx[nn];
27  float exxl[nn];
28  float exxh[nn];
29  float yyx[nn];
30  float eyxl[nn];
31  float eyxh[nn];
32 
33  gStyle->SetOptStat(0);
34  TFile* file = new TFile("./2DPlotsForFittingND.root","READ");
35 
36  TH1D* chiValues = new TH1D("chiValue",";Chi squared per NDF;Fits",100, 0.0, 10.0);
37 
38  TH2D* hist = file->Get("HadE_hist");
39 
40  TCanvas* c7 = new TCanvas("c7","c7");
41  hist->GetXaxis()->CenterTitle();
42  hist->GetYaxis()->CenterTitle();
43  hist->Draw("colz");
44  gPad->Update();
45  gPad->SetRightMargin(0.13);
46  TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
47  palette1->SetBBoxCenterX(625);
48  gPad->Modified();
49  gPad->Update();
50  c7->Print("hadronccMother.pdf");
51 
52  for(int iw = 1; iw <= hist->GetNbinsX(); ++iw){
53 
54  // make projection in Y view
55  TH1D* py = hist->ProjectionY("py",iw,iw,"");
56  int maxbin = py->GetMaximumBin(); //Bin in slice with most entries
57  double maxcont = py->GetBinContent(maxbin); //Height of the bin with most entries
58 
59  if (maxcont == 0){continue; } //If no entries in the slice, don't make a graph point
60 
61  double totalEntries = py->GetEntries();
62  if (totalEntries < 30){continue; } //If not many entries in the slice, don't make a graph point
63 
64  double modex = py->GetBinCenter(maxbin); //Center location of bin with most entries
65 
66  //Find the location of width at half max
67  int yyhighbin = maxbin+1;
68  int yylowbin = maxbin-1;
69  // find upper bin of half max
70  for(int ihigh = maxbin+1; ihigh <= nn; ++ihigh){
71  yyhighbin = ihigh;
72  if(py->GetBinContent(yyhighbin) <= maxcont/2.0){ break; }
73  }
74  // find lower bin of half max
75  for(int ilow = maxbin-1; ilow > 0; --ilow){
76  yylowbin = ilow;
77  if(py->GetBinContent(yylowbin) <= maxcont/2.0){ break; }
78  }
79 
80  double multiplier = 1.5; //This defines how many times wider than width at half max I want to fit over
81  if (iw>100){multiplier = 10.0;}
82  if (iw<3){multiplier = 0.5;}
83  if ((iw<20)&&(iw>2)){multiplier = 1.0;}
84  int veryHighBin = (yyhighbin - maxbin)*multiplier+maxbin;
85  int veryLowBin = maxbin - multiplier*(maxbin-yylowbin);
86 
87  double lowSide = modex-multiplier*(modex-py->GetBinCenter(yylowbin));
88  double highSide = multiplier*(py->GetBinCenter(yyhighbin)-modex)+modex;
89 
90  xxx[iw-1] = (hist->GetXaxis()->GetBinCenter(iw)); //Center of x bin position set as graph point
91 
92  TF1 *f1 = new TF1("customGaus","gaus",0,5);
93  py->Fit("customGaus","OQ","",lowSide,highSide); //O and Q mean: do not plot result, quiet mode
94  double mean = f1->GetParameter(1);
95  double error = f1->GetParError(1);
96  double chi = f1->GetChisquare();
97  double ndf = f1->GetNDF();
98 
99  if (ndf > 0) {chiValues->Fill(chi/ndf);}
100  else if (chi > 0.000001){chiValues->Fill(chi/ndf);}
101 
102  yyx[iw-1] = mean; //Set mean of gaussian fit as y location of graph point
103 
104  //For now, setting the x error bars as zero. Can change if want to denote variable bin widths.
105  exxl[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
106  exxh[iw-1] = 0; //effVwx->GetBinWidth(iw)/2.0;
107  //Setting the error on the mean reported by the fit as error in y
108  eyxl[iw-1] = error;
109  eyxh[iw-1] = error;
110 
111  }//End of loop over bins in 2D histogram
112 
113  //Making graph from points
114  TGraphAsymmErrors* customProfile = new TGraphAsymmErrors(nn,xxx,yyx,exxl,exxh,eyxl,eyxh);
115 
116  customProfile->GetXaxis()->SetNoExponent(kTRUE);
117  customProfile->GetXaxis()->SetTitle("Visible Hadronic E (GeV)");
118  customProfile->GetYaxis()->SetTitle("True Hadronic E (GeV)");
119  customProfile->GetXaxis()->CenterTitle();
120  customProfile->GetYaxis()->CenterTitle();
121  customProfile->SetMarkerStyle(6);
122  customProfile->GetXaxis()->SetLimits(0,2.0);
123  customProfile->GetHistogram()->SetMaximum(5);
124  customProfile->GetHistogram()->SetMinimum(0);
125  customProfile->SetTitle("");
126 
127  //Overlays 2d plot and graph
128  TCanvas* c1 = new TCanvas("c1","c1");
129  customProfile->SetMarkerColor(18);
130  customProfile->SetLineColor(18);
131  hist->Draw("colz");
132  customProfile->Draw("P");
133  gPad->Update();
134  gPad->SetRightMargin(0.13);
135  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
136  //palette1->SetBBoxCenterX(625);
137  gPad->Modified();
138  gPad->Update();
139  c1->Print("hadronccMotherAndGraph.pdf");
140 
141  //Overlays 2d loqz plot and graph
142  TCanvas* c2 = new TCanvas("c2","c2");
143  c2->SetLogz();
144  customProfile->SetMarkerColor(13);
145  customProfile->SetLineColor(13);
146  hist->Draw("colz");
147  customProfile->Draw("P");
148  gPad->Update();
149  gPad->SetRightMargin(0.13);
150  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
151  //palette1->SetBBoxCenterX(625);
152  gPad->Modified();
153  gPad->Update();
154  c2->Print("hadronccMotherLogzAndGraph.pdf");
155 
156  //Just the graph
157  TCanvas* c3 = new TCanvas("c3","c3");
158  customProfile->SetMarkerColor(1);
159  customProfile->SetLineColor(1);
160  customProfile->Draw("AP");
161  gPad->Update();
162  gPad->SetRightMargin(0.13);
163  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
164  //palette1->SetBBoxCenterX(625);
165  gPad->Modified();
166  gPad->Update();
167  c3->Print("hadronccGraph.pdf");
168 
169  //Drawing fit over the 2d plot and graph
170 
171  //Old Tune Values
172  //double stitch1 = 0.225007; // GeV
173  //double offset = 0.202837; // GeV
174  //double slope1 = 1.21606; // unitless
175  //double slope2 = 1.56778; // unitless
176 
177  //New Tune Values
178  //double stitch1 = 0.0145409; // GeV
179  //double offset = 0.0607282; // GeV
180  //double slope1 = 0.642951; // unitless
181  //double slope2 = 2.07193; // unitless
182 
183  //New Tune Values with Error Rounding
184  double stitch1 = 0.015; // GeV
185  double offset = 0.061; // GeV
186  double slope1 = 0.64; // unitless
187  double slope2 = 2.072; // unitless
188 
189  fitLine1 = TLine(0,offset,stitch1,slope1*stitch1+offset);
190  fitLine1.SetLineColor(2);
191 
192  fitLine2 = TLine(stitch1,slope2*stitch1+(slope1-slope2)*stitch1+offset,2.0,slope2*2.0+(slope1-slope2)*stitch1+offset);
193  fitLine2.SetLineColor(2);
194 
195 
196  stitchLine1 = TLine(stitch1, 0, stitch1,5.0);
197  stitchLine1.SetLineColor(2);
198  stitchLine1.SetLineStyle(7);
199 
200 
201  //Graph with fitting line
202  TCanvas* c4 = new TCanvas("c4","c4");
203  customProfile->Draw("AP");
204  fitLine1.Draw();
205  fitLine2.Draw();
206  stitchLine1.Draw();
207  gPad->Update();
208  gPad->SetRightMargin(0.13);
209  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
210  //palette1->SetBBoxCenterX(625);
211  gPad->Modified();
212  gPad->Update();
213  c4->Print("fittingLineOnGraphhadcc.pdf");
214 
215  //2d plot with fitting line
216  TCanvas* c5 = new TCanvas("c5","c5");
217  hist->Draw("colz");
218  fitLine1.Draw();
219  fitLine2.Draw();
220  stitchLine1.Draw();
221  gPad->Update();
222  gPad->SetRightMargin(0.13);
223  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
224  //palette1->SetBBoxCenterX(625);
225  gPad->Modified();
226  gPad->Update();
227  c5->Print("fittingLineOnMotherhadcc.pdf");
228 
229  //Plot of chi squared per ndf
230  TCanvas* c6 = new TCanvas("c6","c6");
231  chiValues->Draw();
232  c6->Print("graphChiSquaredperNDFhadcc.pdf");
233 
234  //2d plot with fitting line on logz
235  TCanvas* c8 = new TCanvas("c8","c8");
236  c8->SetLogz();
237  hist->Draw("colz");
238  fitLine1.Draw();
239  fitLine2.Draw();
240  stitchLine1.Draw();
241  gPad->Update();
242  gPad->SetRightMargin(0.13);
243  //TPaletteAxis *palette1 = (TPaletteAxis*)hist->GetListOfFunctions()->FindObject("palette");
244  //palette1->SetBBoxCenterX(625);
245  gPad->Modified();
246  gPad->Update();
247  c8->Print("fittingLineOnMotherLogzhadcc.pdf");
248 
249  TFile f("hadronccProfile.root","RECREATE");
250  customProfile->Write("1");
251  chiValues->Write();
252 
253 
254 
255 }
256 
257 
258 
c2
Definition: demo5.py:33
Float_t f1
OStream cout
Definition: OStream.cxx:6
void MakeProfileHadEND()
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24