ProtonDedxComp.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 double Mode(TH1D* hist)
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 }
128 
129 double ModeRestrict(TH1D* hist, double minx, double maxx)
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 }
157 
158 
159 void GausFit(TH1D* hist, double& mean, double& err, double minrange, double maxrange, double lowsig, double highsig)
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 }
194 
195 TH1D* PeakProjectionY(TH2F* hist,const char* name, double minrange = -1.0, double maxrange = -1.0, double lowsig = 0.001, double highsig = 0.001)
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 }
226 
227 
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 = "datadedxVxFineFAProtonBirksC.root";
235  const char* filenamemc = "mcdedxVxFineFAProtonBirksC.root";
236  const char* filenamemcsig = "mcsignaldedxVxFineFAProtonBirksC.root";
237 
238  TFile* filedata = new TFile(filenamedata,"READ");
239  TFile* filemc = new TFile(filenamemc,"READ");
240  TFile* filemcsig = new TFile(filenamemcsig,"READ");
241 
242  TH2F* dedxdata = (TH2F*)filedata->Get("DedxVx");
243  TH2F* dedxmc = (TH2F*)filemc->Get("DedxVx");
244  TH2F* dedxmcsig = (TH2F*)filemcsig->Get("DedxVx");
245 
246  //std::cout<<"This is the pointer: "<<dedxdata<<std::endl;
247 
248  gStyle->SetOptStat(00000);
249 
250  TH1D* dedxpeakdata = PeakProjectionY(dedxdata,"dedxpeakdata");
251  TH1D* dedxpeakmc = PeakProjectionY(dedxmc,"dedxpeakmc");
252  TH1D* dedxpeakmcsig = PeakProjectionY(dedxmcsig,"dedxpeakmcsig");
253 
254 
255  TCanvas* c0 = new TCanvas("c0","c0");
256  dedxpeakdata->Draw();
257  dedxpeakmc->Draw("same");
258  dedxpeakmcsig->Draw("same");
259  dedxpeakmc->SetLineColor(kRed);
260  dedxpeakmc->SetMarkerColor(kRed);
261  dedxpeakmcsig->SetLineColor(kGreen+2);
262  dedxpeakmcsig->SetMarkerColor(kGreen+2);
263  dedxpeakdata->GetYaxis()->SetRangeUser(0.0,0.008);
264  dedxpeakdata->GetXaxis()->SetNoExponent();
265  //dedxpeakdata->GetYaxis()->SetNoExponent();
266  TLegend* leg0 = new TLegend(0.65,0.70,0.85,0.85);
267  leg0->AddEntry(dedxpeakdata,"Data","LE");
268  leg0->AddEntry(dedxpeakmc,"MC","LE");
269  leg0->AddEntry(dedxpeakmcsig,"MC: #nu_{#mu} CC QE","LE");
270  leg0->Draw("same");
271  c0->Print("dedxfineProton.pdf");
272 
273  TCanvas* c00 = new TCanvas("c00","c00");
274  dedxdata->Draw("colz");
275  dedxpeakdata->Draw("same");
276  dedxdata->GetYaxis()->SetRangeUser(0.0,0.008);
277  dedxdata->GetXaxis()->SetNoExponent();
278  //dedxdata->GetYaxis()->SetNoExponent();
279  c00->Print("dedxDataMotherProton.pdf");
280 
281  TCanvas* c000 = new TCanvas("c000","c000");
282  //dedxpeakmc->SetLineColor(kBlack);
283  dedxmc->Draw("colz");
284  dedxpeakmc->Draw("same");
285  dedxmc->GetYaxis()->SetRangeUser(0.0,0.008);
286  dedxmc->GetXaxis()->SetNoExponent();
287  //dedxmc->GetYaxis()->SetNoExponent();
288  c000->Print("dedxMCMotherProton.pdf");
289 
290  TCanvas* c000s = new TCanvas("c000s","c000s");
291  //dedxpeakmc->SetLineColor(kBlack);
292  dedxmcsig->Draw("colz");
293  dedxpeakmcsig->Draw("same");
294  dedxmcsig->GetYaxis()->SetRangeUser(0.0,0.008);
295  dedxmcsig->GetXaxis()->SetNoExponent();
296  //dedxmc->GetYaxis()->SetNoExponent();
297  c000s->Print("dedxMCSigMotherProton.pdf");
298 
299 
300  int nsteps = 1000;
301  double minscale = 0.90;
302  double maxscale = 1.05;
303 
304  TH1F* chi2match = new TH1F("chi2match",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
305  TH1F* chi2matchrestrict = new TH1F("chi2matchrestrict",";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
306 
307 
308  for(int ibin = 1; ibin <= chi2match->GetNbinsX(); ++ibin){
309  // Get the scale factor
310  double scalefactor = chi2match->GetXaxis()->GetBinCenter(ibin);
311  //calculate a chi2 between data and scaled mc
312  double chi2 = 0;
313  double chi2rest = 0;
314  for(int jbin = 1; jbin <= dedxpeakdata->GetNbinsX(); ++jbin){
315  // get the dist from the end of the track
316  double dist = dedxpeakdata->GetXaxis()->GetBinCenter(jbin);
317  double expt = dedxpeakdata->GetBinContent(jbin);
318  double errdata = dedxpeakdata->GetBinError(jbin);
319  double errmc = dedxpeakmc->GetBinError(jbin);
320  double err2 = errdata*errdata+errmc*errmc;
321  double obsunscale = dedxpeakmc->GetBinContent(jbin);
322  double obs = scalefactor*dedxpeakmc->GetBinContent(jbin);
323  double diff = expt-obs;
324  //chi2+=(diff*diff/expt);
325  // remove error cases
326  if(expt < 0 || obsunscale < 0){ continue; }
327  chi2+=(diff*diff/err2);
328  // std::cout<<"dist: "<<dist<<std::endl;
329  //if(dist>100 && dist < 200){
330  if(dist>200){
331  chi2rest+=(diff*diff/err2);
332  }
333 
334  }
335 // // make a histogram to store the scaled version of mc dedxpeak histogram
336 // TString name;
337 // name+=ibin;
338 // TH1D* dedxpeakmcscale = new TH1D(name.Data(),"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
339 // dedxpeakmcscale->Add(dedxpeakmc);
340 // dedxpeakmcscale->Scale(scalefactor);
341 // dedxpeakmcscale->
342 // double chi2 = dedxpeakdata->Chi2Test(dedxpeakmcscale,"CHI2");
343 // std::cout<<"scalefactor: "<<scalefactor<<" chi2: "<<chi2<<std::endl;
344  chi2match->SetBinContent(ibin,chi2);
345  chi2matchrestrict->SetBinContent(ibin,chi2rest);
346 
347  }
348 
349  TCanvas* c00a = new TCanvas("c00a","c00a");
350  chi2match->Draw();
351  c00a->Print("chi2fine.pdf");
352 
353  double optscale = chi2match->GetXaxis()->GetBinCenter(chi2match->GetMinimumBin());
354  std::cout<<"optscale: "<<optscale<<std::endl;
355 
356  TCanvas* c000a = new TCanvas("c000a","c000a");
357  chi2matchrestrict->Draw();
358  c000a->Print("chi2restrictfine.pdf");
359 
360  double optscalerestrict = chi2matchrestrict->GetXaxis()->GetBinCenter(chi2matchrestrict->GetMinimumBin());
361  std::cout<<"Restricted optscale: "<<optscalerestrict<<std::endl;
362 
363  TH1D* dedxpeakmcscale = new TH1D("OptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
364  dedxpeakmcscale->Add(dedxpeakmc);
365  dedxpeakmcscale->Scale(optscale);
366 
367  TCanvas* c2 = new TCanvas("c2","c2");
368  dedxpeakdata->Draw();
369  dedxpeakmcscale->Draw("same");
370  dedxpeakmcscale->SetLineColor(kRed);
371  dedxpeakmcscale->SetMarkerColor(kRed);
372  dedxpeakdata->GetYaxis()->SetRangeUser(0.002,0.006);
373  dedxpeakdata->GetXaxis()->SetNoExponent();
374  //dedxpeakdata->GetYaxis()->SetNoExponent();
375  leg0 = new TLegend(0.65,0.70,0.85,0.85);
376  leg0->AddEntry(dedxpeakdata,"Data","LE");
377  leg0->AddEntry(dedxpeakmcscale,"Scaled MC","LE");
378  leg0->Draw("same");
379  c2->Print("dedxoptfine.pdf");
380 
381  double nonoptscale = 0.97375; //1.483/1.55;
382 
383  TH1D* dedxpeakmcnonoptscale = new TH1D("NonOptimalMatch","",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
384  dedxpeakmcnonoptscale->Add(dedxpeakmc);
385  dedxpeakmcnonoptscale->Scale(nonoptscale);
386 
387  TCanvas* c3 = new TCanvas("c3","c3");
388  dedxpeakdata->Draw();
389  dedxpeakmcnonoptscale->Draw("same");
390  dedxpeakmcnonoptscale->SetLineColor(kRed);
391  dedxpeakmcnonoptscale->SetMarkerColor(kRed);
392  dedxpeakdata->GetYaxis()->SetRangeUser(0.002,0.006);
393  dedxpeakdata->GetXaxis()->SetNoExponent();
394  //dedxpeakdata->GetYaxis()->SetNoExponent();
395  leg0 = new TLegend(0.65,0.70,0.85,0.85);
396  leg0->AddEntry(dedxpeakdata,"Data","LE");
397  leg0->AddEntry(dedxpeakmcnonoptscale,"Scaled MC","LE");
398  leg0->Draw("same");
399  c3->Print("dedxnonoptfine.pdf");
400 
401 
402  TCanvas* c4 = new TCanvas("c4","c4");
403  TH1D* dedxratio = MakeRatio(dedxpeakdata,dedxpeakmc,"dedxratio");
404  TH1D* dedxratioSig = MakeRatio(dedxpeakdata,dedxpeakmcsig,"dedxratiosig");
405  TLine* unitline = new TLine(-1,1,1001,1);
406  unitline->SetLineWidth(2);
407  unitline->SetLineStyle(7);
408  dedxratio->Draw();
409  dedxratioSig->Draw("same");
410  unitline->Draw("same");
411  dedxratio->GetYaxis()->SetRangeUser(0.9,1.1);
412 
413  TCanvas* c5 = MakeSplitRatio(c0,c4,"c5");
414  c5->Print("dedxratio.pdf");
415 
416 
417 
418 
419 }
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1D * MakeRatio(TH1D *num, TH1D *denom, const char *name)
const int nsteps
TH1 * ratio(TH1 *h1, TH1 *h2)
double ModeRestrict(TH1D *hist, double minx, double maxx)
double chi2()
double dist
Definition: runWimpSim.h:113
void GausFit(TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig)
c2
Definition: demo5.py:33
TCanvas * MakeSplitRatio(TCanvas *top, TCanvas *bot, const char *name)
void ProtonDedxComp()
expt
Definition: demo5.py:34
OStream cout
Definition: OStream.cxx:6
double Mode(TH1D *hist)
int num
Definition: f2_nu.C:119
enum BeamMode kGreen
TGeoVolume * top
Definition: make_fe_box.C:9
TH1D * PeakProjectionY(TH2F *hist, const char *name, double minrange=-1.0, double maxrange=-1.0, double lowsig=0.001, double highsig=0.001)