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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines