test/Macros/CompareWithSystematicError.cc

0001 #include <iostream>
0002 #include <TH1D.h>
0003 #include <TCanvas.h>
0004 #include <TFile.h>
0005 #include <TLegend.h>
0006 #include <TGraphAsymmErrors.h>
0007
0008 void drawHisto(TString type, TFile * outputFile,
0009            const double & minX, const double & maxX,
0010            const double & minY, const double & maxY)
0011 {
0012   TFile * inputFileData = new TFile("/home/demattia/MuScleFit/Collisions/JPsi/UpTo143201/OniaSelection/ScaleFit/NoCombinedFit/AdjustedLimit/test.root", "READ");
0013   TCanvas * canvasData = (TCanvas*)inputFileData->Get("canvassigmaPtVs"+type);
0014   TGraphAsymmErrors * graphData = (TGraphAsymmErrors*)canvasData->GetPrimitive("Graph_from_sigmaPtVs"+type);
0015
0016
0017   TFile * inputFile = new TFile("test.root", "READ");
0018   TCanvas * canvas = (TCanvas*)inputFile->Get("canvassigmaPtVs"+type);
0019   TGraphAsymmErrors * graph = (TGraphAsymmErrors*)canvas->GetPrimitive("Graph_from_sigmaPtVs"+type);
0020
0021   TFile * inputFile2 = new TFile("ComparedResol.root", "READ");
0022   TCanvas * canvas2 = (TCanvas*)inputFile2->Get("resolPtVS"+type);
0023   TH1D * histo = (TH1D*)canvas2->GetPrimitive("hResolPtGenVSMu_ResoVS"+type+"_resol_after");
0024
0025 //   Double_t x[n]   = {-0.22, 0.05, 0.25, 0.35, 0.5, 0.61,0.7,0.85,0.89,0.95};
0026 //   Double_t y[n]   = {1,2.9,5.6,7.4,9,9.6,8.7,6.3,4.5,1};
0027 //   Double_t exl[n] = {.05,.1,.07,.07,.04,.05,.06,.07,.08,.05};
0028 //   Double_t eyl[n] = {.8,.7,.6,.5,.4,.4,.5,.6,.7,.8};
0029 //   Double_t exh[n] = {.02,.08,.05,.05,.03,.03,.04,.05,.06,.03};
0030 //   Double_t eyh[n] = {.6,.5,.4,.3,.2,.2,.3,.4,.5,.6};
0031 //   gr = new TGraphAsymmErrors(n,x,y,exl,exh,eyl,eyh);
0032
0033   int N = graph->GetN();
0034   Double_t * x = new Double_t[N];
0035   Double_t * y = new Double_t[N];
0036   Double_t * exl = new Double_t[N];
0037   Double_t * eyl = new Double_t[N];
0038   Double_t * exh = new Double_t[N];
0039   Double_t * eyh = new Double_t[N];
0040
0041   Double_t * exlData = new Double_t[N];
0042   Double_t * eylData = new Double_t[N];
0043   Double_t * exhData = new Double_t[N];
0044   Double_t * eyhData = new Double_t[N];
0045   Double_t * xData = new Double_t[N];
0046   Double_t * yData = new Double_t[N];
0047
0048   // Loop on the bins of the MC-truth histogram and compute the systematic error from the fit
0049   for( int i=0; i<N; ++i ) {
0050     double fitValue = graph->GetY()[i];
0051     double fitBinValue = graph->GetX()[i];
0052     double fitErrorPos = graph->GetEYhigh()[i];
0053     double fitErrorNeg = graph->GetEYlow()[i];
0054     int bin = histo->FindBin(fitBinValue);
0055     double binContent = histo->GetBinContent(bin);
0056     std::cout << "fitValue("<<i<<") = " << fitValue << " + " << fitErrorPos << " - " << fitErrorNeg << std::endl;
0057     std::cout << "binContent("<<i<<") = " << binContent << std::endl;
0058     double delta = fitValue - binContent;
0059     std::cout << "diff("<<i<<") = " << delta << std::endl;
0060
0061     x[i] = fitBinValue;
0062     y[i] = fitValue;
0063     exl[i] = graph->GetEXlow()[i];
0064     eyl[i] = fitErrorNeg;
0065     exh[i] = graph->GetEXhigh()[i];
0066     eyh[i] = fitErrorPos;
0067
0068     xData[i] = graphData->GetX()[i];
0069     yData[i] = graphData->GetY()[i];
0070     exlData[i] = graphData->GetEXlow()[i];
0071     eylData[i] = graphData->GetEYlow()[i];
0072     exhData[i] = graphData->GetEXhigh()[i];
0073     eyhData[i] = graphData->GetEYhigh()[i];
0074
0075     if( delta > 0 ) {
0076       std::cout << "before eyl["<<i<<"] = " << eyl[i] << std::endl;
0077       eylData[i] = sqrt(fitErrorNeg*fitErrorNeg + delta*delta);
0078       std::cout << "after eyl["<<i<<"] = " << eyl[i] << std::endl;
0079     }
0080     else {
0081       std::cout << "before eyh["<<i<<"] = " << eyh[i] << std::endl;
0082       eyhData[i] = sqrt(fitErrorPos*fitErrorPos + delta*delta);
0083       std::cout << "after eyh["<<i<<"] = " << eyh[i] << std::endl;
0084     }
0085   }
0086
0087   TGraphAsymmErrors * newGraph = new TGraphAsymmErrors(N,x,y,exl,exh,eyl,eyh);
0088   TGraphAsymmErrors * newGraphData = new TGraphAsymmErrors(N,xData,yData,exlData,exhData,eylData,eyhData);
0089   TGraph * lineGraph = new TGraph(N, xData, yData);
0090   newGraph->SetTitle("TGraphAsymmErrors Example");
0091   // TGraphAsymmErrors * newGraph = graph;
0092
0093   newGraph->GetXaxis()->SetRangeUser(minX, maxX);
0094   newGraph->GetYaxis()->SetRangeUser(minY, maxY);
0095   // canvas->Draw();
0096   TCanvas * newCanvas = new TCanvas("newCanvas", "newCanvas", 1000, 800);
0097   newCanvas->Draw();
0098   newGraphData->SetFillColor(kGray);
0099   newGraphData->GetXaxis()->SetRangeUser(-2.4,2.4);
0100   //  gStyle->SetOptTitle(0);
0101   newGraphData->Draw("A2");
0102   newGraph->Draw("P");
0103   // histo->Draw("SAME");
0104   // graph->Draw("SAME");
0105   lineGraph->Draw("SAME");
0106
0107   TLegend * legend = new TLegend(0.7,0.71,0.98,1.);
0108   legend->SetTextSize(0.02);
0109   legend->SetFillColor(0); // Have a white background
0110   legend->AddEntry(newGraph, "resolution from fit on MC", "lep");
0111   legend->AddEntry(newGraphData, "resolution from fit on Data");
0112
0113   legend->Draw("SAME");
0114
0115   outputFile->cd();
0116   canvas->Write();
0117 }
0118
0119 void CompareWithSystematicError()
0120 {
0121   TFile * outputFile = new TFile("output.root", "RECREATE");
0122   drawHisto("Pt", outputFile, 0., 20., 0., 0.06);
0123   drawHisto("Eta", outputFile, -2.39, 2.39, 0., 0.06);
0124   outputFile->Write();
0125   outputFile->Close();
0126 }