ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/dhidas/OSUAnalysis/Tools/scripts/fit_new.C
Revision: 1.1
Committed: Thu Dec 1 16:28:48 2011 UTC (13 years, 5 months ago) by dhidas
Content type: text/plain
Branch point for: dhidas, MAIN
Log Message:
Initial revision

File Contents

# User Rev Content
1 dhidas 1.1 //
2     // 15 Jul 2010
3     //
4     // * Perform fit to RelIso distribution to estimate QCD in various jet bins.
5     // * Different fitting procedures are employed:
6     // (a) free functional fits in all jet bins
7     // (b) free fit in 1,2j, but constrained fits in 3,4mj, using
8     // fitted results from 1,2j.
9     // (c) plot RelIso without fitting.
10     //
11     // * Configurable options:
12     // (a) fit function: pol3, gaus
13     // (b) fit range
14     // (c) whether to zoom in on y-axis or not
15     //
16     // * Each fit produces one such multi-pad plot:
17     // ----------------
18     // | allj | 0j |
19     // ----------------
20     // | 1j | 2j |
21     // ----------------
22     // | 3j | 4mj |
23     // ----------------
24     //
25     // Tips/reminder:
26     // - To use only MC, use run_mc.sh
27     // - To fit data, use run_data.sh
28     // - Every time new data is used, the integrated luminosity needs to be updated (intlumi).
29     //
30     // Things to look out for:
31     // - Because data is plotted first, sometimes the MC histograms are
32     // cut off in the y-scale. To cure this, set explicitly the upper limit on y
33     // using the YMAX_nj (n=all,1,2 etc) global variables.
34     //
35     // If code crashes, turn on the (global) 'debug' flag and re-run, but hopefully you won't need this!
36     //
37     //----------------------------------------------------
38     #include <iomanip>
39     #include <vector>
40     #include "colour_scheme.h"
41    
42     //-------------
43     bool draw_mc = true;
44     bool draw_data = false; //real data
45    
46    
47     string mcFile = "/Users/Frankie/2010/DataRunsOctober/MC301010/all_mc_mixture.root";
48    
49     string dataFile = "/Users/Frankie/2010/DataRunsOctober/Data2010AB/data_data.root";
50    
51    
52    
53     double intlumi = 10.0;//78e-3; //pb-1 <-- measured in data
54     double mcIntlumi = 10.0; //pb-1
55    
56     //-------------
57    
58     bool limit_gaus_mean_12j = false;
59     double gaus_mean_min = 0.3;
60     const double gaus_mean_max = 0.6;
61    
62     const int nj = 6;
63    
64     const bool debug = 0;// true;
65    
66     //--------------------------------
67     // to customize y-axis max range
68     // if set to -1, it means do not set
69     //--------------------------------
70     double YMAX_allj = -1;
71     double YMAX_0j = -1;
72     double YMAX_1j = -1;
73     double YMAX_2j = -1; //20;
74     double YMAX_3j = -1;
75     double YMAX_4mj = -1;
76    
77     //--------------------------------
78    
79     const double sig_from = 0; //new reliso
80     const double sig_upto = 0.1;
81     const double control_max = 1.6;
82     double fit_from;
83     double fit_upto;
84    
85     double this_bw;
86    
87     double nqcd_actual_sig[nj];
88     double n_extrap[nj];
89     double fit_chi2[nj];
90     int fit_ndf[nj];
91    
92     string do_fit;
93    
94     int nFreePara;
95     double fit_par0[nj];
96     double fit_par1[nj];
97     double fit_par2[nj];
98     double fit_par3[nj];
99    
100     // get integral of nEntries from 0-1.6
101     const double control_max = 1.6;
102     double n_cand_data_observed[nj];
103     double n_cand_mc_predicted[nj];
104    
105     int fit_new() {
106     fit();
107     return 0;
108     }
109    
110     void SetDrawMC(int val) {
111     draw_mc = val;
112     }
113     void SetDrawData(int val) {
114     draw_data = val;
115     }
116    
117     void SetDataFile(string p) {
118     dataFile = p;
119     }
120     void SetMCFile(string p) {
121     mcFile = p;
122     }
123    
124     void SetDataIntLumi(double p) {
125     intlumi = p;
126     }
127    
128     int fit(string do_fit_user = "fix", string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.4,
129     bool zoom_in = true) {
130    
131     // Apply TDR style
132     setTDRStyle();
133     // slight adaptation
134     tdrStyle->SetPadRightMargin(0.05); //originally was 0.02, too narrow!
135     gROOT->ForceStyle();
136    
137     double bw = 0.1;
138     int rb = bw / 0.01;
139    
140     string func = func_user;
141    
142     do_fit = do_fit_user;
143    
144     fit_from = fit_from_user; //0.3; // <------
145     fit_upto = fit_upto_user; //1.0; // <------
146    
147    
148     const int rebin[6] = { 10, 10, 10, 10, 10, 10 }; //<----
149    
150    
151     gStyle->SetTitleH(0.1);
152    
153     gStyle->SetStatH(0.22); //0.24);
154     gStyle->SetStatW(0.22); //0.26);
155     gStyle->SetOptStat(1); //on stat
156     gStyle->SetLineScalePS(2); //D=3
157    
158     if (do_fit == "none")
159     gStyle->SetOptStat(0);//off stat
160    
161    
162     TFile *fdata;
163     TFile *fmc;
164    
165     cout << "data file: " << dataFile << endl;
166     cout << "mc file: " << mcFile << endl;
167    
168     if (draw_data)
169     fdata = new TFile(Form("%s", dataFile.c_str()));
170     else
171     fdata = new TFile(Form("%s", mcFile.c_str()));
172     if (draw_mc)
173     fmc = new TFile(Form("%s", mcFile.c_str()));
174    
175     // fixPara: Fit parameters
176     if (do_fit == "fix") {
177     TF1 ff(func.c_str(), func.c_str(), 0, 1);
178     nFreePara = ff.GetNumberFreeParameters();
179     }
180    
181     TCanvas *c1 = new TCanvas("c1", "c1", 980, 1050);//980,700
182    
183     c1->Divide(2, 2);
184    
185     char *njet[6] = { "allj", "0j", "1j", "2j", "3j", "4mj" };
186     char *njlabel[6] = { "#geq0 jets", "0 jet", "1 jet", "2 jets", "3 jets", "#geq4 jets" };
187    
188     if (draw_data && draw_mc) {
189     double sf = intlumi / mcIntlumi;
190     cout << "-----------------------" << endl;
191     cout << "data int lumi: " << intlumi << endl;
192     cout << "mc int lumi: " << mcIntlumi << endl;
193     cout << "scale factor: " << sf << endl;
194     cout << "------------------------" << endl;
195     }
196    
197     if (debug)
198     cout << "About to enter nj loop" << endl;
199    
200     int npad = 1;
201    
202     for (int j = 2; j < 6; j++) {
203    
204     if (debug) {
205     cout << "\n Njet: " << j << " (" << njlabel[j] << ")" << endl;
206     cout << "----------" << endl;
207     }
208     TH1D *all = (TH1D*) fdata->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__data", njet[j]));
209     if (draw_mc) {
210     TH1D *qcd = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__QCD", njet[j]));
211     TH1D *enri1 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__enri1", njet[j]));
212     TH1D *enri2 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__enri2", njet[j]));
213     TH1D *enri3 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__enri3", njet[j]));
214     TH1D *bce1 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__bce1", njet[j]));
215     TH1D *bce2 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__bce2", njet[j]));
216     TH1D *bce3 = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__bce3", njet[j]));
217    
218     TH1D *tt = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__ttbar", njet[j]));
219     TH1D *wj = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__wj", njet[j]));
220     TH1D *zj = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__zj", njet[j]));
221    
222     TH1D *stop = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__singleTop", njet[j]));
223     TH1D *pj = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__pjet", njet[j]));
224    
225     TH1D *mc_all = (TH1D*) fmc->Get(Form("QCD_estimation/QCDest_CombRelIso_%s__data", njet[j]));
226     }
227    
228     // Scaling
229     if (draw_data && draw_mc) {
230     qcd->Scale(sf);
231     enri1->Scale(sf);
232     enri2->Scale(sf);
233     enri3->Scale(sf);
234     bce1->Scale(sf);
235     bce2->Scale(sf);
236     bce3->Scale(sf);
237     tt->Scale(sf);
238     wj->Scale(sf);
239     zj->Scale(sf);
240     stop->Scale(sf);
241     pj->Scale(sf);
242     mc_all->Scale(sf);
243     }
244    
245     const double binw = 0.01; //original binw = 0.01 = 1.1/110 (old) = 10/1000 (new)
246    
247     const int sig_bin_low = (int) (sig_from / binw + 1);
248     const int sig_bin_up = (int) (sig_upto / binw);
249     const int control_bin_max = (int) (control_max / binw);
250    
251     nqcd_actual_sig[j] = -1;
252     if (draw_mc) {
253     nqcd_actual_sig[j] = qcd->Integral(sig_bin_low, sig_bin_up); //actual number of _QCD_ event in signal region
254    
255     }
256     if (1 || debug) {
257     printf("sig_bin_low = %.2f\n", sig_bin_low);
258     printf("sig_bin_up = %.2f\n", sig_bin_up);
259     printf("qcd->Integral( sig_bin_low, sig_bin_up ) = %.2f\n", nqcd_actual_sig[j]);
260     }
261    
262     // get number of candidates
263     cout << "control_bin_max: " << control_bin_max << endl;
264     if (draw_mc)
265     n_cand_mc_predicted[j] = mc_all->Integral(1, control_bin_max);
266     n_cand_data_observed[j] = all->Integral(1, control_bin_max);
267    
268     if (debug)
269     cout << "amtb 1: Rebin, before" << endl;
270    
271     all->Rebin(rebin[j]);
272     if (draw_mc) {
273     qcd->Rebin(rebin[j]);
274     enri1->Rebin(rebin[j]);
275     enri2->Rebin(rebin[j]);
276     enri3->Rebin(rebin[j]);
277     //bce1->Rebin(rebin[j]);
278     bce2->Rebin(rebin[j]);
279     bce3->Rebin(rebin[j]);
280    
281     tt->Rebin(rebin[j]);
282     wj->Rebin(rebin[j]);
283     zj->Rebin(rebin[j]);
284    
285     stop->Rebin(rebin[j]);
286     pj->Rebin(rebin[j]);
287    
288     mc_all->Rebin(rebin[j]);
289     }
290     if (debug)
291     cout << "amtb 1: Rebin, after" << endl;
292    
293     // QCD component
294     THStack *qcdcomp = new THStack("qcdcomp", "QCD components");
295     if (draw_mc) {
296     enri1->SetFillColor(kAzure + 2); //normal
297     enri2->SetFillColor(kAzure + 1);
298     enri3->SetFillColor(kAzure - 9);
299     //bce1->SetFillColor(kSpring+3); //kGreen+2 darkest
300     bce2->SetFillColor(kSpring + 2); //kGreen 8
301     bce3->SetFillColor(kSpring + 1); //kSpring-9
302    
303     tt->SetFillColor(kRed - 7);
304     wj->SetFillColor(kGreen - 9);
305     zj->SetFillColor(kYellow - 7);
306    
307     stop->SetFillColor(kViolet - 4);
308     pj->SetFillColor(kGray + 1);
309    
310     // add from bottom to top
311     //qcdcomp->Add(bce1);
312     qcdcomp->Add(bce2);
313     qcdcomp->Add(bce3);
314     qcdcomp->Add(enri1);
315     qcdcomp->Add(enri2);
316     qcdcomp->Add(enri3);
317    
318     qcdcomp->Add(pj);
319     qcdcomp->Add(zj);
320     qcdcomp->Add(wj);
321    
322     qcdcomp->Add(stop);
323     qcdcomp->Add(tt);
324     }
325     // double this_bw = binw*rebin[j];
326     this_bw = binw * rebin[j];
327    
328     if (debug)
329     cout << "amtb 2" << endl;
330    
331     all->GetXaxis()->SetRangeUser(0, 1.6 - 0.01);
332     if (draw_mc) {
333     qcd->GetXaxis()->SetRangeUser(0, 1.6 - 0.01);
334     qcd->SetLineColor(kAzure + 2);
335     qcd->SetFillColor(kAzure - 9);
336     }
337    
338     all->SetMarkerStyle(20);
339     all->SetMarkerSize(0.6); //if run in batch mode
340    
341     all->SetLineColor(kRed);
342     all->SetTitle(Form("%s", njlabel[j]));
343     all->SetXTitle("RelIso");
344    
345     all->SetName(Form("%s (%.1f-%.1f)", func.c_str(), fit_from, fit_upto));
346    
347     all->SetLineColor(kBlack);
348    
349     c1->cd(npad);
350     npad++;
351    
352     gStyle->SetOptFit(112);
353    
354     // To zoom in on y-axis
355     if (zoom_in) {
356    
357     // find tallest bin from 3rd to 9th (0.2-1.0)
358     double highest = 0;
359     unsigned int scanfrom = all->FindBin(0.2);
360     unsigned int scanto = all->FindBin(1.0) - 1;//do not count bin 1.0
361     for (short i = scanfrom; i < scanto; i++) {
362     double bin_height = all->GetBinContent(i);
363     if (bin_height > highest)
364     highest = bin_height;
365     }
366     all->GetYaxis()->SetRangeUser(0, highest * 2.3);
367    
368     } else {
369    
370     // customize unzoom plot
371     if (njet[j] == "allj" && YMAX_allj > 0)
372     all->GetYaxis()->SetRangeUser(0, YMAX_allj);
373     if (njet[j] == "0j" && YMAX_0j > 0)
374     all->GetYaxis()->SetRangeUser(0, YMAX_0j);
375     if (njet[j] == "1j" && YMAX_1j > 0)
376     all->GetYaxis()->SetRangeUser(0, YMAX_1j);
377     if (njet[j] == "2j" && YMAX_2j > 0)
378     all->GetYaxis()->SetRangeUser(0, YMAX_2j);
379     if (njet[j] == "3j" && YMAX_3j > 0)
380     all->GetYaxis()->SetRangeUser(0, YMAX_3j);
381     if (njet[j] == "4mj" && YMAX_4mj > 0)
382     all->GetYaxis()->SetRangeUser(0, YMAX_4mj);
383     }
384    
385     if (debug)
386     cout << "amtb: gPad->ls()" << endl;
387    
388    
389     if (debug)
390     cout << "amtb: set margin" << endl;
391    
392    
393     gPad->SetBottomMargin(0.14); //TEST
394    
395     if (debug)
396     cout << "amtb: set offset" << endl;
397    
398     all->GetXaxis()->SetTitleOffset(0.8); //on stat
399    
400     if (debug)
401     cout << "amtb: set title size" << endl;
402    
403     all->GetXaxis()->SetTitleSize(0.07);
404    
405     if (debug)
406     cout << "amtb: draw" << endl;
407    
408     all->Draw("ae");
409     if (draw_mc) {
410     qcdcomp->Draw("ahist same"); //add
411     }
412     all->Draw("ae same");
413     all->Draw("axis same");
414    
415     if (debug)
416     cout << "amtb: update gPad" << endl;
417    
418     gPad->Update();
419    
420     if (debug)
421     cout << "amtb 3: get stat box" << endl;
422    
423     TPaveStats *s1 = (TPaveStats*) all->GetListOfFunctions()->FindObject("stats");
424     if (s1 > 0) {
425     double sp[4] = { 0.35, 0.55, 0.7, 0.99 };
426     s1->SetX1NDC(sp[0]);
427     s1->SetY1NDC(sp[1]);
428     s1->SetX2NDC(sp[2]);
429     s1->SetY2NDC(sp[3]);
430     }
431    
432     if (debug)
433     cout << "amtb 4" << endl;
434    
435     if (do_fit == "none") {
436     continue;
437     }
438    
439     if (debug)
440     cout << "perform the fit" << endl;
441     //--------------------
442     // Perform the Fit
443     //--------------------
444    
445    
446     if (do_fit == "free") {
447    
448     //----------
449     // Free Fit
450     //----------
451     if (debug)
452     cout << "free fit" << endl;
453    
454     //1st opt: "V": verbose, "Q": Quiet, "0": do not draw
455     all->Fit(func.c_str(), "Q0", "ah", fit_from, fit_upto);
456    
457    
458    
459     getFitRes(j, gPad, func, all);
460    
461     } else if (do_fit == "fix") {
462    
463     //-------------------
464     // Constrained Fit
465     //-------------------
466     if (debug)
467     cout << "constrained fit" << endl;
468    
469     TF1 *fitf;
470     if (func == "gaus") {
471     fitf = new TF1("gaus", "gaus", 0, 2);
472     } else if (func == "pol3") {
473     // pol3 = A + Bx + Cx^2+ Dx^3 = a(1 + b x + c x2 + d x3)
474     //
475     fitf = new TF1("pol3", "[0] * ( 1 + [1]*x + [2]*x^2 + [3]*x^3 )", 0, 2);
476     } else if (func == "landau") {
477     fitf = new TF1("landau", "landau", 0, 2);
478     }
479    
480     // (a) free-fit in 1,2 jet bins
481     //
482     // For gaussian, require mean to be whithin 0.2-0.6, the flat or peak region.
483     //
484     //if( j<2 ) { //1,2j - Free Fit OR gaus-mean-constrained //this doesn't work now that 0mj and 0j were added
485     if (j < 4) { //1,2j - Free Fit OR gaus-mean-constrained
486    
487     if (func == "gaus" && limit_gaus_mean_12j && j > 1 && j < 4) {//j=2,3 being 1jet or 2jets
488    
489     cout << "constraining gaus mean in12j to " << gaus_mean_min << "-" << gaus_mean_max << endl;
490     fitf->SetParLimits(1, gaus_mean_min, gaus_mean_max); //mean of gaus (1j)
491    
492     }
493    
494     // 1st opt: "V"=verbose; "Q"=quite
495     all->Fit(fitf, "Q0", "ah", fit_from, fit_upto); //Free Fit (Q)
496    
497     } else { //3,4mj - Constrained Fit
498    
499     string fopt;
500     if (func == "gaus") {
501     fitf->FixParameter(1, (fit_par1[0] + fit_par1[1]) / 2); //mean or p1
502     fitf->FixParameter(2, (fit_par2[0] + fit_par2[1]) / 2); //sigma or p2
503     fopt = "BQ0";
504     } else if (func == "pol3") {
505     // fix b,c,d
506     double avg_b = (fit_par1[0] + fit_par1[1]) / 2;
507     double avg_c = (fit_par2[0] + fit_par2[1]) / 2;
508     double avg_d = (fit_par3[0] + fit_par3[1]) / 2;
509    
510     fitf->FixParameter(1, avg_b);
511     fitf->FixParameter(2, avg_c);
512     fitf->FixParameter(3, avg_d);
513     fopt = "BQ0";
514     } else if (func == "landau") {
515     double mpv_min = min(fit_par1[0], fit_par1[1]);//mpv or p1
516     double mpv_max = max(fit_par1[0], fit_par1[1]);
517     cout << "mpv_min " << mpv_min << endl;
518     cout << "mpv_max " << mpv_max << endl;
519     fitf->SetParLimits(1, mpv_min, mpv_max);
520     fopt = "Q0"; //"B" gives wrong answer
521     }
522    
523     // "B"=use fixed parameter value
524     all->Fit(fitf, fopt.c_str(), "ah", fit_from, fit_upto);
525     }
526     delete fitf;
527    
528     if (debug)
529     cout << "getFitRes 1" << endl;
530     getFitRes(j, gPad, func, all);
531    
532     }
533    
534     if (debug)
535     cout << "amtb go to next nj" << endl;
536    
537     }// njet loop
538    
539     if (debug)
540     cout << "amtb: after nj loop" << endl;
541    
542     // Print
543     printf("\nNumber of electron candidates with RelIso 0-1.6.\n");
544     printf(" mc data\n");
545     printf(" allj %10.2f %10.2f\n", n_cand_mc_predicted[0], n_cand_data_observed[0]);
546     printf(" 0j %10.2f %10.2f\n", n_cand_mc_predicted[1], n_cand_data_observed[1]);
547     printf(" 1j %10.2f %10.2f\n", n_cand_mc_predicted[2], n_cand_data_observed[2]);
548     printf(" 2j %10.2f %10.2f\n", n_cand_mc_predicted[3], n_cand_data_observed[3]);
549     printf(" 3j %10.2f %10.2f\n", n_cand_mc_predicted[4], n_cand_data_observed[4]);
550     printf(" >=4j %10.2f %10.2f\n", n_cand_mc_predicted[5], n_cand_data_observed[5]);
551    
552     cout << endl;
553    
554     // Legend
555     //--------
556     if (debug)
557     cout << "amtb: creating legend" << endl;
558    
559     double coord[4] = { 0.7, 0.4, 0.9, 0.9 };
560    
561     TLegend *leg = new TLegend(coord[0], coord[1], coord[2], coord[3]);
562    
563     TF1 *blue = new TF1("blue", "pol0", 0, 1);
564     TF1 *red = new TF1("red", "pol0", 0, 1);
565     blue->SetLineColor(kBlue);
566     red->SetLineColor(kRed);
567    
568     leg->SetFillColor(0);
569     leg->SetBorderSize(0);
570    
571     if (!draw_data)
572     leg->AddEntry(all, "All MC events", "LP"); //Line+Marker, E(error)
573     else
574     leg->AddEntry(all, "Data", "LP");
575     if (draw_mc) {
576     leg->AddEntry(tt, "t#bar{t}", "f");
577     leg->AddEntry(stop, "Single-Top", "f");
578     leg->AddEntry(wj, "W#rightarrowl#nu", "f");
579     leg->AddEntry(zj, "Z#rightarrowl^{+}l^{-}", "f");
580     leg->AddEntry(pj, "#gamma+jets", "f");
581     leg->AddEntry(enri1, "EME pt 20-30", "f");
582     leg->AddEntry(enri2, "EME pt 30-80", "f");
583     leg->AddEntry(enri3, "EME pt 80-170", "f");
584     //leg->AddEntry(bce1, "BCE pt 20-30", "f");
585     leg->AddEntry(bce2, "BCE pt 30-80", "f");
586     leg->AddEntry(bce3, "BCE pt 80-170", "f");
587     }
588     if (do_fit != "none") {
589     leg->AddEntry(red, "Fit", "l");
590     leg->AddEntry(blue, "Extrapolation", "l");
591     }
592    
593     gStyle->SetOptStat(0);
594     gStyle->SetOptFit(0);
595    
596     if (debug)
597     cout << "amtb: draw legend" << endl;
598     c1->cd(1);
599     leg->Draw();
600     add_cms_label();
601     add_njet_label(1);
602     c1->cd(2);
603     leg->Draw();
604     add_cms_label();
605     add_njet_label(2);
606     c1->cd(3);
607     leg->Draw();
608     add_cms_label();
609     add_njet_label(3);
610     c1->cd(4);
611     leg->Draw();
612     add_cms_label();
613     add_njet_label(4);
614    
615     if (do_fit == "none") {
616    
617     if (debug)
618     cout << "amtb: save plot" << endl;
619     if (!zoom_in) {
620     c1->SaveAs("QCD_reliso.gif");
621     c1->SaveAs("QCD_reliso.eps");
622     gROOT->ProcessLine(".!ps2pdf -dEPSCrop QCD_reliso.eps");
623     } else {
624     c1->SaveAs("QCD_reliso_zoom.gif");
625     c1->SaveAs("QCD_reliso_zoom.eps");
626     gROOT->ProcessLine(".!ps2pdf -dEPSCrop QCD_reliso_zoom.eps");
627     }
628     c1->Close();
629    
630     return 0;
631     }
632    
633     if (debug)
634     cout << "amtb: print fit results" << endl;
635    
636     //TDatime now;
637     //now.Print();
638     cout << "\n-----------------------------------------------" << endl;
639     cout << "Fit " << func << " (range: " << fit_from << "-" << fit_upto << ")" << endl;
640     cout << "-----------------------------------------------" << endl;
641     cout << " " << intlumi << "/pb rb True Estimate Diff" << endl;
642     printf("allj jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[0], nqcd_actual_sig[0], n_extrap[0], (n_extrap[0]
643     / nqcd_actual_sig[0] - 1) * 100);
644     printf(" 0 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[1], nqcd_actual_sig[1], n_extrap[1], (n_extrap[1]
645     / nqcd_actual_sig[1] - 1) * 100);
646     printf(" 1 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[2], nqcd_actual_sig[2], n_extrap[2], (n_extrap[2]
647     / nqcd_actual_sig[2] - 1) * 100);
648     printf(" 2 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[3], nqcd_actual_sig[3], n_extrap[3], (n_extrap[3]
649     / nqcd_actual_sig[3] - 1) * 100);
650     printf(" 3 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[4], nqcd_actual_sig[4], n_extrap[4], (n_extrap[4]
651     / nqcd_actual_sig[4] - 1) * 100);
652     printf(" >=4 jet: %2d %10.1f %10.1f %6.1f %%\n", rebin[5], nqcd_actual_sig[5], n_extrap[5], (n_extrap[5]
653     / nqcd_actual_sig[5] - 1) * 100);
654     cout << "-----------------------------------------------\n" << endl;
655    
656     if (debug)
657     cout << "amtb: print constrained-fit results" << endl;
658    
659     if (do_fit == "fix") {
660    
661     cout << " Fitted parameter values" << endl;
662     cout << "-----------------------------------------------" << endl;
663     if (func == "gaus")
664     cout << " Gaus const mean sigma" << endl;
665     else if (func == "landau")
666     cout << " Landau const mpv sigma" << endl;
667     else
668     cout << " p0 p1 p2 p3" << endl;
669     cout << "-----------------------------------------------" << endl;
670     char *jjlabel[5] = { "0", "1", "2", "3", ">=4" };
671    
672     for (short i = 0; i < 5; ++i) {
673     if (nFreePara == 3)
674     printf(" %3s jet: %10.1f %10.2f %10.2f\n", jjlabel[i], fit_par0[i], fit_par1[i], fit_par2[i]);
675     else if (nFreePara == 4)
676     printf(" %3s jet: %10.1f %10.2f %10.2f %10.2f\n", jjlabel[i], fit_par0[i], fit_par1[i], fit_par2[i],
677     fit_par3[i]);
678     }
679     cout << "-----------------------------------------------" << endl;
680    
681     }
682    
683     if (debug)
684     cout << "amtb: save picture" << endl;
685    
686     string out = Form("%s_r%.1fto%.1f", func.c_str(), fit_from, fit_upto);
687     if (do_fit == "fix")
688     out = Form("fixPara_%s_r%.1fto%.1f", func.c_str(), fit_from, fit_upto);
689    
690     if (!zoom_in) {
691     c1->SaveAs(Form("fit_%s.gif", out));
692     c1->SaveAs(Form("fit_%s.pdf", out));
693    
694     } else {
695     c1->SaveAs(Form("fit_zoom_%s.gif", out));
696     c1->SaveAs(Form("fit_zoom_%s.pdf", out));
697     }
698     c1->Close();
699    
700    
701     if (debug)
702     cout << "amtb: write result to text file" << endl;
703    
704     ofstream myfile;
705    
706     string outText;
707     outText = "est_";
708     if (zoom_in)
709     outText += "zoom_";
710     outText = outText + do_fit + "_" + func; //place where "est_free_gaus.txt" is made
711    
712     if (debug)
713     cout << "amtb: " << outText << endl;
714    
715     myfile.open(Form("%s.txt", outText.c_str()), ios::app);
716    
717     myfile.setf(ios::fixed, ios::floatfield);
718     myfile << endl;
719     for (int k = 2; k < 6; ++k) {
720     myfile << n_extrap[k] / nqcd_actual_sig[k] - 1 << endl;
721     }
722    
723     if (debug)
724     cout << "amtb: close text file" << endl;
725     myfile.close();
726    
727     if (debug)
728     cout << "amtb: complete" << endl;
729    
730     return 0;
731     }
732    
733     //-------------------------------------------------------------------------------
734     // just plot without fitting
735     void fit_none(bool zoom_in = true) {
736     fit("none", "", 0, 1, zoom_in);
737     }
738    
739     //-------------------------------------------------------------------------------
740     // free fit
741     void fit(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6, bool zoom_in = true) {
742     fit("free", func_user, fit_from_user, fit_upto_user, zoom_in);
743     }
744    
745     //-------------------------------------------------------------------------------
746     // constrained fit in 3,4j (Free in 1,2j)
747     void fit_fixPara(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6, bool zoom_in = true) {
748     fit("fix", func_user, fit_from_user, fit_upto_user, zoom_in);
749     }
750     //-------------------------------------------------------------------------------
751     // constrained fit in 3,4j (gaus-mean-constrained in 1,2j > 0.3)
752     void fit_fixPara_mean12j_geq03(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6,
753     bool zoom_in = true) {
754     limit_gaus_mean_12j = true;
755     gaus_mean_min = 0.3;
756     fit("fix", func_user, fit_from_user, fit_upto_user, zoom_in);
757     }
758     //-------------------------------------------------------------------------------
759     // constrained fit in 3,4j (gaus-mean-constrained in 1,2j > 0.4)
760     void fit_fixPara_mean12j_geq04(string func_user = "gaus", double fit_from_user = 0.2, double fit_upto_user = 1.6,
761     bool zoom_in = true) {
762     limit_gaus_mean_12j = true;
763     gaus_mean_min = 0.4;
764     fit("fix", func_user, fit_from_user, fit_upto_user, zoom_in);
765     }
766    
767     //-------------------------------------------------------------------------------
768     void getFitRes(int j, TPad *gpad, string func, TH1D *all) {
769    
770     if (debug)
771     cout << "\n start of getFitRes\n" << endl;
772    
773     // Get fitted function
774     //---------------------
775     if (debug)
776     cout << "get fitted function" << endl;
777    
778     TF1 *myf = all->GetFunction(func.c_str());
779     myf->SetLineColor(kRed);
780    
781     if (debug)
782     cout << "clone myf" << endl;
783     TF1 *myf2 = (TF1*) myf->Clone(); //range 0-0.1
784     myf2->SetLineColor(kBlue);
785     myf2->SetRange(0, 0.1);
786    
787     TF1 *myf3 = (TF1*) myf->Clone(); //range 0.1 to 0.2
788     myf3->SetLineColor(kBlue);
789     myf3->SetLineStyle(kDashed);
790     myf3->SetRange(0.1, fit_from);
791    
792     if (debug)
793     cout << "amtb: go to pad" << endl;
794     gpad->cd();
795    
796     if (debug)
797     cout << "amtb: draw functions" << endl;
798     myf->Draw("same");
799     myf2->Draw("same");
800     myf3->Draw("same");
801    
802     // Get fit results
803     //----------------
804     if (debug)
805     cout << "amtb: get fit results" << endl;
806     fit_chi2[j] = myf->GetChisquare();
807     fit_ndf[j] = myf->GetNDF();
808     n_extrap[j] = myf->Integral(sig_from, sig_upto) / this_bw;
809    
810     if (do_fit == "fix")
811     getParametersForConstrainedFit(j, myf);
812    
813     if (debug)
814     cout << "\n end of getFitRes\n" << endl;
815    
816     }
817    
818     //-------------------------------------------------------------------------------
819     void getParametersForConstrainedFit(int j, TF1 *myf) {
820    
821     // fixPara: Get fit results
822     //--------------------------
823     fit_par0[j] = myf->GetParameter(0);//1st para p0
824     fit_par1[j] = myf->GetParameter(1);//2nd para p1
825     if (nFreePara >= 3) {
826     fit_par2[j] = myf->GetParameter(2);//3rd para p2
827     if (nFreePara >= 4) {
828     fit_par3[j] = myf->GetParameter(3);//4th para p3
829     }
830     }
831     }
832     //-------------------------------------------------------------------------------
833    
834     void add_cms_label() {
835    
836     TPaveText * mytext = new TPaveText(0.3, 0.83, 0.6, 0.93, "NDC");
837     mytext->AddText("CMS Preliminary");
838     mytext->AddText(Form("%.0f nb^{-1} at #sqrt{s} = 7 TeV", intlumi * 1000));
839     mytext->SetFillColor(0);
840     mytext->SetBorderSize(0);
841     mytext->SetTextFont(42);
842     mytext->SetTextAlign(13);
843     //mytext->Draw();
844    
845     }
846     //-------------------------------------------------------------------------------
847    
848     void add_njet_label(int nj) {
849    
850     string nj_label;
851    
852     if (nj == 0)
853     nj_label = "N_{jet}=0";
854     else if (nj == 1)
855     nj_label = "N_{jet}=1";
856     else if (nj == 2)
857     nj_label = "N_{jet}=2";
858     else if (nj == 3)
859     nj_label = "N_{jet}=3";
860     else if (nj == 4)
861     nj_label = "N_{jet}#geq4";
862     else if (nj == 5)
863     nj_label = "N_{jet}#geq0";
864    
865     text = new TLatex(3.57, 23.08, Form("%s", nj_label.c_str()));
866     text->SetNDC();
867     text->SetTextAlign(13);
868     text->SetX(0.4);
869     text->SetY(0.75);
870     text->SetTextFont(42);
871     text->SetTextSize(0.07);
872     text->Draw();
873     }
874    
875     //-------------------------------------------------------------------------------