DedxCompBirk.C
Go to the documentation of this file.
1 #include <iostream>
2 #include <fstream>
3 #include <iomanip>
4 #include <sstream>
5 #include <string>
6 #include <cstring>
7 
8 #include "TFile.h"
9 #include "TChain.h"
10 #include "TH2.h"
11 #include "TTree.h"
12 #include "TCanvas.h"
13 #include "TLegend.h"
14 #include "TProfile.h"
15 #include "TF1.h"
16 #include "TROOT.h"
17 #include "TStyle.h"
18 #include "TGaxis.h"
19 
20 
21 TH1D* MakeRatio(TH1D* num, TH1D* denom,const char* name)
22 {
23 
24  TH1D* ratio = new TH1D(name,";;Data/MC",num->GetNbinsX(),num->GetXaxis()->GetXbins()->GetArray());
25 
26  ratio->GetXaxis()->SetTitle(num->GetXaxis()->GetTitle());
27 
28  ratio->Divide(num,denom);
29 
30  ratio->SetLineColor(denom->GetLineColor());
31  ratio->SetMarkerColor(denom->GetMarkerColor());
32  ratio->SetMarkerSize(denom->GetMarkerSize());
33 
34  return ratio;
35 
36 }
37 
38 
39 TCanvas* MakeSplitRatio(TCanvas* top, TCanvas* bot,const char* name)
40 {
41 
42  TCanvas* can = new TCanvas(name,name,600,550);
43 
44  double pixels = 18.0;
45  double labPix = 15.0;
46  double labOff = 5.0;
47 
48  can->SetBottomMargin(0);
49  can->Divide(1,2,0,0);
50  TVirtualPad* can_1 = can->cd(1);
51  can_1->SetPad(0,0.4,1,1);
52  can_1->SetFillStyle(0);
53  top->DrawClonePad();
54  can_1->SetRightMargin(0.05);
55  can_1->SetTopMargin(0.1);
56  can_1->SetBottomMargin(0.005);
57  can_1->SetLeftMargin(0.13);
58  can_1->Update();
59  double textSizeFac = 1.0/(fabs(can_1->YtoAbsPixel(can_1->GetY2()) - can_1->YtoAbsPixel(can_1->GetY1())));
60  double marginSize = 80;
61  double textsizeA = pixels/(can_1->GetAbsHNDC()*can_1->GetWh());
62  double textsizeAl = labPix/(can_1->GetAbsHNDC()*can_1->GetWh());
63  double textsizeAo = labOff/(can_1->GetWw() * can_1->GetAbsWNDC());
64 
65  TList* primlist = can_1->GetListOfPrimitives();
66  for(int ilist = 0; ilist < primlist->GetSize(); ++ilist){
67  TObject* prim = primlist->At(ilist);
68  if(prim->InheritsFrom("TH1")){
69  TH1* primhist = (TH1*)prim;
70  primhist->GetYaxis()->SetTitleSize(textsizeA);
71  primhist->GetYaxis()->SetLabelSize(textsizeAl);
72  primhist->GetYaxis()->SetLabelOffset(textsizeAo);
73  primhist->GetYaxis()->SetTitleOffset(0.035 / textsizeA);
74  primhist->GetXaxis()->SetNoExponent(kTRUE);
75  }
76  }
77 
78  can_1->Update();
79 
80  TVirtualPad* can_2 = can->cd(2);
81  can_2->SetPad(0,0,1,0.4);
82  bot->DrawClonePad();
83  can_2->SetTopMargin(0.025);
84  can_2->SetBottomMargin(0.22);
85  can_2->SetRightMargin(0.05);
86  can_2->SetLeftMargin(0.13);
87  can_2->Update();
88 
89  double textsize = pixels/(can_2->GetAbsHNDC()*can_2->GetWh());
90  double textsizel = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
91  double textsizeo = labOff/(can_2->GetWw() * can_2->GetAbsWNDC());
92  double textsizex = pixels/(can_2->GetAbsHNDC()*can_2->GetWh());
93  double textsizelx = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
94  textsizeAl = labPix/(can_2->GetAbsHNDC()*can_2->GetWh());
95  double textsizeox = labOff/(can_2->GetWw() * can_2->GetAbsWNDC());
96 
97  TList* primlist_2 = can_2->GetListOfPrimitives();
98  for(int ilist = 0; ilist < primlist_2->GetSize(); ++ilist){
99  TObject* prim = primlist_2->At(ilist);
100  if(prim->InheritsFrom("TH1")){
101  TH1* primhist = (TH1*)prim;
102  primhist->GetYaxis()->SetTitleSize(textsize);
103  primhist->GetYaxis()->SetLabelSize(textsizel);
104  primhist->GetYaxis()->SetLabelOffset(textsizeAo);
105  primhist->GetXaxis()->SetTitleSize(textsizex);
106  primhist->GetXaxis()->SetLabelSize(textsizelx);
107  primhist->GetXaxis()->SetLabelOffset(textsizeox);
108  primhist->GetYaxis()->SetTitleOffset(0.035 / textsize);
109  primhist->GetXaxis()->SetTitleOffset(0.9);
110  primhist->GetXaxis()->SetNoExponent(kTRUE);
111  }
112  }
113 
114  can_2->Update();
115 
116  return can;
117 
118 }
119 
120 void FindTruncate(TH1D* hist, double& low, double& high)
121 {
122  //std::cout<<"BEGINING OF TRUNCATE FUNCATION!"<<std::endl;
123  double percentCut = 0.05;
124  double histMin = hist->GetXaxis()->GetXmin();
125  double histMax = hist->GetXaxis()->GetXmax();
126  int minbin = hist->GetXaxis()->FindFixBin(histMin);
127  int maxbin = hist->GetXaxis()->FindFixBin(histMax)-1;
128  double integral = hist->Integral(minbin,maxbin);
129  //std::cout<<"Total integral is: "<<integral<<std::endl;
130  // if (integral == 63){
131  //TCanvas* cs = new TCanvas("cs","cs");
132  //hist->Draw();
133  //cs->Print("test.pdf");
134  //}
135  double lowCut = integral*percentCut;
136  double highCut = integral*(1-percentCut);
137 
138  int lowbin = 0;
139  int highbin = 0;
140  bool foundLow = false;
141  double ongoingIntegral = 0.0;
142 
143  for(int ibin = minbin; ibin <= maxbin; ++ibin){
144  double bincont = hist->GetBinContent(ibin);
145  //std::cout<<"bin: "<<ibin<<" with content: "<<bincont<<std::endl;
146  if (!foundLow){
147  ongoingIntegral += bincont;
148  //std::cout<<" ongoingIntegral: "<<ongoingIntegral<<std::endl;
149  if (ongoingIntegral > lowCut){
150  lowbin = ibin;
151  foundLow = true;
152  //std::cout<<"Found Low! Set lowbin : "<<lowbin<<std::endl;
153  }
154  }
155  else{
156  if ((ongoingIntegral<highCut)&&((ongoingIntegral+bincont)>=highCut)){
157  ongoingIntegral += bincont;
158  highbin = ibin;
159  //std::cout<<"Found High! Set highbin : "<<highbin<<std::endl;
160  }
161  else if (ongoingIntegral<highCut){
162  ongoingIntegral += bincont;
163  //std::cout<<" ongoingIntegral: "<<ongoingIntegral<<std::endl;
164  }
165  }//Already found low, now looking for high bin
166  //Now need to set low, high to exclude the lowbin and highbin
167  low = hist->GetXaxis()->GetBinUpEdge(lowbin);
168  high = hist->GetXaxis()->GetBinLowEdge(highbin);
169  }
170  return;
171 }
172 
173 double Mode(TH1D* hist)
174 {
175 
176  double modex = hist->GetXaxis()->GetBinCenter(hist->GetMaximumBin());
177  //std::cout<<"In MODE with Modex: "<<modex<<std::endl;
178  //SL -- adding this line, because I think it needs something ....
179  return modex;
180 }
181 
182 double ModeRestrict(TH1D* hist, double minx, double maxx)
183 {
184 
185  double mode = Mode(hist);
186  //std::cout<<"In mode restrict with Mode: "<<mode<<" maxx: "<<maxx<<" minx: "<<minx<<std::endl;
187  if( mode > maxx || mode < minx){
188  // loop over bins between maxx and minx to find the mode in the restricted range;
189  //std::cout<<"HERE HERE HERE!!!"<<std::endl;
190  int minxbin = hist->GetXaxis()->FindFixBin(minx);
191  int maxxbin = hist->GetXaxis()->FindFixBin(maxx);
192  double modeval = hist->GetBinContent(minxbin);
193  int modebin = minxbin;
194 // std::cout<<"minxbin: "<<minxbin<<" maxxbin: "<<maxxbin<<" modeval: "<<modeval<<std::endl;
195  for(int ibin = minxbin+1; ibin <= maxxbin; ++ibin){
196  double bincont = hist->GetBinContent(ibin);
197 // std::cout<<"ibin: "<<ibin<<" bincont: "<<bincont<<std::endl;
198  if(bincont > modeval){
199  modebin = ibin;
200  modeval = bincont;
201  //std::cout<<"modeval: "<<modeval<<" modebin: "<<modebin<<std::endl;
202  }
203  }
204  mode = hist->GetXaxis()->GetBinCenter(modebin);
205  }
206 
207  return mode;
208 
209 }
210 
211 
212 void GausFit(TH1D* hist, double& mean, double& err, double minrange, double maxrange, double lowsig, double highsig, bool useTruncate)
213 {
214  // std::cout<<" We are inside GausFit! "<<std::endl;
215  //std::cout<<" we have minrange = "<<minrange<<" and maxrange = "<<maxrange<<std::endl;
216  //std::cout<<" we have lowsig = "<<lowsig<<" and highsig = "<<highsig<<std::endl;
217  double mode = ModeRestrict(hist,minrange,maxrange);
218  //std::cout<<" we have mode = "<<mode<<std::endl;
219  double fitminx = -1.0;
220  double fitmaxx = -1.0;
221  if (useTruncate){
222  fitminx = lowsig;
223  fitmaxx = highsig;
224  }
225  else{
226  fitminx = mode-lowsig;
227  fitmaxx = mode+highsig;
228  }
229  //std::cout<<" we have fitminx = "<<fitminx<<" and fitmaxx = "<<fitmaxx<<std::endl;
230  if(fitminx < minrange){ fitminx = minrange; }
231  if(fitmaxx > maxrange){ fitmaxx = maxrange; }
232  //std::cout<<" we have fitminx = "<<fitminx<<" and fitmaxx = "<<fitmaxx<<std::endl;
233  TF1* fgaus = new TF1("gaus","gaus(0)",fitminx,fitmaxx);
234  //std::cout<<"fitminx: "<<fitminx<<" fitmaxx: "<<fitmaxx<<" mode: "<<mode<<std::endl;
235  fgaus->SetParameters(hist->GetMaximum(),mode,(fitmaxx-fitminx)/2.0);
236  int minbin = hist->GetXaxis()->FindFixBin(fitminx);
237  int maxbin = hist->GetXaxis()->FindFixBin(fitmaxx);
238  double integral = hist->Integral(minbin,maxbin);
239  mean = -1;
240  err = 0.000001;
241  if(integral > 100){
242  //std::cout<<" We are inside the fit loop! "<<std::endl;
243  int fitstatus = hist->Fit(fgaus,"RIQ");
244  if(fitstatus == 0){
245  mean = fgaus->GetParameter(1);
246  err = fgaus->GetParError(1);
247  if(mean > maxrange || mean < minrange){
248  mean = -1;
249  err = 0.000001;
250  }
251  //err = fgaus->GetParameter(2);
252  }
253  }
254 
255  // std::cout<<"mean: "<<mean<<" err: "<<err<<std::endl;
256 
257 }
258 
259 TH1D* PeakProjectionY(TH2F* hist,const char* name, double minrange = -1.0, double maxrange = -1.0, double lowsig = 0.001, double highsig = 0.001, bool useTruncate = true)
260 {
261  //std::cout<<"In peak projection!"<<std::endl;
262  TH1D* peakhist = new TH1D(name,"",hist->GetNbinsX(),hist->GetXaxis()->GetXbins()->GetArray());
263  peakhist->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
264  peakhist->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
265  //std::cout<<"And the axis title is: "<<hist->GetXaxis()->GetTitle()<<std::endl;
266  double fitmin = minrange;
267  double fitmax = maxrange;
268  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
269  if(minrange == -1.0){ fitmin = hist->GetXaxis()->GetXmin(); }
270  if(maxrange == -1.0){ fitmax = hist->GetXaxis()->GetXmax(); }
271  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
272  //std::cout<<"Number of x bins: "<<hist->GetNbinsX()<<std::endl;
273  for(int ixbin = 1; ixbin <= hist->GetNbinsX(); ++ixbin){
274  // get the projection of the 2d hist
275  TH1D* projection = hist->ProjectionY("",ixbin,ixbin,"o");
276  //std::cout<<" X bin: "<<ixbin<<" Entries: "<<projection->GetEntries()<<std::endl;
277  if(projection->GetEntries() == 0){ continue; }
278  // fit the projection to a Gaussian around the peak
279  double mean = 0.00175 ,err = 0.001;
280  double lowTruncate = -1.0;
281  double highTruncate = -1.0;
282  if (useTruncate){
283  FindTruncate(projection, lowTruncate, highTruncate);
284  if ((lowTruncate != -1.0) && (highTruncate != -1.0)){
285  lowsig = lowTruncate;
286  highsig = highTruncate;
287  }
288  }
289  //std::cout<<"Low truncate: "<<lowTruncate<<" High truncate: "<<highTruncate<<std::endl;
290  GausFit(projection,mean,err,fitmin,fitmax,lowsig,highsig,useTruncate);
291  //std::cout<<" Mean: "<<mean<<" Error: "<<err<<std::endl;
292  peakhist->SetBinContent(ixbin,mean);
293  peakhist->SetBinError(ixbin,err);
294 
295  }
296 
297  return peakhist;
298 
299 }
300 
301 
303 {
304 
305  //const char* filenamedata = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/datadedxVxFineS15-05-22.root";
306  //const char* filenamemc = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/mcdedxVxFineS15-05-22.root";
307 
308  const char* filenamedata = "/local/nova/users/lein/datadedxVxFineFA.root";
309  const char* filenamemc = "/local/nova/users/lein/mcdedxVxFineFA.root";
310  const char* filenamemcbirkB = "/local/nova/users/lein/mcdedxVxFineFABirksB.root";
311  const char* filenamemcbirkC = "/local/nova/users/lein/mcdedxVxFineFABirksC.root";
312 
313  TFile* filedata = new TFile(filenamedata,"READ");
314  TFile* filemc = new TFile(filenamemc,"READ");
315  TFile* filemcbirkB = new TFile(filenamemcbirkB,"READ");
316  TFile* filemcbirkC = new TFile(filenamemcbirkC,"READ");
317 
318  TH2F* dedxdata = (TH2F*)filedata->Get("DedxVx");
319  TH2F* dedxmc = (TH2F*)filemc->Get("DedxVx");
320  TH2F* dedxmcbirkB = (TH2F*)filemcbirkB->Get("DedxVx");
321  TH2F* dedxmcbirkC = (TH2F*)filemcbirkC->Get("DedxVx");
322 
323  //std::cout<<"This is the pointer: "<<dedxdata<<std::endl;
324 
325  gStyle->SetOptStat(00000);
326 
327  TH1D* dedxpeakdata = PeakProjectionY(dedxdata,"dedxpeakdata");
328  TH1D* dedxpeakmc = PeakProjectionY(dedxmc,"dedxpeakmc");
329  TH1D* dedxpeakmcbirkB = PeakProjectionY(dedxmcbirkB,"dedxpeakmcbirkB");
330  TH1D* dedxpeakmcbirkC = PeakProjectionY(dedxmcbirkC,"dedxpeakmcbirkC");
331 
332  //If you want to run WITHOUT the truncation and just limit the fit range to +- 0.001 around the mode, do these two lines
333  //instead of the two lines above
334  //TH1D* dedxpeakdata = PeakProjectionY(dedxdata,"dedxpeakdata",-1.0,-1.0,0.001,0.001,false);
335  //TH1D* dedxpeakmc = PeakProjectionY(dedxmc,"dedxpeakmc",-1.0,-1.0,0.001,0.001,false);
336 
337 
338  TCanvas* c0 = new TCanvas("c0","c0");
339  dedxpeakdata->Draw();
340  dedxpeakmc->Draw("same");
341  dedxpeakmcbirkB->Draw("same");
342  dedxpeakmcbirkC->Draw("same");
343  dedxpeakmc->SetLineColor(kRed);
344  dedxpeakmc->SetMarkerColor(kRed);
345  dedxpeakmcbirkB->SetLineColor(kGreen+2);
346  dedxpeakmcbirkB->SetMarkerColor(kGreen+2);
347  dedxpeakmcbirkC->SetLineColor(kBlue-4);
348  dedxpeakmcbirkC->SetMarkerColor(kBlue-4);
349  dedxpeakdata->GetYaxis()->SetRangeUser(0.0,0.003);
350  dedxpeakdata->GetXaxis()->SetNoExponent();
351  //dedxpeakdata->GetYaxis()->SetNoExponent();
352  TLegend* leg0 = new TLegend(0.65,0.60,0.85,0.85);
353  leg0->AddEntry(dedxpeakdata,"Data","LE");
354  leg0->AddEntry(dedxpeakmc,"MC","LE");
355  leg0->AddEntry(dedxpeakmcbirkB,"MC BirksB","LE");
356  leg0->AddEntry(dedxpeakmcbirkC,"MC BirksC","LE");
357  leg0->Draw("same");
358  c0->Print("dedxfine.pdf");
359 
360  TCanvas* c00 = new TCanvas("c00","c00");
361  dedxdata->Draw("colz");
362  dedxpeakdata->Draw("same");
363  dedxdata->GetYaxis()->SetRangeUser(0.0,0.003);
364  dedxdata->GetXaxis()->SetNoExponent();
365  //dedxdata->GetYaxis()->SetNoExponent();
366  c00->Print("dedxDataMother.pdf");
367 
368  TCanvas* c000 = new TCanvas("c000","c000");
369  //dedxpeakmc->SetLineColor(kBlack);
370  dedxmc->Draw("colz");
371  dedxpeakmc->Draw("same");
372  dedxmc->GetYaxis()->SetRangeUser(0.0,0.003);
373  dedxmc->GetXaxis()->SetNoExponent();
374  //dedxmc->GetYaxis()->SetNoExponent();
375  c000->Print("dedxMCMother.pdf");
376 
377  TCanvas* c000b = new TCanvas("c000b","c000b");
378  //dedxpeakmc->SetLineColor(kBlack);
379  dedxmcbirkB->Draw("colz");
380  dedxpeakmcbirkB->Draw("same");
381  dedxmcbirkB->GetYaxis()->SetRangeUser(0.0,0.003);
382  dedxmcbirkB->GetXaxis()->SetNoExponent();
383  //dedxmc->GetYaxis()->SetNoExponent();
384  c000b->Print("dedxMCBirkBMother.pdf");
385  TCanvas* c000c = new TCanvas("c000c","c000c");
386  //dedxpeakmc->SetLineColor(kBlack);
387  dedxmcbirkC->Draw("colz");
388  dedxpeakmcbirkC->Draw("same");
389  dedxmcbirkC->GetYaxis()->SetRangeUser(0.0,0.003);
390  dedxmcbirkC->GetXaxis()->SetNoExponent();
391  //dedxmc->GetYaxis()->SetNoExponent();
392  c000c->Print("dedxMCBirkCMother.pdf");
393 
394 
395  int nsteps = 1000;
396  double minscale = 0.95;
397  double maxscale = 1.05;
398 
399  TH1F* chi2match = new TH1F("chi2match",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
400  TH1F* chi2matchrestrict = new TH1F("chi2matchrestrict",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
401 
402 
403  for(int ibin = 1; ibin <= chi2match->GetNbinsX(); ++ibin){
404  // Get the scale factor
405  double scalefactor = chi2match->GetXaxis()->GetBinCenter(ibin);
406  //calculate a chi2 between data and scaled mc
407  double chi2 = 0;
408  double chi2rest = 0;
409  for(int jbin = 1; jbin <= dedxpeakdata->GetNbinsX(); ++jbin){
410  // get the dist from the end of the track
411  double dist = dedxpeakdata->GetXaxis()->GetBinCenter(jbin);
412  double expt = dedxpeakdata->GetBinContent(jbin);
413  double errdata = dedxpeakdata->GetBinError(jbin);
414  double errmc = dedxpeakmc->GetBinError(jbin);
415  double err2 = errdata*errdata+errmc*errmc;
416  double obsunscale = dedxpeakmc->GetBinContent(jbin);
417  double obs = scalefactor*dedxpeakmc->GetBinContent(jbin);
418  double diff = expt-obs;
419  //chi2+=(diff*diff/expt);
420  // remove error cases
421  if(expt < 0 || obsunscale < 0){ continue; }
422  chi2+=(diff*diff/err2);
423  // std::cout<<"dist: "<<dist<<std::endl;
424  //if(dist>100 && dist < 200){
425  if(dist>200){
426  chi2rest+=(diff*diff/err2);
427  }
428 
429  }
430 // // make a histogram to store the scaled version of mc dedxpeak histogram
431 // TString name;
432 // name+=ibin;
433 // TH1D* dedxpeakmcscale = new TH1D(name.Data(),"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
434 // dedxpeakmcscale->Add(dedxpeakmc);
435 // dedxpeakmcscale->Scale(scalefactor);
436 // dedxpeakmcscale->
437 // double chi2 = dedxpeakdata->Chi2Test(dedxpeakmcscale,"CHI2");
438 // std::cout<<"scalefactor: "<<scalefactor<<" chi2: "<<chi2<<std::endl;
439  chi2match->SetBinContent(ibin,chi2);
440  chi2matchrestrict->SetBinContent(ibin,chi2rest);
441 
442  }
443 
444  TCanvas* c00a = new TCanvas("c00a","c00a");
445  chi2match->Draw();
446  c00a->Print("chi2fine.pdf");
447 
448  double optscale = chi2match->GetXaxis()->GetBinCenter(chi2match->GetMinimumBin());
449  std::cout<<"optscale: "<<optscale<<std::endl;
450 
451  TCanvas* c000a = new TCanvas("c000a","c000a");
452  chi2matchrestrict->Draw();
453  c000a->Print("chi2restrictfine.pdf");
454 
455  double optscalerestrict = chi2matchrestrict->GetXaxis()->GetBinCenter(chi2matchrestrict->GetMinimumBin());
456  std::cout<<"Restricted optscale: "<<optscalerestrict<<std::endl;
457 
458  TH1D* dedxpeakmcscale = new TH1D("OptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
459  dedxpeakmcscale->Add(dedxpeakmc);
460  dedxpeakmcscale->Scale(optscale);
461 
462  TCanvas* c2 = new TCanvas("c2","c2");
463  dedxpeakdata->Draw();
464  dedxpeakmcscale->Draw("same");
465  dedxpeakmcscale->SetLineColor(kRed);
466  dedxpeakmcscale->SetMarkerColor(kRed);
467  dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
468  dedxpeakdata->GetXaxis()->SetNoExponent();
469  //dedxpeakdata->GetYaxis()->SetNoExponent();
470  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
471  leg0->AddEntry(dedxpeakdata,"Data","LE");
472  leg0->AddEntry(dedxpeakmcscale,"MC","LE");
473  leg0->Draw("same");
474  c2->Print("dedxoptfine.pdf");
475 
476  double nonoptscale = 0.97375; //1.483/1.55;
477 
478  TH1D* dedxpeakmcnonoptscale = new TH1D("NonOptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
479  dedxpeakmcnonoptscale->Add(dedxpeakmc);
480  dedxpeakmcnonoptscale->Scale(nonoptscale);
481 
482  TCanvas* c3 = new TCanvas("c3","c3");
483  dedxpeakdata->Draw();
484  dedxpeakmcnonoptscale->Draw("same");
485  dedxpeakmcnonoptscale->SetLineColor(kRed);
486  dedxpeakmcnonoptscale->SetMarkerColor(kRed);
487  dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
488  dedxpeakdata->GetXaxis()->SetNoExponent();
489  //dedxpeakdata->GetYaxis()->SetNoExponent();
490  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
491  leg0->AddEntry(dedxpeakdata,"Data","LE");
492  leg0->AddEntry(dedxpeakmcnonoptscale,"MC","LE");
493  leg0->Draw("same");
494  c3->Print("dedxnonoptfine.pdf");
495 
496 
497  TCanvas* c4 = new TCanvas("c4","c4");
498  TH1D* dedxratio = MakeRatio(dedxpeakdata,dedxpeakmc,"dedxratio");
499  TH1D* dedxratiobirkB = MakeRatio(dedxpeakdata,dedxpeakmcbirkB,"dedxratiobirkB");
500  TH1D* dedxratiobirkC = MakeRatio(dedxpeakdata,dedxpeakmcbirkC,"dedxratiobirkC");
501  TLine* unitline = new TLine(-1,1,1001,1);
502  unitline->SetLineWidth(2);
503  unitline->SetLineStyle(7);
504  dedxratio->Draw();
505  unitline->Draw("same");
506  dedxratiobirkB->Draw("same");
507  dedxratiobirkC->Draw("same");
508  dedxratio->GetYaxis()->SetRangeUser(0.9,1.1);
509 
510  TCanvas* c5 = MakeSplitRatio(c0,c4,"c5");
511  c5->Print("dedxratio.pdf");
512 
513 
514 
515 
516 }
const XML_Char * name
Definition: expat.h:151
void DedxCompBirk()
Definition: DedxCompBirk.C:302
enum BeamMode kRed
double Mode(TH1D *hist)
Definition: DedxCompBirk.C:173
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
void GausFit(TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig, bool useTruncate)
Definition: DedxCompBirk.C:212
double ModeRestrict(TH1D *hist, double minx, double maxx)
Definition: DedxCompBirk.C:182
const int nsteps
TH1 * ratio(TH1 *h1, TH1 *h2)
TH1D * MakeRatio(TH1D *num, TH1D *denom, const char *name)
Definition: DedxCompBirk.C:21
double chi2()
double dist
Definition: runWimpSim.h:113
c2
Definition: demo5.py:33
expt
Definition: demo5.py:34
OStream cout
Definition: OStream.cxx:6
TH1D * PeakProjectionY(TH2F *hist, const char *name, double minrange=-1.0, double maxrange=-1.0, double lowsig=0.001, double highsig=0.001, bool useTruncate=true)
Definition: DedxCompBirk.C:259
int num
Definition: f2_nu.C:119
enum BeamMode kGreen
enum BeamMode kBlue
TGeoVolume * top
Definition: make_fe_box.C:9
void FindTruncate(TH1D *hist, double &low, double &high)
Definition: DedxCompBirk.C:120
TCanvas * MakeSplitRatio(TCanvas *top, TCanvas *bot, const char *name)
Definition: DedxCompBirk.C:39