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.5 by jengbou, Wed May 6 05:55:45 2009 UTC vs.
Revision 1.6 by jengbou, Fri Oct 30 17:55:46 2009 UTC

# Line 21 | Line 21 | using namespace std;
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";
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 = 1.75; // 1.82 for loose 2.1 for tight 1.81875/1.6625 for full; 1.9875 for Fast
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"; // Only New for Landau
39 < bool drawtail = true;
39 > bool drawtail = false;
40   bool debug = false;
41 +
42   void fitLandau()
43   {
44    gROOT->SetStyle("CMS");
45 <  //   gStyle->SetOptStat("n");
45 >  //   gStyle->SetOptStat(1);
46    //   gStyle->SetOptTitle(1);
47 <  gStyle->SetOptFit(01111);
47 >  gStyle->SetOptFit(11111);
48    
49    TFile *file = TFile::Open(fileName);
50    TTree *tree = (TTree*)file->Get("myTree");
# Line 55 | Line 57 | void fitLandau()
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;
# Line 69 | Line 72 | void fitLandau()
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.053; // default: 0.053; loose:0.11
82 <    loosecut = 0.11;
83 <    bx0 = 0.2; // 0.2
84 <    bxf = 2.; // 1.2
81 >    axf = 1./19.; // ~0.053
82 >    loosecut = 1./9.; // ~0.11
83 >    bx0 = 0.2;
84 >    bxf = 1.2;
85    }
86    
87    int niter = 1;
# Line 87 | Line 90 | void fitLandau()
90    double ns[2] = {0.,0.};
91    double ns2 = 0.;
92    
93 <  TCanvas *cv = new TCanvas("cv","cv",600,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 95 | Line 98 | void fitLandau()
98    TH1D *h_qcd = (TH1D*)h1->Clone();
99    h_all->SetLineWidth(2);
100    h_all->SetMinimum(0);
98 –  if( jbin == 0)
99 –    h_all->SetTitle(isoname+" RelIso distribution (#geq 1 jet)");
100 –  else if( jbin == 4)
101 –    h_all->SetTitle(isoname+" RelIso distribution (#geq 4 jets)");
102 –  else if ( jbin < 4 ) {
103 –    TString title = isoname+" RelIso distribution (";
104 –    title += jbin;
105 –    if (jbin == 1)
106 –      title += " jet)";
107 –    else
108 –      title += " jets)";
109 –    h_all->SetTitle(title);
110 –  }
111 –  else{
112 –    cout << "Wrong bin number!" << endl;
113 –  }
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)
# Line 157 | Line 145 | void fitLandau()
145    h_all->GetYaxis()->SetTitleSize(0.035);
146    h_all->GetYaxis()->SetTitleOffset(1.6);
147    h_all->GetXaxis()->SetRangeUser(0.,1.6);
160 –  //   h_all->Draw();
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);
# Line 175 | Line 164 | void fitLandau()
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 212 | Line 203 | void fitLandau()
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) {
208          bound[0] = bx0;
209          bound[1] = bxf;
# Line 249 | Line 241 | void fitLandau()
241      }
242    }
243    if(debug){
244 <    cout << ">>> First dynamic fit norm. chi2 = " << pass[0] << endl;
245 <    cout << ">>> Min dynamic norm. chi2       = " << pass[1] << endl;
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] << ", " << bound[1] << ")"<< endl;
247 >    cout << ">>> Best fitting region (bx0, bxf) = (" << bound[0] << ", ";
248 >    cout << bound[1] << ")"<< endl;
249    }
250 <  // Get number of events within fitted (control) region
250 >  
251 >  // Get number of events within fitted region
252    TAxis *axis = h_all->GetXaxis();
253    int bmin = axis->FindBin(bound[0]);
254    int bmax = axis->FindBin(bound[1]);
# Line 265 | Line 259 | void fitLandau()
259    
260    // ////////////////////////
261    // Histos and extrapolation
262 <  
262 >  // ////////////////////////  
263    TF1 *f1 = (TF1*)f0->Clone();
264    f1->SetRange(bound[0],bound[1]);
265    for (int i=0;i<npar;i++) {
# Line 274 | Line 268 | void fitLandau()
268    }
269    f1->SetLineColor(2);
270    
277 –  TPaveStats *p0 = (TPaveStats*)h_qcd->GetListOfFunctions()->FindObject("stats");
278 –  h_qcd->GetListOfFunctions()->Remove(p0);
279 –  h_qcd->SetStats(0);
271    h_qcd->SetLineWidth(2);
272    h_qcd->SetLineColor(40);
273    h_qcd->SetFillStyle(3018);
274    h_qcd->SetFillColor(38);
275 <  h_qcd->Draw("sames hist");
276 <  //   h_qcd->Draw("hist sames");
286 <  f1->Draw("sames");
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,bound[0]);
281    fin->SetLineStyle(2);
282    fin->SetLineColor(8);
283 <  fin->Draw("sames");
284 <
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("sames");
290 >    fine->Draw("same");
291    }
292    
293    TF1 *f2 = (TF1*)f1->Clone();
294    f2->SetRange(ax0,axf);
295    f2->SetLineColor(4);
296 <  f2->Draw("sames");
296 >  f2->Draw("same");
297    int np = 100;
298    double *x=new double[np];
299    double *w=new double[np];
300    f2->CalcGaussLegendreSamplingPoints(np,x,w,1e-15);
301 <
302 <  if (fileName.Contains("0053") || fileName.Contains("TandL")) {
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")) {
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;
# Line 328 | Line 321 | void fitLandau()
321      cout << ">>> Loose Ns by usual integral    = " << ns[1] << endl;
322    else
323      cout << ">>> Tight Ns by usual integral    = " << ns[0] << endl;
324 <
324 >  
325    delete [] x;
326    delete [] w;
327    
# Line 337 | 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 <  if (fileName.Contains("0053") || fileName.Contains("TandL"))
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"))
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]));
# Line 357 | Line 353 | void fitLandau()
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);
# Line 385 | Line 383 | void fitLandau()
383    // Print results
384    cout << "\n";
385    cout << ">>>>> Jet bin " << jbin << ":" << endl;
386 <  cout << ">>>>> Bin number @ peak = " << nbMax << 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 norm. chi2       = " << pass[1] << endl;
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    }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines