ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/scripts/fitLandau.C
(Generate patch)

Comparing UserCode/Jeng/scripts/fitLandau.C (file contents):
Revision 1.3 by jengbou, Mon Apr 27 22:37:00 2009 UTC vs.
Revision 1.6 by jengbou, Fri Oct 30 17:55:46 2009 UTC

# Line 17 | Line 17 | using namespace std;
17   #include "Math/GSLIntegrator.h"
18   #include "Math/WrappedTF1.h"
19  
20 < TString fileName = "skimmed/QCDbin4_005.root";
20 > // Full Sim
21 > // TString fileName = "skimmed226/QCDbin4_NewRelIso_rebin1_TandL.root";
22 > // TString fileName = "skimmed226/QCDbin3_NewRelIso_0053.root";
23 > // TString fileName = "skimmed226/QCDbin2_NewRelIso_wErrors_TandL.root";
24 > TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL_0.root";
25 > // Fast Sim
26 > // TString fileName = "skimmed226/QCDbin3_NewRelIso_FastSim_0053.root";
27 > // TString fileName = "skimmed226/QCDbin4_NewRelIso_FastSim_0053.root";
28 > // TString fileName = "skimmed226/QCDbin4_NewRelIso_wErrors_TandL_FastSim_1.root";
29   // Dynamic(true) or static(false) range
30 +
31   bool dynamic = true;
32   // Dynamic fit region weight
33 < double wxmin = 2.; // 1.82 for loose 2.1 for tight
34 < double wxmax = 0.8; // 0.6 for loose 0.8 for tight
33 > double wxmin = 1.825; // 1.82 for loose 2.1 for tight 1.81875/  1.825 for full; 2.0 for Fast
34 > double wxmax = wxmin/2.; // 0.6 for loose 0.8 for tight
35   // Select fitting function:
36   TString s0 = "landau";
37   // Select New or Old reliso:
38 < TString isoname = "New";
38 > TString isoname = "New"; // Only New for Landau
39 > bool drawtail = false;
40 > bool debug = false;
41  
42   void fitLandau()
43   {
44 <  gStyle->SetOptFit(01111);
44 >  gROOT->SetStyle("CMS");
45 >  //   gStyle->SetOptStat(1);
46 >  //   gStyle->SetOptTitle(1);
47 >  gStyle->SetOptFit(11111);
48    
49    TFile *file = TFile::Open(fileName);
50    TTree *tree = (TTree*)file->Get("myTree");
51    TH1D *h0 = (TH1D*)file->Get("histos/h_RelIso_all");
52    TH1D *h1 = (TH1D*)file->Get("histos/h_RelIso_qcd");
53    Int_t jbin;
54 <  Double_t na;
54 >  Double_t na,na1;
55    tree->SetBranchAddress("jbin",&jbin);
56    tree->SetBranchAddress("na",&na);
57 <  //   Int_t nentries = (Int_t)tree->GetEntries();
57 >  if (fileName.Contains("TandL"))
58 >    tree->SetBranchAddress("na1",&na1);  
59    tree->GetEntry(0);
60 +  
61    // define fit region:
62 <  double ax0 = 0.;
63 <  double axf = 0.;
64 <  double bx0= 0.;
65 <  double bxf = 0.;
62 >  double ax0 = 0.0;
63 >  double axf = 0.0;
64 >  double bx0 = 0.0;
65 >  double bxf = 0.0;
66 >  double loosecut = 0.;
67    
68    cout<<"Use "<<isoname<<" RelIso"<<endl;
69    
70    if ( isoname.Contains("Old") ) {
71      // old
72 +    loosecut = 0.9;
73      ax0 = 0.95;
74      axf = 1.0;
75 <    bx0 = 0.0;
76 <    bxf = 0.7;
75 >    bx0 = 0.3; // Must be smaller than x at peak
76 >    bxf = 0.9; // Should not overlap signal region
77    }
78    else if ( isoname.Contains("New") ) {
79      // new
80      ax0 = 0.0;
81 <    axf = 0.05; // default: 0.053; loose:0.11
81 >    axf = 1./19.; // ~0.053
82 >    loosecut = 1./9.; // ~0.11
83      bx0 = 0.2;
84 <    bxf = 2.;
84 >    bxf = 1.2;
85    }
86    
87    int niter = 1;
88    double pass[3] = {0.,0.,0.};
89 <  double ns, ns2;
90 <  ns = ns2 = 0.;
89 >  double bound[2] = {0.,0.};
90 >  double ns[2] = {0.,0.};
91 >  double ns2 = 0.;
92    
93 <  TCanvas *cv = new TCanvas("cv","cv",800,600);
93 >  TCanvas *cv = new TCanvas("cv","cv",800,800);
94    
95    // Fit QCD in background region
96    cout<<"Fitting with "<<s0<<" function!"<<endl;
# Line 78 | Line 98 | void fitLandau()
98    TH1D *h_qcd = (TH1D*)h1->Clone();
99    h_all->SetLineWidth(2);
100    h_all->SetMinimum(0);
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    if (jbin == 1)
89      title += " jet)";
90    else
91      title += " jets)";
92    h_all->SetTitle(title);
93  }
94  else{
95    cout << "Wrong bin number!" << endl;
96  }
97  h_all->Draw();
101    
102 +  // Get num of bin at peak and bin width
103    Int_t nbxMax = h_all->GetXaxis()->FindBin(2.0);
104    h_all->GetXaxis()->SetRange(2,h_all->GetXaxis()->FindBin(bxf));
105    Int_t nbMax = h_all->GetMaximumBin();
106 <  cout << "Bin number at peak= " << nbMax << endl;
106 >  cout << ">>>>> Bin number at peak= " << nbMax << endl;
107    h_all->GetXaxis()->SetRange(1,nbxMax);
108    TAxis *axis0 = h_all->GetXaxis();
109    double bwidth = axis0->GetBinWidth(nbMax);
110 <  cout << "bwidth = " << bwidth << endl;
110 >  cout << ">>>>> bwidth = " << bwidth << endl;
111 >  
112 >  TString title = "RelIso'";
113 >  if( jbin == 4)
114 >    title += " (#geq 4 jets)";
115 >  else if ( jbin == 1 ) {
116 >    title += " (";
117 >    title += jbin;
118 >    title += " jet)";
119 >  }
120 >  else if ( jbin < 4 ) {
121 >    title += " (";
122 >    title += jbin;
123 >    title += " jets)";
124 >    h_RelIso_all->SetTitle(title);
125 >  }
126 >  else if ( jbin == 5 ) {
127 >    h_RelIso_all->SetTitle( isoname+" RelIso distribution (#geq 1 jets (up to 3 jets))");
128 >  }
129 >  else {
130 >    cout<<"Don't mess around!!! jbin should be less than 6!"<<endl;
131 >    return;
132 >  }
133 >  //   gStyle->SetErrorX(0);
134 >  h_all->GetXaxis()->SetTitle(title);
135 >  h_all->GetXaxis()->SetLabelFont(42);
136 >  h_all->GetXaxis()->SetLabelSize(0.03);
137 >  h_all->GetXaxis()->SetTitleFont(42);
138 >  h_all->GetXaxis()->SetTitleSize(0.035);
139 >  h_all->GetXaxis()->SetTitleOffset(1.2);
140 >  
141 >  h_all->GetYaxis()->SetTitle("Events ( L_{Int}= 20 pb^{-1} )");
142 >  h_all->GetYaxis()->SetLabelFont(42);
143 >  h_all->GetYaxis()->SetLabelSize(0.03);
144 >  h_all->GetYaxis()->SetTitleFont(42);
145 >  h_all->GetYaxis()->SetTitleSize(0.035);
146 >  h_all->GetYaxis()->SetTitleOffset(1.6);
147 >  h_all->GetXaxis()->SetRangeUser(0.,1.6);
148 >  h_all->Draw("e");
149    
150    cout << "\n>>>>> Iteration step: "<< niter << endl;
151 <  h_all->Fit(s0,"0","",bx0,bxf);
151 >  cout << "\n>>>>> Start fitting control region: " << bx0 << " to " << bxf << endl;
152 >  
153 >  h_all->Fit(s0,"","",bx0,bxf);
154    
155    // Very first fit function
156    TF1 *f0 = (TF1*)h_all->GetFunction(s0);
157    Int_t npar = f0->GetNpar();
158    Double_t chi2 = f0->GetChisquare();
159    Int_t ndof = f0->GetNDF();
160 <  pass[0] = chi2/ndof;
161 <  pass[1] = pass[0];
162 <  cout << ">> Norm chi2 = " << pass[0] << endl;
163 <    
160 >  pass[0] = chi2/ndof; // Very first fit norm chi2. Fixed range.
161 >  cout << ">>>> Norm chi2 = " << pass[0] << endl;
162 >  if (!dynamic)
163 >    pass[1] = pass[0];
164 >  else
165 >    pass[1] = 999.;
166 >  
167 >  cout << " >>> Peak at " << f0->GetMaximumX(bx0,bxf) << endl;
168 >  
169    double * par0 = new double [npar];
170    double * parErr0 = new double [npar];
171    
# Line 129 | Line 178 | void fitLandau()
178      par0[i] = f0->GetParameter(i);
179      parErr0[i] = f0->GetParError(i);
180    }
181 <  
181 >  bound[0] = bx0;
182 >  bound[1] = bxf;  
183    if(dynamic) {
184      bx0 = f0->GetParameter(1) - wxmin*f0->GetParameter(2);
185      bxf = f0->GetParameter(1) + wxmax*f0->GetParameter(2);
186    }
187 +  
188    double delta = pass[1]-pass[2];
189    
190 <  while (dynamic && niter <= 20 && abs(pass[1]-pass[2]) > 0.00001) {
190 >  while (dynamic && niter <= 20 && abs(delta) > 0.00001) {
191      if (niter > 1)
192 <      pass[1]=pass[2];
193 <    if (delta >= 0){
192 >      pass[1] = pass[2];
193 >    if (delta != 0){
194        niter++;
195        cout << "\n>>>>> Iteration step: "<< niter << endl;
196 <      cout << ">> Pass[1] = " << pass[1] << endl;
196 >      cout << ">> Previous step norm. chi2 = " << pass[1] << endl;
197 >      cout << ">>>> Starting fitting new range (bx0, bxf) = (" << bx0;
198 >      cout << ", " << bxf << ")" << endl;
199        
200        h_all->Fit(s0,"0","",bx0,bxf);
201 <      TF1 *fi = (TF1*)h_all->GetFunction(s0);    
201 >      TF1 *fi = (TF1*)h_all->GetFunction(s0);
202        pass[2] = fi->GetChisquare()/fi->GetNDF();
203 <      cout << ">> Norm chi2 = pass[2] = " << pass[2] << endl;
203 >      cout << ">> Current step norm. chi2  = " << pass[2] << endl;
204        delta = pass[1]-pass[2];
205        cout << ">> Delta = " << delta << endl;
206        
207 <      if (delta >= 0) {
207 >      if (delta > 0) {
208 >        bound[0] = bx0;
209 >        bound[1] = bxf;
210          printf("Function has %i parameters. Chisquare = %g\n",
211                 npar,
212                 fi->GetChisquare());
# Line 167 | Line 222 | void fitLandau()
222          }
223          bx0 = fi->GetParameter(1) - wxmin*fi->GetParameter(2);
224          bxf = fi->GetParameter(1) + wxmax*fi->GetParameter(2);
225 <        cout << ">> bx0, bxf = " << bx0 << ", " << bxf << endl;
225 >        cout << ">> Range for next iteration (bx0, bxf) = (" << bx0;
226 >        cout << ", " << bxf << ")" << endl;
227 >      }
228 >      else if (delta < 0) {
229 >        cout << " Use previous fit parameters for extrapolation fit" << endl;
230 >        delta = 0;
231        }
232        else {
233 <        cout << ">> Use previous fit!" << endl;
233 >        //      bound[0] = (bound[0]<bx0)?bound[0]:bx0;
234 >        //      bound[1] = (bound[1]>bxf)?bound[1]:bxf;
235 >        cout << ">> Fit converges." << endl;
236 >        cout << ">> Best fitting region = (" << bound[0];
237 >        cout << ", " << bound[1] << ")" << endl;
238        }
239 +      if (niter == 2)
240 +        pass[0]=pass[2]; // First dynamic fit norm. chi2
241      }
242    }
243 <  
244 <  cout << ">>> Final fitting resion (bx0, bxf) = " << bx0 << ", " << bxf << endl;
243 >  if(debug){
244 >    cout << ">>> First dynamic fit norm. chi2   = " << pass[0] << endl;
245 >    cout << ">>> Min dynamic norm. chi2         = " << pass[1] << endl;
246 >    cout << ">>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
247 >    cout << ">>> Best fitting region (bx0, bxf) = (" << bound[0] << ", ";
248 >    cout << bound[1] << ")"<< endl;
249 >  }
250    
251    // Get number of events within fitted region
252    TAxis *axis = h_all->GetXaxis();
253 <  int bmin = axis->FindBin(bx0);
254 <  int bmax = axis->FindBin(bxf);
253 >  int bmin = axis->FindBin(bound[0]);
254 >  int bmax = axis->FindBin(bound[1]);
255    double nfac1 = h_all->Integral(bmin,bmax);
256 <  nfac1 -= (h_all->GetBinContent(bmin))*(bx0-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
257 <  nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bxf)/axis->GetBinWidth(bmax);
256 >  nfac1 -= (h_all->GetBinContent(bmin))*(bound[0]-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
257 >  nfac1 -= (h_all->GetBinContent(bmax))*(axis->GetBinUpEdge(bmax)-bound[1])/axis->GetBinWidth(bmax);
258    cout << ">>> Final Nfac = " << nfac1 << endl;
259    
260    // ////////////////////////
261    // Histos and extrapolation
262 <  
262 >  // ////////////////////////  
263    TF1 *f1 = (TF1*)f0->Clone();
264 <  f1->SetRange(bx0,bxf);
264 >  f1->SetRange(bound[0],bound[1]);
265    for (int i=0;i<npar;i++) {
266      cout<<" >> Par["<<i<<"]="<<par0[i]<<"; Error="<<parErr0[i]<<endl;
267      f1->SetParameter(i,par0[i]);
268    }
198  
269    f1->SetLineColor(2);
270 +  
271 +  h_qcd->SetLineWidth(2);
272    h_qcd->SetLineColor(40);
273    h_qcd->SetFillStyle(3018);
274    h_qcd->SetFillColor(38);
275 <  h_qcd->Draw("same");
275 >  h_qcd->Draw("same hist");
276    f1->Draw("same");
277    
278    // Connect fitted and extrapolated region
279    TF1 *fin = (TF1*)f1->Clone();
280 <  fin->SetRange(axf,bx0);
280 >  fin->SetRange(axf,bound[0]);
281    fin->SetLineStyle(2);
282    fin->SetLineColor(8);
283    fin->Draw("same");
284    
285 +  if (drawtail){
286 +    TF1 *fine = (TF1*)f1->Clone();
287 +    fine->SetRange(bound[1],2.0);
288 +    fine->SetLineStyle(2);
289 +    fine->SetLineColor(8);
290 +    fine->Draw("same");
291 +  }
292 +  
293    TF1 *f2 = (TF1*)f1->Clone();
294    f2->SetRange(ax0,axf);
295    f2->SetLineColor(4);
# Line 218 | Line 298 | void fitLandau()
298    double *x=new double[np];
299    double *w=new double[np];
300    f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
301 <  ns = f2->IntegralFast(np,x,w,ax0,axf);
302 <  ns/=(f2->Integral(bx0,bxf));
303 <  ns*=nfac1;
304 <  cout<<">>> Ns by usual integral = "<<ns<<endl;
301 >  
302 >  if (fileName.Contains("0053") ||
303 >      fileName.Contains("TandL") ||
304 >      fileName.Contains("Tight")) {
305 >    ns[0] = f2->IntegralFast(np,x,w,ax0,axf);
306 >    ns[0]/=(f2->Integral(bound[0],bound[1]));
307 >    ns[0]*=nfac1;
308 >  }
309 >  else if (fileName.Contains("011") ||
310 >           fileName.Contains("Loose")) {
311 >    ns[0] = f2->IntegralFast(np,x,w,ax0,loosecut);
312 >    ns[0]/=(f2->Integral(bound[0],bound[1]));
313 >    ns[0]*=nfac1;
314 >  }
315 >  
316 >  ns[1] = f2->IntegralFast(np,x,w,ax0,loosecut);
317 >  ns[1]/=(f2->Integral(bound[0],bound[1]));
318 >  ns[1]*=nfac1;
319 >  
320 >  if (fileName.Contains("TandL") || fileName.Contains("011"))
321 >    cout << ">>> Loose Ns by usual integral    = " << ns[1] << endl;
322 >  else
323 >    cout << ">>> Tight Ns by usual integral    = " << ns[0] << endl;
324 >  
325    delete [] x;
326    delete [] w;
327    
# Line 230 | Line 330 | void fitLandau()
330    ROOT::Math::GSLIntegrator ig(1.E-8,1.E-8,1000);
331    ROOT::Math::WrappedTF1 wf(*g);
332    ig.SetFunction(wf);
333 <  ns2 = ig.Integral(ax0,axf);
334 <  ns2/=(g->Integral(bx0,bxf));
333 >  if (fileName.Contains("0053") ||
334 >      fileName.Contains("TandL") ||
335 >      fileName.Contains("Tight"))
336 >    ns2 = ig.Integral(ax0,axf);
337 >  else if (fileName.Contains("011") ||
338 >           fileName.Contains("Loose"))
339 >    ns2 = ig.Integral(ax0,loosecut);
340 >  
341 >  ns2/=(g->Integral(bound[0],bound[1]));
342    ns2*=nfac1;
343 <  cout<<">>> Ns by MathMore integral = "<<ns2<<endl;
344 <    
345 <  cv->Update();
346 <  
347 <  TPaveStats *p1 = (TPaveStats*)h_all->GetListOfFunctions()->FindObject("stats");
348 <  p1->SetTextColor(kBlue);
349 <  p1->SetX1NDC(0.6);
343 >  cout<<">>> Ns by MathMore integral       = "<<ns2<<endl;
344 >  
345 >  if(!dynamic){
346 >    cv->Update();
347 >    TPaveStats *p1 = (TPaveStats*)cv->GetPrimitive("stats");
348 >  }
349 >  else {
350 >    h_all->Fit(s0,"0","",bound[0],bound[1]);
351 >    cv->Update();
352 >    TPaveStats *p1 = (TPaveStats*)cv->GetPrimitive("stats");
353 >  }
354 >  p1->SetName("Landau fit");
355 >  TList *list = p1->GetListOfLines();
356 >  //   TText *tconst = p1->GetLineWith("Constant");
357 >  //   list->Remove(tconst);
358 >  TLatex *myt = new TLatex(0,0,"Landau fit");
359 >  list->AddFirst(myt);
360 >  h_all->SetStats(0);
361 >  cv->Modified();
362 >  
363 >  p1->SetTextFont(42);
364 >  p1->SetFillColor(0);
365 >  p1->SetFillStyle(0);
366 >  p1->SetBorderSize(0);
367 >  p1->SetX1NDC(0.65);
368    p1->SetX2NDC(0.88);
369 <  p1->SetY1NDC(0.62);
370 <  p1->SetY2NDC(0.88);
371 <  p1->Draw();
369 >  p1->SetY1NDC(0.65);
370 >  p1->SetY2NDC(0.85);
371 >  gPad->Update();
372    
373    // Label histo
374 <  TLegend * leg1 = new TLegend(0.3,0.15,0.68,0.45);
374 >  TLegend * leg1 = new TLegend(0.25,0.65,0.5,0.85);
375 >  leg1->SetFillColor(0);
376 >  leg1->SetFillStyle(0);
377 >  leg1->AddEntry(f1,"Fit of all events in control region","l");
378 >  leg1->AddEntry(f2,"Extrapolation to signal region","l");
379    leg1->AddEntry(h_all,"All events (S+B)");
380    leg1->AddEntry(h_qcd,"QCD events");
252  leg1->AddEntry(f1,"Fit of all events in control region");
253  leg1->AddEntry(f2,"Extrapolation to signal region");
381    leg1->Draw();
382 +  
383    // Print results
384    cout << "\n";
385    cout << ">>>>> Jet bin " << jbin << ":" << endl;
386 <  cout << ">>>>> Observed Ns = " << ns << endl;
386 >  cout << ">>>>> Bin number at peak = " << nbMax << endl;
387 >  cout << ">>>>> bwidth = " << bwidth << endl;
388 >  if (dynamic) {
389 >    cout << ">>>>> First dynamic fit norm. chi2   = " << pass[0] << endl;
390 >    cout << ">>>>> Min dynamic fit norm. chi2     = " << pass[1] << endl;
391 >    cout << ">>>>> Last fitting region (bx0, bxf) = (" << bx0 << ", " << bxf << ")"<< endl;
392 >    cout << ">>>>> Best fitting region (bx0, bxf) = (" << bound[0] << ", " << bound[1] << ")"<< endl;
393 >  }
394 >  cout << ">>>>> Observed Ns = " << ns[0] << endl;
395    cout << ">>>>> Expected Na = " << na << endl;
396 <  cout << "Deviation = " << (ns-na)/na*100.0 << " %" <<endl;
397 <
396 >  cout << "Deviation = " << (ns[0]-na)/na*100.0 << " %" <<endl;
397 >  
398 >  if (fileName.Contains("TandL")) {
399 >    cout << "====================================" << endl;
400 >    cout << ">>>>> Observed Loose Ns = " << ns[1] << endl;
401 >    cout << ">>>>> Expected Loose Na = " << na1 << endl;
402 >    cout << "Deviation = " << (ns[1]-na1)/na1*100.0 << " %" <<endl;
403 >  }
404   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines