Functions
ProtonDedxCompBirk.C File Reference
#include <iostream>
#include <fstream>
#include <iomanip>
#include <sstream>
#include <string>
#include <cstring>
#include "TFile.h"
#include "TChain.h"
#include "TH2.h"
#include "TTree.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "TProfile.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGaxis.h"

Go to the source code of this file.

Functions

TH1D * MakeRatio (TH1D *num, TH1D *denom, const char *name)
 
TCanvas * MakeSplitRatio (TCanvas *top, TCanvas *bot, const char *name)
 
double Mode (TH1D *hist)
 
double ModeRestrict (TH1D *hist, double minx, double maxx)
 
void GausFit (TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig)
 
TH1D * PeakProjectionY (TH2F *hist, const char *name, double minrange=-1.0, double maxrange=-1.0, double lowsig=0.001, double highsig=0.001)
 
void ProtonDedxCompBirk ()
 

Function Documentation

void GausFit ( TH1D *  hist,
double &  mean,
double &  err,
double  minrange,
double  maxrange,
double  lowsig,
double  highsig 
)

Definition at line 159 of file ProtonDedxCompBirk.C.

References febshutoff_auto::integral, submit_nova_art::mode, and ModeRestrict().

Referenced by PeakProjectionY().

160 {
161  //std::cout<<" We are inside GausFit! "<<std::endl;
162  double mode = ModeRestrict(hist,minrange,maxrange);
163  double fitminx = mode-lowsig;
164  double fitmaxx = mode+highsig;
165  if(fitminx < minrange){ fitminx = minrange; }
166  if(fitmaxx > maxrange){ fitmaxx = maxrange; }
167 
168  TF1* fgaus = new TF1("gaus","gaus(0)",fitminx,fitmaxx);
169  //std::cout<<"fitminx: "<<fitminx<<" fitmaxx: "<<fitmaxx<<" mode: "<<mode<<std::endl;
170  fgaus->SetParameters(hist->GetMaximum(),mode,(fitmaxx-fitminx)/2.0);
171  int minbin = hist->GetXaxis()->FindFixBin(fitminx);
172  int maxbin = hist->GetXaxis()->FindFixBin(fitmaxx);
173  double integral = hist->Integral(minbin,maxbin);
174  mean = -1;
175  err = 0.000001;
176  //std::cout<<"And the integral is: "<<integral<<std::endl;
177  if(integral > 100){
178  //std::cout<<" We are inside the fit loop! "<<std::endl;
179  int fitstatus = hist->Fit(fgaus,"RIQ");
180  if(fitstatus == 0){
181  mean = fgaus->GetParameter(1);
182  err = fgaus->GetParError(1);
183  if(mean > maxrange || mean < minrange){
184  mean = -1;
185  err = 0.000001;
186  }
187  //err = fgaus->GetParameter(2);
188  }
189  }
190 
191  //std::cout<<"mean: "<<mean<<" err: "<<err<<std::endl;
192 
193 }
double ModeRestrict(TH1D *hist, double minx, double maxx)
TH1D* MakeRatio ( TH1D *  num,
TH1D *  denom,
const char *  name 
)

Definition at line 21 of file ProtonDedxCompBirk.C.

References ratio().

Referenced by ProtonDedxCompBirk().

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 }
const XML_Char * name
Definition: expat.h:151
TH1 * ratio(TH1 *h1, TH1 *h2)
int num
Definition: f2_nu.C:119
TCanvas* MakeSplitRatio ( TCanvas *  top,
TCanvas *  bot,
const char *  name 
)

Definition at line 39 of file ProtonDedxCompBirk.C.

References allTimeWatchdog::can, and stan::math::fabs().

Referenced by ProtonDedxCompBirk().

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 }
const XML_Char * name
Definition: expat.h:151
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TGeoVolume * top
Definition: make_fe_box.C:9
double Mode ( TH1D *  hist)

Definition at line 120 of file ProtonDedxCompBirk.C.

Referenced by ModeRestrict().

121 {
122 
123  double modex = hist->GetXaxis()->GetBinCenter(hist->GetMaximumBin());
124  //std::cout<<"In MODE with Modex: "<<modex<<std::endl;
125  //SL -- adding this line, because I think it needs something ....
126  return modex;
127 }
double ModeRestrict ( TH1D *  hist,
double  minx,
double  maxx 
)

Definition at line 129 of file ProtonDedxCompBirk.C.

References make_syst_table_plots::ibin, Mode(), and submit_nova_art::mode.

Referenced by GausFit().

130 {
131 
132  double mode = Mode(hist);
133  //std::cout<<"In mode restrict with Mode: "<<mode<<" maxx: "<<maxx<<" minx: "<<minx<<std::endl;
134  if( mode > maxx || mode < minx){
135  // loop over bins between maxx and minx to find the mode in the restricted range;
136  //std::cout<<"HERE HERE HERE!!!"<<std::endl;
137  int minxbin = hist->GetXaxis()->FindFixBin(minx);
138  int maxxbin = hist->GetXaxis()->FindFixBin(maxx);
139  double modeval = hist->GetBinContent(minxbin);
140  int modebin = minxbin;
141 // std::cout<<"minxbin: "<<minxbin<<" maxxbin: "<<maxxbin<<" modeval: "<<modeval<<std::endl;
142  for(int ibin = minxbin+1; ibin <= maxxbin; ++ibin){
143  double bincont = hist->GetBinContent(ibin);
144 // std::cout<<"ibin: "<<ibin<<" bincont: "<<bincont<<std::endl;
145  if(bincont > modeval){
146  modebin = ibin;
147  modeval = bincont;
148  //std::cout<<"modeval: "<<modeval<<" modebin: "<<modebin<<std::endl;
149  }
150  }
151  mode = hist->GetXaxis()->GetBinCenter(modebin);
152  }
153 
154  return mode;
155 
156 }
double Mode(TH1D *hist)
TH1D* PeakProjectionY ( TH2F *  hist,
const char *  name,
double  minrange = -1.0,
double  maxrange = -1.0,
double  lowsig = 0.001,
double  highsig = 0.001 
)

Definition at line 195 of file ProtonDedxCompBirk.C.

References GausFit(), and extractScale::mean.

Referenced by ProtonDedxCompBirk().

196 {
197  //std::cout<<"In peak projection!"<<std::endl;
198  TH1D* peakhist = new TH1D(name,"",hist->GetNbinsX(),hist->GetXaxis()->GetXbins()->GetArray());
199  peakhist->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
200  peakhist->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
201  //std::cout<<"And the axis title is: "<<hist->GetXaxis()->GetTitle()<<std::endl;
202  double fitmin = minrange;
203  double fitmax = maxrange;
204  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
205  if(minrange == -1.0){ fitmin = hist->GetXaxis()->GetXmin(); }
206  if(maxrange == -1.0){ fitmax = hist->GetXaxis()->GetXmax(); }
207  //std::cout<<"Minrange: "<<fitmin<<" and maxrange: "<<fitmax<<std::endl;
208  //std::cout<<"Number of x bins: "<<hist->GetNbinsX()<<std::endl;
209  for(int ixbin = 1; ixbin <= hist->GetNbinsX(); ++ixbin){
210  // get the projection of the 2d hist
211  TH1D* projection = hist->ProjectionY("",ixbin,ixbin,"o");
212  //std::cout<<" X bin: "<<ixbin<<" Entries: "<<projection->GetEntries()<<std::endl;
213  if(projection->GetEntries() == 0){ continue; }
214  // fit the projection to a Gaussian around the peak
215  double mean = 0.00175 ,err = 0.001;
216  GausFit(projection,mean,err,fitmin,fitmax,lowsig,highsig);
217  //std::cout<<" Mean: "<<mean<<" Error: "<<err<<std::endl;
218  peakhist->SetBinContent(ixbin,mean);
219  peakhist->SetBinError(ixbin,err);
220 
221  }
222 
223  return peakhist;
224 
225 }
const XML_Char * name
Definition: expat.h:151
void GausFit(TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig)
void ProtonDedxCompBirk ( )

Definition at line 228 of file ProtonDedxCompBirk.C.

References demo5::c2, chisquared::c3, chisquared::c4, chi2(), om::cout, release_diff::diff, dist, allTimeWatchdog::endl, demo5::expt, make_syst_table_plots::ibin, MakeRatio(), MakeSplitRatio(), nsteps, and PeakProjectionY().

229 {
230 
231  //const char* filenamedata = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/datadedxVxFineS15-05-22.root";
232  //const char* filenamemc = "/local/nova/users/raddatz/S15-05-04_reco_pid_tuning/ReMId/DataMCDedx/mcdedxVxFineS15-05-22.root";
233 
234  const char* filenamedata = "/local/nova/users/lein/datadedxVxFineFAProton.root";
235  const char* filenamemc = "/local/nova/users/lein/mcdedxVxFineFAProton.root";
236  //const char* filenamemcsig = "/local/nova/users/lein/mcSignaldedxVxFineFAProton.root";
237  const char* filenamemcbirkB = "/local/nova/users/lein/mcdedxVxFineFAProtonBirksB.root";
238  const char* filenamemcbirkC = "/local/nova/users/lein/mcdedxVxFineFAProtonBirksC.root";
239 
240  TFile* filedata = new TFile(filenamedata,"READ");
241  TFile* filemc = new TFile(filenamemc,"READ");
242  //TFile* filemcsig = new TFile(filenamemcsig,"READ");
243  TFile* filemcbirkB = new TFile(filenamemcbirkB,"READ");
244  TFile* filemcbirkC = new TFile(filenamemcbirkC,"READ");
245 
246 
247  TH2F* dedxdata = (TH2F*)filedata->Get("DedxVx");
248  TH2F* dedxmc = (TH2F*)filemc->Get("DedxVx");
249  //TH2F* dedxmcsig = (TH2F*)filemcsig->Get("DedxVx");
250  TH2F* dedxmcbirkB = (TH2F*)filemcbirkB->Get("DedxVx");
251  TH2F* dedxmcbirkC = (TH2F*)filemcbirkC->Get("DedxVx");
252 
253  //std::cout<<"This is the pointer: "<<dedxdata<<std::endl;
254 
255  gStyle->SetOptStat(00000);
256 
257  TH1D* dedxpeakdata = PeakProjectionY(dedxdata,"dedxpeakdata");
258  TH1D* dedxpeakmc = PeakProjectionY(dedxmc,"dedxpeakmc");
259  //TH1D* dedxpeakmcsig = PeakProjectionY(dedxmcsig,"dedxpeakmcsig");
260  TH1D* dedxpeakmcbirkB = PeakProjectionY(dedxmcbirkB,"dedxpeakmcbirkB");
261  TH1D* dedxpeakmcbirkC = PeakProjectionY(dedxmcbirkC,"dedxpeakmcbirkC");
262 
263 
264  TCanvas* c0 = new TCanvas("c0","c0");
265  dedxpeakdata->Draw();
266  dedxpeakmc->Draw("same");
267  //dedxpeakmcsig->Draw("same");
268  dedxpeakmcbirkB->Draw("same");
269  dedxpeakmcbirkC->Draw("same");
270  dedxpeakmc->SetLineColor(kRed);
271  dedxpeakmc->SetMarkerColor(kRed);
272  //dedxpeakmcsig->SetLineColor(kGreen+2);
273  //dedxpeakmcsig->SetMarkerColor(kGreen+2);
274  dedxpeakmcbirkB->SetLineColor(kGreen+2);
275  dedxpeakmcbirkB->SetMarkerColor(kGreen+2);
276  dedxpeakmcbirkC->SetLineColor(kBlue-4);
277  dedxpeakmcbirkC->SetMarkerColor(kBlue-4);
278  dedxpeakdata->GetYaxis()->SetRangeUser(0.0,0.008);
279  dedxpeakdata->GetXaxis()->SetNoExponent();
280  //dedxpeakdata->GetYaxis()->SetNoExponent();
281  TLegend* leg0 = new TLegend(0.65,0.60,0.85,0.85);
282  leg0->AddEntry(dedxpeakdata,"Data","LE");
283  leg0->AddEntry(dedxpeakmc,"MC","LE");
284  //leg0->AddEntry(dedxpeakmcsig,"MC: #nu_{#mu} CC QE","LE");
285  leg0->AddEntry(dedxpeakmcbirkB,"MC BirksB","LE");
286  leg0->AddEntry(dedxpeakmcbirkC,"MC BirksC","LE");
287  leg0->Draw("same");
288  c0->Print("dedxfineProton.pdf");
289 
290  TCanvas* c00 = new TCanvas("c00","c00");
291  dedxdata->Draw("colz");
292  dedxpeakdata->Draw("same");
293  dedxdata->GetYaxis()->SetRangeUser(0.0,0.008);
294  dedxdata->GetXaxis()->SetNoExponent();
295  //dedxdata->GetYaxis()->SetNoExponent();
296  c00->Print("dedxDataMotherProton.pdf");
297 
298  TCanvas* c000 = new TCanvas("c000","c000");
299  //dedxpeakmc->SetLineColor(kBlack);
300  dedxmc->Draw("colz");
301  dedxpeakmc->Draw("same");
302  dedxmc->GetYaxis()->SetRangeUser(0.0,0.008);
303  dedxmc->GetXaxis()->SetNoExponent();
304  //dedxmc->GetYaxis()->SetNoExponent();
305  c000->Print("dedxMCMotherProton.pdf");
306  /*
307  TCanvas* c000s = new TCanvas("c000s","c000s");
308  //dedxpeakmc->SetLineColor(kBlack);
309  dedxmcsig->Draw("colz");
310  dedxpeakmcsig->Draw("same");
311  dedxmcsig->GetYaxis()->SetRangeUser(0.0,0.008);
312  dedxmcsig->GetXaxis()->SetNoExponent();
313  //dedxmc->GetYaxis()->SetNoExponent();
314  c000s->Print("dedxMCSigMotherProton.pdf");
315  */
316  TCanvas* c000b = new TCanvas("c000b","c000b");
317  //dedxpeakmc->SetLineColor(kBlack);
318  dedxmcbirkB->Draw("colz");
319  dedxpeakmcbirkB->Draw("same");
320  dedxmcbirkB->GetYaxis()->SetRangeUser(0.0,0.008);
321  dedxmcbirkB->GetXaxis()->SetNoExponent();
322  //dedxmc->GetYaxis()->SetNoExponent();
323  c000b->Print("dedxMCBirkBMotherProton.pdf");
324  TCanvas* c000c = new TCanvas("c000c","c000c");
325  //dedxpeakmc->SetLineColor(kBlack);
326  dedxmcbirkC->Draw("colz");
327  dedxpeakmcbirkC->Draw("same");
328  dedxmcbirkC->GetYaxis()->SetRangeUser(0.0,0.008);
329  dedxmcbirkC->GetXaxis()->SetNoExponent();
330  //dedxmc->GetYaxis()->SetNoExponent();
331  c000c->Print("dedxMCBirkCMotherProton.pdf");
332 
333 
334 
335  int nsteps = 1000;
336  double minscale = 0.90;
337  double maxscale = 1.05;
338 
339  TH1F* chi2match = new TH1F("chi2match",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
340  TH1F* chi2matchrestrict = new TH1F("chi2matchrestrict",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
341 
342 
343  for(int ibin = 1; ibin <= chi2match->GetNbinsX(); ++ibin){
344  // Get the scale factor
345  double scalefactor = chi2match->GetXaxis()->GetBinCenter(ibin);
346  //calculate a chi2 between data and scaled mc
347  double chi2 = 0;
348  double chi2rest = 0;
349  for(int jbin = 1; jbin <= dedxpeakdata->GetNbinsX(); ++jbin){
350  // get the dist from the end of the track
351  double dist = dedxpeakdata->GetXaxis()->GetBinCenter(jbin);
352  double expt = dedxpeakdata->GetBinContent(jbin);
353  double errdata = dedxpeakdata->GetBinError(jbin);
354  double errmc = dedxpeakmc->GetBinError(jbin);
355  double err2 = errdata*errdata+errmc*errmc;
356  double obsunscale = dedxpeakmc->GetBinContent(jbin);
357  double obs = scalefactor*dedxpeakmc->GetBinContent(jbin);
358  double diff = expt-obs;
359  //chi2+=(diff*diff/expt);
360  // remove error cases
361  if(expt < 0 || obsunscale < 0){ continue; }
362  chi2+=(diff*diff/err2);
363  // std::cout<<"dist: "<<dist<<std::endl;
364  //if(dist>100 && dist < 200){
365  if(dist>200){
366  chi2rest+=(diff*diff/err2);
367  }
368 
369  }
370 // // make a histogram to store the scaled version of mc dedxpeak histogram
371 // TString name;
372 // name+=ibin;
373 // TH1D* dedxpeakmcscale = new TH1D(name.Data(),"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
374 // dedxpeakmcscale->Add(dedxpeakmc);
375 // dedxpeakmcscale->Scale(scalefactor);
376 // dedxpeakmcscale->
377 // double chi2 = dedxpeakdata->Chi2Test(dedxpeakmcscale,"CHI2");
378 // std::cout<<"scalefactor: "<<scalefactor<<" chi2: "<<chi2<<std::endl;
379  chi2match->SetBinContent(ibin,chi2);
380  chi2matchrestrict->SetBinContent(ibin,chi2rest);
381 
382  }
383 
384  TCanvas* c00a = new TCanvas("c00a","c00a");
385  chi2match->Draw();
386  c00a->Print("chi2fine.pdf");
387 
388  double optscale = chi2match->GetXaxis()->GetBinCenter(chi2match->GetMinimumBin());
389  std::cout<<"optscale: "<<optscale<<std::endl;
390 
391  TCanvas* c000a = new TCanvas("c000a","c000a");
392  chi2matchrestrict->Draw();
393  c000a->Print("chi2restrictfine.pdf");
394 
395  double optscalerestrict = chi2matchrestrict->GetXaxis()->GetBinCenter(chi2matchrestrict->GetMinimumBin());
396  std::cout<<"Restricted optscale: "<<optscalerestrict<<std::endl;
397 
398  TH1D* dedxpeakmcscale = new TH1D("OptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
399  dedxpeakmcscale->Add(dedxpeakmc);
400  dedxpeakmcscale->Scale(optscale);
401 
402  TCanvas* c2 = new TCanvas("c2","c2");
403  dedxpeakdata->Draw();
404  dedxpeakmcscale->Draw("same");
405  dedxpeakmcscale->SetLineColor(kRed);
406  dedxpeakmcscale->SetMarkerColor(kRed);
407  dedxpeakdata->GetYaxis()->SetRangeUser(0.002,0.006);
408  dedxpeakdata->GetXaxis()->SetNoExponent();
409  //dedxpeakdata->GetYaxis()->SetNoExponent();
410  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
411  leg0->AddEntry(dedxpeakdata,"Data","LE");
412  leg0->AddEntry(dedxpeakmcscale,"Scaled MC","LE");
413  leg0->Draw("same");
414  c2->Print("dedxoptfine.pdf");
415 
416  double nonoptscale = 0.97375; //1.483/1.55;
417 
418  TH1D* dedxpeakmcnonoptscale = new TH1D("NonOptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
419  dedxpeakmcnonoptscale->Add(dedxpeakmc);
420  dedxpeakmcnonoptscale->Scale(nonoptscale);
421 
422  TCanvas* c3 = new TCanvas("c3","c3");
423  dedxpeakdata->Draw();
424  dedxpeakmcnonoptscale->Draw("same");
425  dedxpeakmcnonoptscale->SetLineColor(kRed);
426  dedxpeakmcnonoptscale->SetMarkerColor(kRed);
427  dedxpeakdata->GetYaxis()->SetRangeUser(0.002,0.006);
428  dedxpeakdata->GetXaxis()->SetNoExponent();
429  //dedxpeakdata->GetYaxis()->SetNoExponent();
430  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
431  leg0->AddEntry(dedxpeakdata,"Data","LE");
432  leg0->AddEntry(dedxpeakmcnonoptscale,"Scaled MC","LE");
433  leg0->Draw("same");
434  c3->Print("dedxnonoptfine.pdf");
435 
436 
437  TCanvas* c4 = new TCanvas("c4","c4");
438  TH1D* dedxratio = MakeRatio(dedxpeakdata,dedxpeakmc,"dedxratio");
439  //TH1D* dedxratioSig = MakeRatio(dedxpeakdata,dedxpeakmcsig,"dedxratiosig");
440  TH1D* dedxratiobirkB = MakeRatio(dedxpeakdata,dedxpeakmcbirkB,"dedxratiobirkB");
441  TH1D* dedxratiobirkC = MakeRatio(dedxpeakdata,dedxpeakmcbirkC,"dedxratiobirkC");
442  TLine* unitline = new TLine(-1,1,1001,1);
443  unitline->SetLineWidth(2);
444  unitline->SetLineStyle(7);
445  dedxratio->Draw();
446  //dedxratioSig->Draw("same");
447  dedxratiobirkB->Draw("same");
448  dedxratiobirkC->Draw("same");
449  unitline->Draw("same");
450  dedxratio->GetYaxis()->SetRangeUser(0.8,1.2);
451 
452  TCanvas* c5 = MakeSplitRatio(c0,c4,"c5");
453  c5->Print("dedxratio.pdf");
454 
455 
456 
457 
458 }
const int nsteps
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 * MakeRatio(TH1D *num, TH1D *denom, const char *name)
TCanvas * MakeSplitRatio(TCanvas *top, TCanvas *bot, const char *name)
TH1D * PeakProjectionY(TH2F *hist, const char *name, double minrange=-1.0, double maxrange=-1.0, double lowsig=0.001, double highsig=0.001)