24 TH1D*
ratio =
new TH1D(name,
";;Data/MC",num->GetNbinsX(),num->GetXaxis()->GetXbins()->GetArray());
26 ratio->GetXaxis()->SetTitle(num->GetXaxis()->GetTitle());
28 ratio->Divide(num,denom);
30 ratio->SetLineColor(denom->GetLineColor());
31 ratio->SetMarkerColor(denom->GetMarkerColor());
32 ratio->SetMarkerSize(denom->GetMarkerSize());
42 TCanvas*
can =
new TCanvas(name,name,600,550);
48 can->SetBottomMargin(0);
50 TVirtualPad* can_1 = can->cd(1);
51 can_1->SetPad(0,0.4,1,1);
52 can_1->SetFillStyle(0);
54 can_1->SetRightMargin(0.05);
55 can_1->SetTopMargin(0.1);
56 can_1->SetBottomMargin(0.005);
57 can_1->SetLeftMargin(0.13);
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());
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);
80 TVirtualPad* can_2 = can->cd(2);
81 can_2->SetPad(0,0,1,0.4);
83 can_2->SetTopMargin(0.025);
84 can_2->SetBottomMargin(0.22);
85 can_2->SetRightMargin(0.05);
86 can_2->SetLeftMargin(0.13);
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());
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);
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);
135 double lowCut = integral*percentCut;
136 double highCut = integral*(1-percentCut);
140 bool foundLow =
false;
141 double ongoingIntegral = 0.0;
144 double bincont = hist->GetBinContent(
ibin);
147 ongoingIntegral += bincont;
149 if (ongoingIntegral > lowCut){
156 if ((ongoingIntegral<highCut)&&((ongoingIntegral+bincont)>=highCut)){
157 ongoingIntegral += bincont;
161 else if (ongoingIntegral<highCut){
162 ongoingIntegral += bincont;
167 low = hist->GetXaxis()->GetBinUpEdge(lowbin);
168 high = hist->GetXaxis()->GetBinLowEdge(highbin);
176 double modex = hist->GetXaxis()->GetBinCenter(hist->GetMaximumBin());
187 if( mode > maxx || mode < minx){
190 int minxbin = hist->GetXaxis()->FindFixBin(minx);
191 int maxxbin = hist->GetXaxis()->FindFixBin(maxx);
192 double modeval = hist->GetBinContent(minxbin);
193 int modebin = minxbin;
196 double bincont = hist->GetBinContent(
ibin);
198 if(bincont > modeval){
204 mode = hist->GetXaxis()->GetBinCenter(modebin);
212 void GausFit(TH1D*
hist,
double&
mean,
double& err,
double minrange,
double maxrange,
double lowsig,
double highsig,
bool useTruncate)
219 double fitminx = -1.0;
220 double fitmaxx = -1.0;
226 fitminx = mode-lowsig;
227 fitmaxx = mode+highsig;
230 if(fitminx < minrange){ fitminx = minrange; }
231 if(fitmaxx > maxrange){ fitmaxx = maxrange; }
233 TF1* fgaus =
new TF1(
"gaus",
"gaus(0)",fitminx,fitmaxx);
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);
243 if (integral == 320){
244 TCanvas*
cs =
new TCanvas(
"cs",
"cs");
245 hist->Fit(fgaus,
"RIQ");
246 cs->Print(
"test.pdf");
248 int fitstatus = hist->Fit(fgaus,
"RI");
251 mean = fgaus->GetParameter(1);
252 err = fgaus->GetParError(1);
253 if(mean > maxrange || mean < minrange){
265 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)
268 TH1D* peakhist =
new TH1D(name,
"",hist->GetNbinsX(),hist->GetXaxis()->GetXbins()->GetArray());
269 peakhist->GetYaxis()->SetTitle(hist->GetYaxis()->GetTitle());
270 peakhist->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
272 double fitmin = minrange;
273 double fitmax = maxrange;
275 if(minrange == -1.0){ fitmin = hist->GetXaxis()->GetXmin(); }
276 if(maxrange == -1.0){ fitmax = hist->GetXaxis()->GetXmax(); }
279 for(
int ixbin = 1; ixbin <= hist->GetNbinsX(); ++ixbin){
281 TH1D* projection = hist->ProjectionY(
"",ixbin,ixbin,
"o");
283 if(projection->GetEntries() == 0){
continue; }
285 double mean = 0.00175 ,err = 0.001;
286 double lowTruncate = -1.0;
287 double highTruncate = -1.0;
290 if ((lowTruncate != -1.0) && (highTruncate != -1.0)){
291 lowsig = lowTruncate;
292 highsig = highTruncate;
296 GausFit(projection,mean,err,fitmin,fitmax,lowsig,highsig,useTruncate);
298 peakhist->SetBinContent(ixbin,mean);
299 peakhist->SetBinError(ixbin,err);
314 const char* filenamedata =
"/local/nova/users/lein/datadedxVgammaFineS15-05-22.root";
315 const char* filenamemc =
"/local/nova/users/lein/mcdedxVgammaFineS15-05-22.root";
317 const char* filenamedataP =
"/local/nova/users/lein/datadedxVgammaFineS15-05-22Proton.root";
318 const char* filenamemcP =
"/local/nova/users/lein/mcdedxVgammaFineS15-05-22Proton.root";
320 TFile* filedata =
new TFile(filenamedata,
"READ");
321 TFile* filemc =
new TFile(filenamemc,
"READ");
323 TFile* filedataP =
new TFile(filenamedataP,
"READ");
324 TFile* filemcP =
new TFile(filenamemcP,
"READ");
326 TH2F* dedxdata = (TH2F*)filedata->Get(
"DedxVx");
327 TH2F* dedxmc = (TH2F*)filemc->Get(
"DedxVx");
329 TH2F* dedxdataP = (TH2F*)filedataP->Get(
"DedxVx");
330 TH2F* dedxmcP = (TH2F*)filemcP->Get(
"DedxVx");
334 gStyle->SetOptStat(00000);
349 TCanvas* c0 =
new TCanvas(
"c0",
"c0");
350 dedxpeakdata->Draw();
351 dedxpeakmc->Draw(
"same");
352 dedxpeakmcP->Draw(
"same");
353 dedxpeakdataP->Draw(
"same");
354 dedxpeakmc->SetLineColor(
kRed);
355 dedxpeakmc->SetMarkerColor(
kRed);
356 dedxpeakmcP->SetLineColor(
kViolet);
357 dedxpeakmcP->SetMarkerColor(
kViolet);
358 dedxpeakdataP->SetLineColor(
kGreen+2);
359 dedxpeakdataP->SetMarkerColor(
kGreen+2);
360 dedxpeakdata->GetYaxis()->SetRangeUser(0.0,0.006);
361 dedxpeakdata->GetXaxis()->SetNoExponent();
363 TLegend* leg0 =
new TLegend(0.65,0.70,0.85,0.85);
364 leg0->AddEntry(dedxpeakdata,
"Muon Data",
"LE");
365 leg0->AddEntry(dedxpeakmc,
"Muon MC",
"LE");
366 leg0->AddEntry(dedxpeakdataP,
"Proton Data",
"LE");
367 leg0->AddEntry(dedxpeakmcP,
"Proton MC",
"LE");
369 c0->Print(
"dedxfine.pdf");
370 c0->Print(
"dedxfine.root");
372 TCanvas* c00 =
new TCanvas(
"c00",
"c00");
373 dedxdata->Draw(
"colz");
374 dedxpeakdata->Draw(
"same");
375 dedxdata->GetYaxis()->SetRangeUser(0.0,0.006);
376 dedxdata->GetXaxis()->SetNoExponent();
378 c00->Print(
"dedxDataMother.pdf");
380 TCanvas* c000 =
new TCanvas(
"c000",
"c000");
382 dedxmc->Draw(
"colz");
383 dedxpeakmc->Draw(
"same");
384 dedxmc->GetYaxis()->SetRangeUser(0.0,0.006);
385 dedxmc->GetXaxis()->SetNoExponent();
387 c000->Print(
"dedxMCMother.pdf");
389 TCanvas* c00P =
new TCanvas(
"c00P",
"c00P");
390 dedxdataP->Draw(
"colz");
391 dedxpeakdataP->Draw(
"same");
392 dedxdataP->GetYaxis()->SetRangeUser(0.0,0.006);
393 dedxdataP->GetXaxis()->SetNoExponent();
395 c00P->Print(
"dedxDataProtonMother.pdf");
397 TCanvas* c000P =
new TCanvas(
"c000P",
"c000P");
399 dedxmcP->Draw(
"colz");
400 dedxpeakmcP->Draw(
"same");
401 dedxmcP->GetYaxis()->SetRangeUser(0.0,0.006);
402 dedxmcP->GetXaxis()->SetNoExponent();
404 c000P->Print(
"dedxMCProtonMother.pdf");
408 double minscale = 0.95;
409 double maxscale = 1.05;
411 TH1F* chi2match =
new TH1F(
"chi2match",
";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
412 TH1F* chi2matchrestrict =
new TH1F(
"chi2matchrestrict",
";Scale Factor;#chi^{2}",nsteps,minscale,maxscale);
415 for(
int ibin = 1;
ibin <= chi2match->GetNbinsX(); ++
ibin){
417 double scalefactor = chi2match->GetXaxis()->GetBinCenter(
ibin);
421 for(
int jbin = 1; jbin <= dedxpeakdata->GetNbinsX(); ++jbin){
423 double dist = dedxpeakdata->GetXaxis()->GetBinCenter(jbin);
424 double expt = dedxpeakdata->GetBinContent(jbin);
425 double errdata = dedxpeakdata->GetBinError(jbin);
426 double errmc = dedxpeakmc->GetBinError(jbin);
427 double err2 = errdata*errdata+errmc*errmc;
428 double obsunscale = dedxpeakmc->GetBinContent(jbin);
429 double obs = scalefactor*dedxpeakmc->GetBinContent(jbin);
430 double diff = expt-obs;
433 if(expt < 0 || obsunscale < 0){
continue; }
434 chi2+=(diff*diff/err2);
436 if(dist>100 && dist < 200){
437 chi2rest+=(diff*diff/err2);
450 chi2match->SetBinContent(
ibin,chi2);
451 chi2matchrestrict->SetBinContent(
ibin,chi2rest);
455 TCanvas* c00a =
new TCanvas(
"c00a",
"c00a");
457 c00a->Print(
"chi2fine.pdf");
459 double optscale = chi2match->GetXaxis()->GetBinCenter(chi2match->GetMinimumBin());
462 TCanvas* c000a =
new TCanvas(
"c000a",
"c000a");
463 chi2matchrestrict->Draw();
464 c000a->Print(
"chi2restrictfine.pdf");
466 double optscalerestrict = chi2matchrestrict->GetXaxis()->GetBinCenter(chi2matchrestrict->GetMinimumBin());
469 TH1D* dedxpeakmcscale =
new TH1D(
"OptimalMatch",
"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
470 dedxpeakmcscale->Add(dedxpeakmc);
471 dedxpeakmcscale->Scale(optscale);
473 TCanvas*
c2 =
new TCanvas(
"c2",
"c2");
474 dedxpeakdata->Draw();
475 dedxpeakmcscale->Draw(
"same");
476 dedxpeakmcscale->SetLineColor(
kRed);
477 dedxpeakmcscale->SetMarkerColor(
kRed);
478 dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
479 dedxpeakdata->GetXaxis()->SetNoExponent();
481 TLegend* leg0 =
new TLegend(0.65,0.70,0.85,0.85);
482 leg0->AddEntry(dedxpeakdata,
"Data",
"LE");
483 leg0->AddEntry(dedxpeakmcscale,
"MC",
"LE");
485 c2->Print(
"dedxoptfine.pdf");
487 double nonoptscale = 0.97375;
489 TH1D* dedxpeakmcnonoptscale =
new TH1D(
"NonOptimalMatch",
"",dedxpeakmc->GetNbinsX(),dedxpeakmc->GetXaxis()->GetXbins()->GetArray());
490 dedxpeakmcnonoptscale->Add(dedxpeakmc);
491 dedxpeakmcnonoptscale->Scale(nonoptscale);
493 TCanvas*
c3 =
new TCanvas(
"c3",
"c3");
494 dedxpeakdata->Draw();
495 dedxpeakmcnonoptscale->Draw(
"same");
496 dedxpeakmcnonoptscale->SetLineColor(
kRed);
497 dedxpeakmcnonoptscale->SetMarkerColor(
kRed);
498 dedxpeakdata->GetYaxis()->SetRangeUser(0.0014,0.0025);
499 dedxpeakdata->GetXaxis()->SetNoExponent();
501 TLegend* leg0 =
new TLegend(0.65,0.70,0.85,0.85);
502 leg0->AddEntry(dedxpeakdata,
"Data",
"LE");
503 leg0->AddEntry(dedxpeakmcnonoptscale,
"MC",
"LE");
505 c3->Print(
"dedxnonoptfine.pdf");
508 TCanvas*
c4 =
new TCanvas(
"c4",
"c4");
509 TH1D* dedxratio =
MakeRatio(dedxpeakdata,dedxpeakmc,
"dedxratio");
510 TLine* unitline =
new TLine(-1,1,1001,1);
511 unitline->SetLineWidth(2);
512 unitline->SetLineStyle(7);
514 unitline->Draw(
"same");
515 dedxratio->GetYaxis()->SetRangeUser(0.9,1.1);
518 c5->Print(
"dedxratio.pdf");
double ModeRestrict(TH1D *hist, double minx, double maxx)
void GausFit(TH1D *hist, double &mean, double &err, double minrange, double maxrange, double lowsig, double highsig, bool useTruncate)
TCanvas * MakeSplitRatio(TCanvas *top, TCanvas *bot, const char *name)
fvar< T > fabs(const fvar< T > &x)
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
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)
void FindTruncate(TH1D *hist, double &low, double &high)
TH1D * MakeRatio(TH1D *num, TH1D *denom, const char *name)