| 17 |
|
#include "Math/GSLIntegrator.h" |
| 18 |
|
#include "Math/WrappedTF1.h" |
| 19 |
|
|
| 20 |
< |
TString fileName = "QCDbin4.root"; |
| 21 |
< |
// Histo range |
| 22 |
< |
double xMin = 0.0; |
| 23 |
< |
double xMax = 2.0; |
| 24 |
< |
// fit region weight |
| 25 |
< |
double wxmin = 2.1; |
| 26 |
< |
double wxmax = 0.8; |
| 20 |
> |
TString fileName = "skimmed/QCDbin4_005.root"; |
| 21 |
> |
// Dynamic(true) or static(false) range |
| 22 |
> |
bool dynamic = true; |
| 23 |
> |
// Dynamic fit region weight |
| 24 |
> |
double wxmin = 2.; // 1.82 for loose 2.1 for tight |
| 25 |
> |
double wxmax = 0.8; // 0.6 for loose 0.8 for tight |
| 26 |
|
// Select fitting function: |
| 27 |
|
TString s0 = "landau"; |
| 28 |
|
// Select New or Old reliso: |
| 31 |
|
void fitLandau() |
| 32 |
|
{ |
| 33 |
|
gStyle->SetOptFit(01111); |
| 34 |
< |
|
| 34 |
> |
|
| 35 |
|
TFile *file = TFile::Open(fileName); |
| 36 |
|
TTree *tree = (TTree*)file->Get("myTree"); |
| 37 |
|
TH1D *h0 = (TH1D*)file->Get("histos/h_RelIso_all"); |
| 54 |
|
// old |
| 55 |
|
ax0 = 0.95; |
| 56 |
|
axf = 1.0; |
| 57 |
< |
bx0 = 0.2; |
| 58 |
< |
bxf = 0.8; |
| 57 |
> |
bx0 = 0.0; |
| 58 |
> |
bxf = 0.7; |
| 59 |
|
} |
| 60 |
|
else if ( isoname.Contains("New") ) { |
| 61 |
|
// new |
| 62 |
|
ax0 = 0.0; |
| 63 |
< |
axf = 0.053; |
| 64 |
< |
bx0 = 0.4; |
| 65 |
< |
bxf = 2.0; |
| 63 |
> |
axf = 0.05; // default: 0.053; loose:0.11 |
| 64 |
> |
bx0 = 0.2; |
| 65 |
> |
bxf = 2.; |
| 66 |
|
} |
| 67 |
|
|
| 68 |
|
int niter = 1; |
| 78 |
|
TH1D *h_qcd = (TH1D*)h1->Clone(); |
| 79 |
|
h_all->SetLineWidth(2); |
| 80 |
|
h_all->SetMinimum(0); |
| 81 |
< |
if( jbin == 4) |
| 81 |
> |
if( jbin == 0) |
| 82 |
> |
h_all->SetTitle(isoname+" RelIso distribution (#geq 1 jet)"); |
| 83 |
> |
else if( jbin == 4) |
| 84 |
|
h_all->SetTitle(isoname+" RelIso distribution (#geq 4 jets)"); |
| 85 |
|
else if ( jbin < 4 ) { |
| 86 |
|
TString title = isoname+" RelIso distribution ("; |
| 87 |
|
title += jbin; |
| 88 |
< |
title += " jets)"; |
| 88 |
> |
if (jbin == 1) |
| 89 |
> |
title += " jet)"; |
| 90 |
> |
else |
| 91 |
> |
title += " jets)"; |
| 92 |
|
h_all->SetTitle(title); |
| 93 |
|
} |
| 94 |
|
else{ |
| 96 |
|
} |
| 97 |
|
h_all->Draw(); |
| 98 |
|
|
| 99 |
+ |
Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0); |
| 100 |
+ |
h_all->GetXaxis()->SetRange(2,h_all->GetXaxis()->FindBin(bxf)); |
| 101 |
+ |
Int_t nbMax = h_all->GetMaximumBin(); |
| 102 |
+ |
cout << "Bin number at peak= " << nbMax << endl; |
| 103 |
+ |
h_all->GetXaxis()->SetRange(1,nbxMax); |
| 104 |
+ |
TAxis *axis0 = h_all->GetXaxis(); |
| 105 |
+ |
double bwidth = axis0->GetBinWidth(nbMax); |
| 106 |
+ |
cout << "bwidth = " << bwidth << endl; |
| 107 |
+ |
|
| 108 |
|
cout << "\n>>>>> Iteration step: "<< niter << endl; |
| 109 |
|
h_all->Fit(s0,"0","",bx0,bxf); |
| 110 |
|
|
| 116 |
|
pass[0] = chi2/ndof; |
| 117 |
|
pass[1] = pass[0]; |
| 118 |
|
cout << ">> Norm chi2 = " << pass[0] << endl; |
| 119 |
< |
|
| 107 |
< |
bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2); |
| 108 |
< |
bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2); |
| 109 |
< |
|
| 119 |
> |
|
| 120 |
|
double * par0 = new double [npar]; |
| 121 |
|
double * parErr0 = new double [npar]; |
| 122 |
< |
|
| 122 |
> |
|
| 123 |
> |
for (Int_t i=0;i<npar;i++) { |
| 124 |
> |
printf("%s = %g +- %g\n", |
| 125 |
> |
f0->GetParName(i), |
| 126 |
> |
f0->GetParameter(i), |
| 127 |
> |
f0->GetParError(i) |
| 128 |
> |
); |
| 129 |
> |
par0[i] = f0->GetParameter(i); |
| 130 |
> |
parErr0[i] = f0->GetParError(i); |
| 131 |
> |
} |
| 132 |
> |
|
| 133 |
> |
if(dynamic) { |
| 134 |
> |
bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2); |
| 135 |
> |
bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2); |
| 136 |
> |
} |
| 137 |
|
double delta = pass[1]-pass[2]; |
| 138 |
|
|
| 139 |
< |
while (niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) { |
| 139 |
> |
while (dynamic && niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) { |
| 140 |
|
if (niter > 1) |
| 141 |
|
pass[1]=pass[2]; |
| 142 |
|
if (delta >= 0){ |
| 202 |
|
h_qcd->SetFillColor(38); |
| 203 |
|
h_qcd->Draw("same"); |
| 204 |
|
f1->Draw("same"); |
| 205 |
< |
|
| 205 |
> |
|
| 206 |
|
// Connect fitted and extrapolated region |
| 207 |
|
TF1 *fin = (TF1*)f1->Clone(); |
| 208 |
|
fin->SetRange(axf,bx0); |
