Macros/RooFit/CompareWithExpected.cc

0001 /**
0002  * This macro can be used to compare the BiasCheck_0 peak of a CB fit (produced by CompareBias.cc)
0003  * with the expected peak from the model convoluted with a gaussian.
0004  */
0005 {
0006   inputFile = new TFile("BiasCheck_0.root", "READ");
0007   inputFile->cd("MassVsEta");
0008   // gDirectory->ls();
0009   sigmaCanvas = (TCanvas*)gDirectory->Get("sigmaCanvas");
0010   sigmaHisto = (TH1D*)sigmaCanvas->GetPrimitive("sigmaHisto");
0011
0012   meanCanvas = (TCanvas*)gDirectory->Get("meanCanvas");
0013   meanHisto = (TH1D*)meanCanvas->GetPrimitive("meanHisto");
0014
0015   func = new TF1("linear", "[0] + [1]*x", 20., 82.);
0016   func->SetParameter(0, 3.09625);
0017   func->SetParameter(1, -0.0000380253);
0018
0019   expectedMeanHisto = (TH1D*)meanHisto->Clone("expectedMeanHisto");
0020
0021   std::cout << "NbinsX = " << meanHisto->GetNbinsX() << std::endl;
0022   for( int i=1; i<=meanHisto->GetNbinsX(); ++i ) {
0023     double sigma = sigmaHisto->GetBinContent(i);
0024     expectedMeanHisto->SetBinContent(i, 0);
0025     if( sigma != 0 ) {
0026       std::cout << "mean("<<i<<") = " << meanHisto->GetBinContent(i) << std::endl;
0027       // std::cout << "sigma("<<i<<") = " << sigma << std::endl;
0028       std::cout << "expected mean("<<i<<") = " << func->Eval(sigma*1000) << std::endl;
0029       expectedMeanHisto->SetBinContent(i, func->Eval(sigma*1000));
0030     }
0031   }
0032
0033   newCanvas = new TCanvas("newCanvas", "newCanvas", 1000, 800);
0034   newCanvas->cd();
0035   newCanvas->Draw();
0036   meanHisto->Draw();
0037   expectedMeanHisto->SetMarkerColor(kRed);
0038   expectedMeanHisto->SetLineColor(kRed);
0039   expectedMeanHisto->Draw("same");
0040   // sigmaHisto->Draw();
0041   leg = new TLegend(0.1,0.7,0.48,0.9);
0042   leg->SetHeader("CrystalBall fit peak");
0043   leg->AddEntry(meanHisto, "Reconstructed MC", "l");
0044   leg->AddEntry(expectedMeanHisto, "Expected from model", "l");
0045   leg->Draw("same");
0046
0047   meanHisto->GetXaxis()->SetRangeUser(-2.4, 2.4);
0048   meanHisto->GetYaxis()->SetRangeUser(3.092, 3.1);
0049   meanHisto->GetXaxis()->SetTitle("muon #eta");
0050   meanHisto->GetYaxis()->SetTitle("peak (GeV)");
0051 }