MCJet/test/DrawL3.C

0001 #include "Settings.h"
0002 void DrawL3()
0003 {
0004   gROOT->SetStyle("Plain");
0005   gStyle->SetOptStat(0000);
0006   gStyle->SetOptFit(000);
0007   gStyle->SetPalette(1);
0008
0009   char name[100];
0010   int i;
0011   double x,y,e;
0012   TFile *inf = new TFile(L3OutputROOTFilename,"r");
0013   TGraphErrors *g_Cor, *g_Resp;
0014   TF1 *CorFit, *RespFit;
0015   TMatrixD *COV_Cor, *COV_Resp;
0016   TH1F *hCorUncertainty, *hRespUncertainty, *hCorFrac, *hRespFrac;
0017   TPaveText *pave = new TPaveText(0.3,0.7,0.5,0.85,"NDC");
0018   pave->AddText(Version);
0019   pave->AddText(Algorithm);
0020   pave->SetLineColor(0);
0021   pave->SetBorderSize(0);
0022   pave->SetFillColor(0);
0023   pave->SetBorderSize(0);
0024   ///////////////////////////////////////////////////////////////
0025   g_Cor = (TGraphErrors*)inf->Get("Correction_vs_CaloPt");
0026   COV_Cor = (TMatrixD*)inf->Get("CovMatrix_Correction");
0027   CorFit = (TF1*)g_Cor->GetFunction("CorFit");
0028   CorFit->SetRange(1,5000);
0029   g_Resp = (TGraphErrors*)inf->Get("Response_vs_RefPt");
0030   COV_Resp = (TMatrixD*)inf->Get("CovMatrix_Resp");
0031   RespFit = (TF1*)g_Resp->GetFunction("RespFit");
0032   RespFit->SetRange(5,5000);
0033   hCorUncertainty = new TH1F("CorrectionUncertainty","CorrectionUncertainty",5000,1,5001);
0034   hRespUncertainty = new TH1F("ResponseUncertainty","ResponseUncertainty",5000,5,5005);
0035   hCorFrac = new TH1F("FractionalCorrectionUncertainty","FractionalCorrectionUncertainty",5000,1,5001);
0036   hRespFrac = new TH1F("FractionalResponseUncertainty","FractionalResponseUncertainty",5000,5,5005);
0037   for(i=0;i<5000;i++)
0038     {
0039       x = hCorUncertainty->GetBinCenter(i+1);
0040       y = CorFit->Eval(x);
0041       e = FitUncertainty(false,CorFit,COV_Cor,x);
0042       hCorUncertainty->SetBinContent(i+1,y);
0043       hCorUncertainty->SetBinError(i+1,e);
0044       hCorFrac->SetBinContent(i+1,100*e/y);
0045
0046       x = hRespUncertainty->GetBinCenter(i+1);
0047       y = RespFit->Eval(x);
0048       e = FitUncertainty(true,RespFit,COV_Resp,x);
0049       hRespUncertainty->SetBinContent(i+1,y);
0050       hRespUncertainty->SetBinError(i+1,e);
0051       hRespFrac->SetBinContent(i+1,100*e/y);
0052     }
0053   ////////////////////// Correction ///////////////////////////////////////
0054   TCanvas *c_Correction = new TCanvas("Correction","Correction",900,600);
0055   c_Correction->cd();
0056   gPad->SetLogx();
0057   hCorUncertainty->SetTitle("");
0058   g_Cor->SetMarkerStyle(20);
0059   g_Cor->SetMarkerColor(1);
0060   g_Cor->SetLineColor(1);
0061   hCorUncertainty->SetMaximum(4);
0062   hCorUncertainty->SetMinimum(1);
0063   hCorUncertainty->GetXaxis()->SetTitle("p_{T} (GeV)");
0064   hCorUncertainty->GetYaxis()->SetTitle("L3Correction factor");
0065   CorFit->SetLineColor(2);
0066   CorFit->SetLineWidth(2);
0067   CorFit->SetParNames("a0","a1","a2","a3","a4");
0068   hCorUncertainty->SetLineColor(5);
0069   hCorUncertainty->SetFillColor(5);
0070   hCorUncertainty->SetMarkerColor(5);
0071   hCorUncertainty->Draw("E3");
0072   g_Cor->Draw("Psame");
0073   pave->Draw();
0074   TLegend *leg = new TLegend(0.6,0.65,0.89,0.89);
0075   leg->AddEntry(g_Cor,"measurement","LP");
0076   leg->AddEntry(CorFit,"fit","L");
0077   leg->AddEntry(hCorUncertainty,"fit uncertainty","F");
0078   leg->SetFillColor(0);
0079   leg->SetLineColor(0);
0080   leg->Draw();
0081   ////////////////////// Response ///////////////////////////////////////
0082   TCanvas *c_Response = new TCanvas("Response","Response",900,600);
0083   c_Response->cd();
0084   gPad->SetLogx();
0085   hRespUncertainty->SetTitle("");
0086   g_Resp->SetMarkerStyle(20);
0087   g_Resp->SetMarkerColor(1);
0088   g_Resp->SetLineColor(1);
0089   hRespUncertainty->SetMaximum(1);
0090   hRespUncertainty->SetMinimum(0);
0091   hRespUncertainty->GetXaxis()->SetTitle("p_{T}^{gen} (GeV)");
0092   hRespUncertainty->GetYaxis()->SetTitle("Response");
0093   RespFit->SetLineColor(2);
0094   hRespUncertainty->SetLineColor(5);
0095   hRespUncertainty->SetFillColor(5);
0096   hRespUncertainty->SetMarkerColor(5);
0097   hRespUncertainty->Draw("E3");
0098   g_Resp->Draw("Psame");
0099   pave->Draw();
0100   TLegend *leg = new TLegend(0.6,0.15,0.89,0.39);
0101   leg->AddEntry(g_Resp,"measurement","LP");
0102   leg->AddEntry(RespFit,"fit","L");
0103   leg->AddEntry(hRespUncertainty,"fit uncertainty","F");
0104   leg->SetFillColor(0);
0105   leg->SetLineColor(0);
0106   leg->Draw();
0107   ////////////////////// Correction - Response closure ///////////////////////////////////////
0108   TH1F *hClosure = new TH1F("hClosure","hClosure",1000,5,5005);
0109   for(int i=0;i<1000;i++)
0110     {
0111       double dx = 5;
0112       double x = 5+dx*i;
0113       double y = Closure(x,CorFit,RespFit);
0114       hClosure->SetBinContent(i+1,y);
0115       hClosure->SetBinError(i+1,0.);
0116     }
0117   TCanvas *can = new TCanvas("Closure","Closure",900,600);
0118   gPad->SetLogx();
0119   hClosure->SetTitle("");
0120   hClosure->GetXaxis()->SetTitle("p_{T}");
0121   hClosure->GetYaxis()->SetTitle("C(p_{T})#times R(p_{T}C(p_{T}))");
0122   hClosure->Draw();
0123   ////////////////////// Fractional Correction Fit Uncertainty ///////////////////////////////////////
0124   TCanvas *can = new TCanvas("FracCorrUnc","FracCorrUnc",900,600);
0125   gPad->SetLogx();
0126   hCorFrac->SetTitle("");
0127   hCorFrac->GetXaxis()->SetTitle("p_{T} (GeV)");
0128   hCorFrac->GetYaxis()->SetTitle("Fractional Correction Fitting Uncertainty (%)");
0129   hCorFrac->Draw();
0130   ////////////////////// Fractional Response Fit Uncertainty ///////////////////////////////////////
0131   TCanvas *can = new TCanvas("FracRespUnc","FracRespUnc",900,600);
0132   gPad->SetLogx();
0133   hRespFrac->SetTitle("");
0134   hRespFrac->GetXaxis()->SetTitle("p_{T}^{gen} (GeV)");
0135   hRespFrac->GetYaxis()->SetTitle("Fractional Response Fit Uncertainty (%)");
0136   hRespFrac->Draw();
0137 }
0138 ///////////////////////////////////////////////////////////////////////
0139 double FitUncertainty(bool IsResponse, TF1* f, TMatrixD* COV, double x)
0140 {
0141   int i,j,dim,N,npar;
0142   double df,sum,y,z,x;
0143   double PartialDerivative[10],Parameter[10];
0144   if (IsResponse)
0145     npar = 5;
0146   else
0147     npar = 4;
0148   N = f->GetNumberFreeParameters();
0149   dim = COV->GetNrows();
0150   if (dim != npar || N != npar)
0151     {
0152       cout<<"ERROR: wrong number of parameters !!!!"<<endl;
0153       return(-1);
0154     }
0155   for(i=0;i<npar;i++)
0156     Parameter[i] = f->GetParameter(i);
0157   z = pow(log10(x),Parameter[2]);
0158   PartialDerivative[0] = 1.;
0159   PartialDerivative[1] = 1./(z+Parameter[3]);
0160   PartialDerivative[3] = -Parameter[1]/pow(z+Parameter[3],2);
0161   PartialDerivative[2] = PartialDerivative[3]*log(log10(x))*z;
0162   if (IsResponse)
0163     {
0164       PartialDerivative[1] = -1./(z+Parameter[3]);
0165       PartialDerivative[3] = Parameter[1]/pow(z+Parameter[3],2);
0166       PartialDerivative[4] = 1./x;
0167     }
0168   sum = 0.;
0169   for(i=0;i<npar;i++)
0170     for(j=0;j<npar;j++)
0171       {
0172         y = PartialDerivative[i]*PartialDerivative[j]*COV(i,j);
0173         sum+=y;
0174       }
0175   df = sqrt(sum);
0176   return df;
0177 }
0178 double Closure(double x, TF1 *f1, TF1 *f2)
0179 {
0180   double y1,y,tmp;
0181   y1 = f1->Eval(x);
0182   y = x*y1;
0183   tmp = y1*f2->Eval(y);
0184   return tmp;
0185 }
0186
0187