ra4_macros  bede988c286599a3a84b77a4d788ac0a971e89f9
csv_fit.cxx
Go to the documentation of this file.
1 #include "tdrstyle.C"
2 #include "TH1.h"
3 #include "TFile.h"
4 #include "TString.h"
5 #include "TObjArray.h"
6 #include "TLegend.h"
7 #include "TFractionFitter.h"
8 #include "TCanvas.h"
9 #include "TLatex.h"
10 #include "Fit/Fitter.h"
11 #include "TROOT.h"
12 #include <iostream>
13 
14 int main(int argc, char *argv[])
15 {
16  bool fitCharmWithLight=false;
17  bool excludeHighCSV=false;
18 
19  TString luminosity = "2.7";
20  int maxbin=22;
21 
22  // by default, want to look at low njet region
23  TString fitType="low_njet";
24  // otherwise parse second parameter
25  if(argc>1) {
26  fitType=Form("%s", argv[1]);
27  // options: low_njet, low_njet_low_mj, low_njet_high_mj, high_njet, exclude_high_csv
28  if(fitType!="low_njet" &&
29  fitType!="high_njet" &&
30  fitType!="low_njet_low_mj" &&
31  fitType!="low_njet_high_mj" &&
32  fitType!="exclude_high_csv") {
33  std::cout << "Invalid fit type: " << fitType << std::endl;
34  return 1;
35  }
36  if(fitType=="exclude_high_csv") {
37  maxbin=16;
38  // still want to fit low_njet region
39  fitType="low_njet";
40  }
41  }
42  setTDRStyle();
43  gROOT->ForceStyle();
44 
45  gStyle->SetPadTopMargin(0.08);
46  gStyle->SetPadLeftMargin(0.18);
47  gStyle->SetTitleOffset(1.43,"y");
48 
49  TFile *f = TFile::Open("csv_newmethod.root");
50 
51  std::cout << "Getting weighted histograms" << std::endl;
52  TH1F *qcd_b = static_cast<TH1F*>(f->Get(Form("QCD_b_%s", fitType.Data())));
53  TH1F *qcd_c = static_cast<TH1F*>(f->Get(Form("QCD_c_%s", fitType.Data())));
54  TH1F *qcd_l = static_cast<TH1F*>(f->Get(Form("QCD_l_%s", fitType.Data())));
55  TH1F *other = static_cast<TH1F*>(f->Get(Form("other_%s", fitType.Data())));
56  TH1F *qcd_cl = static_cast<TH1F*>(qcd_c->Clone());
57  qcd_cl->Add(qcd_l);
58 
59  std::cout << "Getting raw histograms" << std::endl;
60  TH1F *data = static_cast<TH1F*>(f->Get(Form("data_%s", fitType.Data())));
61  TH1F *qcd_b_raw = static_cast<TH1F*>(f->Get(Form("QCD_b_%s_raw", fitType.Data())));
62  TH1F *qcd_c_raw = static_cast<TH1F*>(f->Get(Form("QCD_c_%s_raw", fitType.Data())));
63  TH1F *qcd_l_raw = static_cast<TH1F*>(f->Get(Form("QCD_l_%s_raw", fitType.Data())));
64  TH1F *other_raw = static_cast<TH1F*>(f->Get(Form("other_%s_raw", fitType.Data())));
65  TH1F *qcd_cl_raw = static_cast<TH1F*>(qcd_c_raw->Clone());
66  qcd_cl_raw->Add(qcd_l_raw);
67 
68  std::cout << "Getting weights" << std::endl;
69  TH1F *qcd_b_weights = static_cast<TH1F*>(f->Get(Form("QCD_b_%s_weights", fitType.Data())));
70  TH1F *qcd_c_weights = static_cast<TH1F*>(f->Get(Form("QCD_c_%s_weights", fitType.Data())));
71  TH1F *qcd_l_weights = static_cast<TH1F*>(f->Get(Form("QCD_l_%s_weights", fitType.Data())));
72  TH1F *other_weights = static_cast<TH1F*>(f->Get(Form("other_%s_weights", fitType.Data())));
73  TH1F *qcd_cl_weights = static_cast<TH1F*>(qcd_c_weights->Clone());
74  qcd_cl_weights->Add(qcd_l_weights);
75 
76 
77  float qcd_yield = qcd_b->Integral()+qcd_c->Integral()+qcd_l->Integral();
78  float scaling = (data->Integral()-other->Integral())/qcd_yield;
79  std::cout << "Scaling qcd by a factor: " << scaling << std::endl;
80  qcd_b->Scale(scaling);
81  qcd_c->Scale(scaling);
82  qcd_l->Scale(scaling);
83  qcd_cl->Scale(scaling);
84 
85  qcd_b_weights->Scale(scaling);
86  qcd_c_weights->Scale(scaling);
87  qcd_l_weights->Scale(scaling);
88  qcd_cl_weights->Scale(scaling);
89 
90  float total_yield = qcd_b->Integral()+qcd_c->Integral()+qcd_l->Integral()+other->Integral();
91  float qcd_b_fracbefore = qcd_b->Integral()/total_yield;
92  float qcd_c_fracbefore = qcd_c->Integral()/total_yield;
93  float qcd_l_fracbefore = qcd_l->Integral()/total_yield;
94  float qcd_cl_fracbefore = (qcd_c->Integral()+qcd_l->Integral())/total_yield;
95  float other_fracbefore = other->Integral()/total_yield;
96 
97  std::cout << "fractions before fit:\n"
98  << "b: " << qcd_b_fracbefore << std::endl
99  << "c: " << qcd_c_fracbefore << std::endl
100  << "l: " << qcd_l_fracbefore << std::endl
101  << "cl: " << qcd_cl_fracbefore << std::endl
102  << "other: " << other_fracbefore << std::endl;
103 
104  TObjArray *mc = new TObjArray(4);
105  mc->Add(qcd_b_raw);
106  mc->Add(other_raw);
107  if(fitCharmWithLight) {
108  mc->Add(qcd_cl_raw);
109  }
110  else {
111  mc->Add(qcd_c_raw);
112  mc->Add(qcd_l_raw);
113  }
114 
115  TFractionFitter *fit = new TFractionFitter(data, mc,"V");
116  fit->Constrain(0, 1e-4, 1.0);
117  fit->Constrain(1, 1e-4, 1.0);
118  fit->Constrain(2, 1e-4, 1.0);
119  if(!fitCharmWithLight) fit->Constrain(3, 1e-4, 1.0);
120 
121  fit->SetWeight(0, qcd_b_weights);
122  fit->SetWeight(1, other_weights);
123  if(fitCharmWithLight) {
124  fit->SetWeight(2, qcd_cl_weights);
125  }
126  else {
127  fit->SetWeight(2, qcd_c_weights);
128  fit->SetWeight(3, qcd_l_weights);
129  }
130  // allow variation of range
131  fit->SetRangeX(1, maxbin);
132 
133  ROOT::Fit::Fitter* fitter = fit->GetFitter();
134  fitter->Config().ParSettings(1).SetValue(other->Integral()/total_yield);
135  fitter->Config().ParSettings(1).Fix();
136  if(!fitCharmWithLight) {
137  fitter->Config().ParSettings(3).SetValue(qcd_l->Integral()/total_yield);
138  fitter->Config().ParSettings(3).Fix();
139  }
140 
141  Int_t status = fit->Fit();
142  if(status==0) {
143  // get degrees of freedom (the fitter does not calculate the number of degrees of freedom correctly
144  // when there are fixed parameters
145  const std::vector<ROOT::Fit::ParameterSettings>& settings = fitter->Config().ParamsSettings();
146  int ndof = maxbin;
147  for (auto isetting : settings) {
148  if(!isetting.IsFixed()) ndof--;
149  }
150  std::cout << "chi^2/ndof: " << fit->GetChisquare() << "/" << ndof << ", prob = " << fit->GetProb() << std::endl;
151  data->SetMinimum(0);
152  data->GetXaxis()->SetNdivisions(505);
153  data->SetMarkerSize(1);
154  data->SetLineColor(kBlack);
155  data->SetMarkerStyle(kFullCircle);
156  data->GetYaxis()->SetLabelSize(0.04);
157  data->SetTitle(";CSV;Events / 0.005");
158  TCanvas *c = new TCanvas;
159  data->Draw("e");
160 
161  TH1F* result_b = static_cast<TH1F*>(fit->GetMCPrediction(0));
162  TH1F* result_other = static_cast<TH1F*>(fit->GetMCPrediction(1));
163  TH1F* result_c = 0;
164  TH1F* result_l = 0;
165  TH1F* result_cl = 0;
166 
167  if(fitCharmWithLight) {
168  result_cl = static_cast<TH1F*>(fit->GetMCPrediction(2));
169  }
170  else {
171  result_c = static_cast<TH1F*>(fit->GetMCPrediction(2));
172  result_l = static_cast<TH1F*>(fit->GetMCPrediction(3));
173  }
174 
175 
176  result_b->SetLineColor(kRed);
177  result_b->SetLineWidth(3);
178  result_b->SetLineStyle(kDashed);
179  result_b->SetMarkerSize(0);
180  result_b->Multiply(qcd_b_weights);
181  result_b->Draw("same");
182  result_b->Draw("same,hist");
183  if(fitCharmWithLight) {
184  result_cl->SetLineColor(kBlack);
185  result_cl->SetLineWidth(3);
186  result_cl->SetLineStyle(kDashed);
187  result_cl->SetMarkerSize(0);
188  result_cl->Draw("same");
189  result_cl->Draw("same,hist");
190  }
191  else {
192  result_c->SetLineColor(kBlack);
193  result_c->SetLineWidth(3);
194  result_c->SetLineStyle(kDashed);
195  result_c->SetMarkerSize(0);
196  result_c->Multiply(qcd_c_weights);
197  result_c->Draw("same");
198  result_c->Draw("same,hist");
199  result_l->SetLineColor(kMagenta);
200  result_l->SetLineWidth(3);
201  result_l->SetLineStyle(kDashed);
202  result_l->SetMarkerSize(0);
203  result_l->Multiply(qcd_l_weights);
204  result_l->Draw("same");
205  result_l->Draw("same,hist");
206  }
207  result_other->SetLineColor(kGreen);
208  result_other->SetLineWidth(3);
209  result_other->SetLineStyle(kDashed);
210  result_other->Multiply(other_weights);
211  result_other->SetMarkerSize(0);
212  result_other->Draw("same");
213  result_other->Draw("same,hist");
214  TH1F* sum = static_cast<TH1F*>( fit->GetPlot());
215 
216  //Original calculation: not actually correct; see documentation for TFractionFitter::GetMCPrediction()
217  /*TH1F *sum = static_cast<TH1F*>(result_b->Clone());
218  if(fitCharmWithLight) {
219  sum->Add(result_cl);
220  }
221  else {
222  sum->Add(result_c);
223  sum->Add(result_l);
224  }
225  sum->Add(result_other);*/
226  sum->SetLineColor(kBlue);
227  sum->SetLineStyle(kSolid);
228  sum->SetLineWidth(3);
229  sum->SetMarkerSize(0);
230  sum->Draw("same");
231  sum->Draw("same,hist");
232  data->Draw("e,same");
233 
234  TLegend *leg1 = new TLegend(0.25, 0.5, 0.65, 0.87);
235  leg1->SetFillStyle(0);
236  leg1->SetBorderSize(0);
237  leg1->AddEntry(data, "Data", "ELP");
238  leg1->AddEntry(sum, "Total fit", "L");
239  leg1->AddEntry(result_b, "b events", "PL");
240  leg1->AddEntry(result_c, "c events", "PL");
241  leg1->AddEntry(result_l, "Light-parton events", "PL");
242  leg1->AddEntry(result_other, "Non-QCD events", "PL");
243  leg1->Draw();
244 
245  TLatex tla;
246  tla.SetTextSize(0.038);
247 
248 
249  TString cmslabel = "#font[62]{CMS} #scale[0.8]{#font[52]{Preliminary}}";
250  TString lumilabel = TString::Format("%1.1f",luminosity.Atof())+" fb^{-1} (13 TeV)";
251 
252  tla.DrawLatexNDC(0.18,0.93,cmslabel);
253  tla.SetTextFont(42);
254  tla.DrawLatexNDC(0.71,0.93,lumilabel);
255 
256  double qcd_b_fracafter, qcd_b_fracafter_err;
257  double qcd_c_fracafter, qcd_c_fracafter_err;
258  fit->GetResult(0, qcd_b_fracafter, qcd_b_fracafter_err);
259  fit->GetResult(2, qcd_c_fracafter, qcd_c_fracafter_err);
260  // float qcd_b_fracafter = fitter->Config().ParSettings(0).Value();
261  float other_fracafter = fitter->Config().ParSettings(1).Value();
262  float qcd_l_fracafter = -1.0;
263  float qcd_cl_fracafter = -1.0;
264  if(fitCharmWithLight) {
265  qcd_cl_fracafter = fitter->Config().ParSettings(2).Value();
266  }
267  else {
268  qcd_c_fracafter = fitter->Config().ParSettings(2).Value();
269  qcd_l_fracafter = fitter->Config().ParSettings(3).Value();
270  }
271  float qcd_b_ratio=qcd_b_fracafter/qcd_b_fracbefore;
272  float qcd_b_ratio_err=qcd_b_fracafter_err/qcd_b_fracbefore;
273  float other_ratio=other_fracafter/other_fracbefore;
274  float qcd_l_ratio, qcd_c_ratio, qcd_cl_ratio;
275  float qcd_l_ratio_err=0.0;
276  qcd_c_ratio = qcd_l_ratio = qcd_cl_ratio = -1;
277  if(fitCharmWithLight) {
278  qcd_cl_ratio = qcd_cl_fracafter/qcd_cl_fracbefore;
279  }
280  else {
281  qcd_c_ratio = qcd_c_fracafter/qcd_c_fracbefore;
282  qcd_l_ratio = qcd_l_fracafter/qcd_l_fracbefore;
283  }
284  float qcd_c_ratio_err=qcd_c_fracafter_err/qcd_c_fracbefore;
285 
286  std::cout << "fractions before fit:\n"
287  << "b: " << qcd_b_fracbefore << std::endl
288  << "c: " << qcd_c_fracbefore << std::endl
289  << "l: " << qcd_l_fracbefore << std::endl
290  << "cl: " << qcd_cl_fracbefore << std::endl
291  << "other: " << other_fracbefore << std::endl;
292  std::cout << "fractions after fit:\n"
293  << "b: " << qcd_b_fracafter << std::endl
294  << "c: " << qcd_c_fracafter << std::endl
295  << "l: " << qcd_l_fracafter << std::endl
296  << "cl: " << qcd_cl_fracafter << std::endl
297  << "other: " << other_fracafter << std::endl;
298  std::cout << "after/before fit ratio:\n"
299  << "b: " << qcd_b_ratio << " +/ " << qcd_b_ratio_err << std::endl
300  << "c: " << qcd_c_ratio << " +/ " << qcd_c_ratio_err << std::endl
301  << "l: " << qcd_l_ratio << std::endl
302  << "cl: " << qcd_cl_ratio << std::endl
303  << "other: " << other_ratio << std::endl;
304  c->Print(Form("plots/csvfit_%s.pdf", fitType.Data()));
305 
306  std::cout << "chi^2/ndof from TFractionFitter: " << fit->GetChisquare() << "/" << ndof << ", prob = " << fit->GetProb() << std::endl;
307  float chi2 =0.;
308  for(int ib =1; ib<=sum->GetNbinsX();ib++){
309  /*std::cout<<"Data yield: "<<data->GetBinContent(ib)<<std::endl;
310  std::cout<<"Sum yield: "<<sum->GetBinContent(ib)<<std::endl;
311  std::cout<<"Uncertainty on difference : "<<sqrt(pow(data->GetBinError(ib),2)+pow(sum->GetBinError(ib),2))<<std::endl;
312  std::cout<<"Sigma "<<(data->GetBinContent(ib) - sum->GetBinContent(ib))/sqrt(pow(sum->GetBinError(ib),2)+pow(data->GetBinError(ib),2))<<std::endl;*/
313  chi2+=pow((data->GetBinContent(ib) - sum->GetBinContent(ib))/sqrt(pow(sum->GetBinError(ib),2)+pow(data->GetBinError(ib),2)),2);
314  }
315  std::cout<<"Total chi2, by hand, is "<<chi2<<std::endl;
316 
317  // qcd_b->Scale(qcd_b_ratio);
318  // qcd_c->Scale(qcd_c_ratio);
319 
320  float chi2_b =0.;
321  for(int ib =1; ib<=sum->GetNbinsX();ib++){
322  /* std::cout<<"b result yield: "<<result_b->GetBinContent(ib)<<std::endl;
323  std::cout<<" b scaled prefit yield: "<<qcd_b->GetBinContent(ib)<<std::endl;
324  std::cout<<"Uncertainty on difference : "<<sqrt(pow(result_b->GetBinError(ib),2)+pow(qcd_b->GetBinError(ib),2))<<std::endl;
325  std::cout<<"Sigma "<<(result_b->GetBinContent(ib) - qcd_b->GetBinContent(ib))/sqrt(pow(qcd_b->GetBinError(ib),2)+pow(result_b->GetBinError(ib),2))<<std::endl;*/
326  chi2_b+=pow((result_b->GetBinContent(ib) - qcd_b->GetBinContent(ib))/sqrt(pow(qcd_b->GetBinError(ib),2)+pow(result_b->GetBinError(ib),2)),2);
327  }
328  std::cout<<"Total chi2 for b hist is "<<chi2_b<<std::endl;
329 
330  float chi2_c =0.;
331  for(int ib =1; ib<=sum->GetNbinsX();ib++){
332  /* std::cout<<"charm Result yield: "<<result_c->GetBinContent(ib)<<std::endl;
333  std::cout<<"charm Scaled prefit yield: "<<qcd_c->GetBinContent(ib)<<std::endl;
334  std::cout<<"Uncertainty on difference : "<<sqrt(pow(result_c->GetBinError(ib),2)+pow(qcd_c->GetBinError(ib),2))<<std::endl;
335  std::cout<<"Sigma "<<(result_c->GetBinContent(ib) - qcd_c->GetBinContent(ib))/sqrt(pow(qcd_c->GetBinError(ib),2)+pow(result_c->GetBinError(ib),2))<<std::endl;*/
336  chi2_c+=pow((result_c->GetBinContent(ib) - qcd_c->GetBinContent(ib))/sqrt(pow(qcd_c->GetBinError(ib),2)+pow(result_c->GetBinError(ib),2)),2);
337  }
338  std::cout<<"Total chi2 for c hist is "<<chi2_c<<std::endl;
339 
340  std::cout<<"Total chi2 for b,c, and sum"<<chi2_c+chi2_b+chi2<<std::endl;
341  // don't want to recreate files for variations
342  TFile *out;
343  if(fitType=="low_njet" && !excludeHighCSV) {
344  out = new TFile(Form("data/%s.root", fitType.Data()), "recreate");
345  TH1F *csv_weight = new TH1F("csv_weight", "csv_weight", 3, 0, 3);
346  csv_weight->SetBinContent(1, qcd_b_ratio);
347  csv_weight->SetBinError(1, qcd_b_ratio_err);
348  csv_weight->GetXaxis()->SetBinLabel(1, "b");
349  csv_weight->SetBinContent(2, qcd_c_ratio);
350  csv_weight->SetBinError(2, qcd_c_ratio_err);
351  csv_weight->GetXaxis()->SetBinLabel(2, "c");
352  csv_weight->SetBinContent(3, qcd_l_ratio);
353  csv_weight->SetBinError(3, qcd_l_ratio_err);
354  csv_weight->GetXaxis()->SetBinLabel(3, "l");
355  csv_weight->Write();
356  out->Close();
357  }
358  }
359 
360  return 0;
361 }
main
int main(int argc, char *argv[])
Definition: csv_fit.cxx:14
rpv::luminosity
TString luminosity
Definition: utilities_macros_rpv.hpp:3
efficiency_map.tla
tuple tla
Definition: efficiency_map.py:324
f
TFile * f
Definition: make_variations.cxx:31
efficiency_map.c
tuple c
Definition: efficiency_map.py:306
parse_card.data
tuple data
Definition: parse_card.py:260
Generated by   doxygen 1.8.11