Calibrations/MC-Calibrations/Jet_Resp_respFitter.C

0001 #include "sPhenixStyle.h"
0002 #include "sPhenixStyle.C"
0003 #include "respFitterDef.h"
0004 // This macro analyzes output from the jet validation sub module of the
0005 // JS-Jet module. To first order, it simply computes the jet energy scale (JES)
0006 // and jet energy resolution (JER), but also compute the jet energy linearity
0007 // necessary for the numerical inversion which is used to calibrate the MC
0008 // to itself.
0009
0010 void Jet_Resp_respFitter(const char* fin = "", int isLin = 0)
0011 {
0012   //old definitions we need to make things work
0013   vector<float> etaBins;
0014   for(float i = etaStart; i <= etaEnd; i+=etaGap)
0015     {
0016       etaBins.push_back(i);
0017     }
0018   const int nEtaBins = etaBins.size() - 1 + extraBins; //one inclusive bin
0019
0020
0021   gSystem -> CopyFile(fin,Form("calibHists_%s",fin),kTRUE);
0022   TFile *f = new TFile(Form("calibHists_%s",fin),"UPDATE");
0023   TH2D *h_MC_Reso_pt[nEtaBins];
0024   for(int i = 0; i < nEtaBins; i++)
0025     {
0026       h_MC_Reso_pt[i] = (TH2D*)f -> Get(Form("h_MC_Reso_eta%d",i));
0027     }
0028
0029
0030   Float_t pt_range_truth[pt_N_truth];
0031   for(int i = 0; i < pt_N_truth+1; i++)
0032   {
0033     pt_range_truth[i] = 5+75./pt_N_truth * i;
0034   }
0035
0036   TLegend *leg = new TLegend(.25,.2,.6,.5);
0037   leg->SetFillStyle(0);
0038   gStyle->SetOptFit(1);
0039
0040   TGraphErrors *g_jes[nEtaBins];
0041   TGraphErrors *g_jer[nEtaBins];
0042   for(int i = 0; i < nEtaBins; i++)
0043     {
0044       g_jes[i] = new TGraphErrors(pt_N);
0045       g_jer[i] = new TGraphErrors(pt_N);
0046     }
0047
0048   int stop;//Sets the number of bins to look at, which are different depending on whether or not you're doing the calibration
0049   float fitEnd;//The fit is super generic so it really doesn't care if you're fitting from 0-2 (JES) or 0-80 (linearity), but just to be thorough
0050   float fitStart;
0051   f->cd();
0052   if(isLin)
0053     {
0054       stop = pt_N_truth-1;
0055       fitEnd = 80;
0056       fitStart = 1.;
0057     }
0058   else
0059     {
0060       stop = pt_N;
0061       fitEnd = 2.;
0062       fitStart = 0.1;
0063     }
0064   for (int i = 0; i < stop; i++)
0065     {
0066       for(int j = 0; j < nEtaBins; j++)
0067     {
0068
0069       TF1 *func = new TF1("func","gaus",0,fitEnd);
0070       //h_MC_Reso_pt[j]->GetXaxis()->SetRange(i+1,i+1);
0071       TH1F *h_temp = (TH1F*) h_MC_Reso_pt[j]->ProjectionY("proj",i+1,i+1);
0072
0073       func -> SetParameter(1, h_temp -> GetBinCenter(h_temp -> GetMaximumBin()));
0074       h_temp->Fit(func,"Lm","",fitStart,fitEnd);
0075       if(isLin)h_temp -> SetTitle(Form("%g < p_{T}^{Truth} < %g",pt_range_truth[i],pt_range_truth[i+1]));
0076       else h_temp -> SetTitle(Form("%g < p_{T}^{Truth} < %g", pt_range[i],  pt_range[i+1]));
0077       h_temp -> Write(Form("hFit_pt%d_eta%d_Fit0",i,j));
0078       float meanGuess = func -> GetParameter(1);
0079       float rangeCut = 1.0;
0080       if(isLin) rangeCut = nSigmaLin;
0081       else rangeCut = nSigmaRat;
0082       func -> SetParameter(1, h_temp -> GetBinCenter(h_temp -> GetMaximumBin()));
0083       func -> SetParameter(2,1.0*meanGuess);
0084       h_temp->Fit(func,"","",func->GetParameter(1)-rangeCut*func->GetParameter(2),func->GetParameter(1)+rangeCut*func->GetParameter(2));
0085       func->SetLineColor(kRed);
0086
0087       h_temp -> Write(Form("hFit_pt%d_eta%d_Fit1",i,j));
0088
0089       leg->Clear();
0090       float dsigma = func->GetParError(2);
0091       float denergy = func->GetParError(1);
0092       float sigma = func->GetParameter(2);
0093       float energy = func->GetParameter(1);
0094
0095       float djer = dsigma/energy + sigma*denergy/pow(energy,2);//correct way to compute error on the resolution.
0096       if(!isLin)
0097         {
0098           g_jes[j]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(1));
0099           g_jes[j]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),func->GetParError(1));
0100
0101           g_jer[j]->SetPoint(i,0.5*(pt_range[i]+pt_range[i+1]),func->GetParameter(2)/func->GetParameter(1));
0102           g_jer[j]->SetPointError(i,0.5*(pt_range[i+1]-pt_range[i]),djer);
0103         }
0104       else
0105         {
0106           g_jes[j]->SetPoint(i,0.5*(pt_range_truth[i]+pt_range_truth[i+1]),func->GetParameter(1));
0107           g_jes[j]->SetPointError(i,0.5*(pt_range_truth[i+1]-pt_range_truth[i]),func->GetParError(1));
0108           g_jer[j]->SetPoint(i,0.5*(pt_range_truth[i]+pt_range_truth[i+1]),func->GetParameter(2)/func->GetParameter(1));
0109           g_jer[j]->SetPointError(i,0.5*(pt_range_truth[i+1]-pt_range_truth[i]),func->GetParError(2));
0110         }
0111
0112     }
0113     }
0114   for(int i = 0; i < nEtaBins; i++)
0115     {
0116       g_jes[i]->Write(Form("g_jes_eta%d",i));
0117       g_jer[i]->Write(Form("g_jer_eta%d",i));
0118     }
0119   std::cout << "All done!" << std::endl;
0120 }
0121
0122