CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
GenericTnPFitter.h
Go to the documentation of this file.
1 #include "RooGaussian.h"
2 #include "RooChebychev.h"
3 #include "RooVoigtian.h"
4 #include "RooExponential.h"
5 #include "RooRealVar.h"
6 #include "RooFormulaVar.h"
7 #include "RooDataHist.h"
8 #include "RooAddPdf.h"
9 #include "RooGlobalFunc.h"
10 #include "RooCategory.h"
11 #include "RooSimultaneous.h"
12 #include "RooChi2Var.h"
13 #include "TFile.h"
14 #include "TH1F.h"
15 #include "TH2F.h"
16 #include "TH3F.h"
17 #include "TProfile2D.h"
18 #include "TCanvas.h"
19 #include "RooPlot.h"
20 
21 namespace dqmTnP {
22 
23  class AbstractFitter {
24  protected:
25  RooRealVar mass;
26  RooRealVar mean;
27  double expectedMean;
28  RooRealVar sigma;
29  double expectedSigma;
30  RooRealVar efficiency;
31  RooRealVar nSignalAll;
32  RooFormulaVar nSignalPass;
33  RooFormulaVar nSignalFail;
34  RooRealVar nBackgroundFail;
35  RooRealVar nBackgroundPass;
36  RooCategory category;
37  RooSimultaneous simPdf;
38  RooDataHist* data;
39  double chi2;
40  bool verbose;
41 
42  public:
43  AbstractFitter(bool verbose_ = false)
44  : mass("mass", "mass", 0., 100., "GeV"),
45  mean("mean", "mean", 0., 100., "GeV"),
46  sigma("sigma", "sigma", 0., 100., "GeV"),
47  efficiency("efficiency", "efficiency", 0.5, 0.0, 1.0),
48  nSignalAll("nSignalAll", "nSignalAll", 0., 1e10),
49  nSignalPass("nSignalPass", "nSignalAll*efficiency", RooArgList(nSignalAll, efficiency)),
50  nSignalFail("nSignalFail", "nSignalAll*(1-efficiency)", RooArgList(nSignalAll, efficiency)),
51  nBackgroundFail("nBackgroundFail", "nBackgroundFail", 0., 1e10),
52  nBackgroundPass("nBackgroundPass", "nBackgroundPass", 0., 1e10),
53  category("category", "category"),
54  simPdf("simPdf", "simPdf", category),
55  data(nullptr),
56  verbose(verbose_) {
57  //turn on/off default messaging of roofit
58  RooMsgService::instance().setSilentMode(!verbose ? kTRUE : kFALSE);
59  for (int i = 0; i < RooMsgService::instance().numStreams(); i++) {
60  RooMsgService::instance().setStreamStatus(i, verbose ? kTRUE : kFALSE);
61  }
62  category.defineType("pass");
63  category.defineType("fail");
64  };
65  virtual ~AbstractFitter() = default;
66  ;
67  void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_) {
68  expectedMean = expectedMean_;
69  expectedSigma = expectedSigma_;
70  mass.setRange(massLow, massHigh);
71  mean.setRange(massLow, massHigh);
72  }
73  virtual void fit(TH1* num, TH1* den) = 0;
74  double getEfficiency() { return efficiency.getVal(); }
75  double getEfficiencyError() { return efficiency.getError(); }
76  double getChi2() { return chi2; }
77  void savePlot(const TString& name) {
78  using namespace RooFit;
79  RooPlot* frame = mass.frame(Name(name), Title("Failing and Passing Probe Distributions"));
80  data->plotOn(frame, Cut("category==category::pass"), LineColor(kGreen), MarkerColor(kGreen));
81  data->plotOn(frame, Cut("category==category::fail"), LineColor(kRed), MarkerColor(kRed));
82  simPdf.plotOn(frame, Slice(category, "pass"), ProjWData(category, *data), LineColor(kGreen));
83  simPdf.plotOn(frame, Slice(category, "fail"), ProjWData(category, *data), LineColor(kRed));
84  simPdf.paramOn(frame, Layout(0.58, 0.99, 0.99));
85  data->statOn(frame, Layout(0.70, 0.99, 0.5));
86  frame->Write();
87  delete frame;
88  }
89 
90  TString calculateEfficiency(
91  TH3* pass, TH3* all, int massDimension, TProfile2D*& eff, TProfile2D*& effChi2, const TString& plotName = "") {
92  //sort out the TAxis
93  TAxis *par1Axis, *par2Axis, *massAxis;
94  int par1C, par2C, massC;
95  if (massDimension == 1) {
96  massAxis = all->GetXaxis();
97  massC = 1;
98  par1Axis = all->GetYaxis();
99  par1C = all->GetXaxis()->GetNbins() + 2;
100  par2Axis = all->GetZaxis();
101  par2C = (all->GetXaxis()->GetNbins() + 2) * (all->GetYaxis()->GetNbins() + 2);
102  } else if (massDimension == 2) {
103  par1Axis = all->GetXaxis();
104  par1C = 1;
105  massAxis = all->GetYaxis();
106  massC = all->GetXaxis()->GetNbins() + 2;
107  par2Axis = all->GetZaxis();
108  par2C = (all->GetXaxis()->GetNbins() + 2) * (all->GetYaxis()->GetNbins() + 2);
109  } else if (massDimension == 3) {
110  par1Axis = all->GetXaxis();
111  par1C = 1;
112  par2Axis = all->GetYaxis();
113  par2C = all->GetXaxis()->GetNbins() + 2;
114  massAxis = all->GetZaxis();
115  massC = (all->GetXaxis()->GetNbins() + 2) * (all->GetYaxis()->GetNbins() + 2);
116  } else {
117  return "massDimension > 3 !, skipping...";
118  }
119 
120  //create eff and effChi2 TProfiles
121  if (!par1Axis || !par2Axis)
122  return "No par1Axis or par2Axis!";
123  if (par1Axis->GetXbins()->GetSize() == 0 && par2Axis->GetXbins()->GetSize() == 0) {
124  eff = new TProfile2D("efficiency",
125  "efficiency",
126  par1Axis->GetNbins(),
127  par1Axis->GetXmin(),
128  par1Axis->GetXmax(),
129  par2Axis->GetNbins(),
130  par2Axis->GetXmin(),
131  par2Axis->GetXmax());
132  } else if (par1Axis->GetXbins()->GetSize() == 0) {
133  eff = new TProfile2D("efficiency",
134  "efficiency",
135  par1Axis->GetNbins(),
136  par1Axis->GetXmin(),
137  par1Axis->GetXmax(),
138  par2Axis->GetNbins(),
139  par2Axis->GetXbins()->GetArray());
140  } else if (par2Axis->GetXbins()->GetSize() == 0) {
141  eff = new TProfile2D("efficiency",
142  "efficiency",
143  par1Axis->GetNbins(),
144  par1Axis->GetXbins()->GetArray(),
145  par2Axis->GetNbins(),
146  par2Axis->GetXmin(),
147  par2Axis->GetXmax());
148  } else {
149  eff = new TProfile2D("efficiency",
150  "efficiency",
151  par1Axis->GetNbins(),
152  par1Axis->GetXbins()->GetArray(),
153  par2Axis->GetNbins(),
154  par2Axis->GetXbins()->GetArray());
155  }
156  eff->SetTitle("");
157  eff->SetXTitle(par1Axis->GetTitle());
158  eff->SetYTitle(par2Axis->GetTitle());
159  eff->SetStats(kFALSE);
160  effChi2 = (TProfile2D*)eff->Clone("efficiencyChi2");
161  eff->SetZTitle("Efficiency");
162  eff->SetOption("colztexte");
163  eff->GetZaxis()->SetRangeUser(-0.001, 1.001);
164  effChi2->SetZTitle("Chi^2/NDF");
165  effChi2->SetOption("colztext");
166 
167  //create the 1D mass distribution container histograms
168  TH1D* all1D = (massAxis->GetXbins()->GetSize() == 0)
169  ? new TH1D("all1D", "all1D", massAxis->GetNbins(), massAxis->GetXmin(), massAxis->GetXmax())
170  : new TH1D("all1D", "all1D", massAxis->GetNbins(), massAxis->GetXbins()->GetArray());
171  auto* pass1D = (TH1D*)all1D->Clone("pass1D");
172 
173  //for each parameter bin fit the mass distributions
174  for (int par1 = 1; par1 <= par1Axis->GetNbins(); par1++) {
175  for (int par2 = 1; par2 <= par2Axis->GetNbins(); par2++) {
176  for (int mass = 1; mass <= massAxis->GetNbins(); mass++) {
177  int index = par1 * par1C + par2 * par2C + mass * massC;
178  all1D->SetBinContent(mass, all->GetBinContent(index));
179  pass1D->SetBinContent(mass, pass->GetBinContent(index));
180  }
181  fit(pass1D, all1D);
182  int index = par1 + par2 * (par1Axis->GetNbins() + 2);
183  eff->SetBinContent(index, getEfficiency());
184  eff->SetBinEntries(index, 1);
185  eff->SetBinError(index,
186  sqrt(getEfficiency() * getEfficiency() + getEfficiencyError() * getEfficiencyError()));
187  effChi2->SetBinContent(index, getChi2());
188  effChi2->SetBinEntries(index, 1);
189  if (plotName != "") {
190  savePlot(TString::Format("%s_%d_%d", plotName.Data(), par1, par2));
191  }
192  }
193  }
194  delete all1D;
195  delete pass1D;
196  return ""; //OK
197  }
198 
199  TString calculateEfficiency(
200  TH2* pass, TH2* all, int massDimension, TProfile*& eff, TProfile*& effChi2, const TString& plotName = "") {
201  //sort out the TAxis
202  TAxis *par1Axis, *massAxis;
203  int par1C, massC;
204  if (massDimension == 1) {
205  massAxis = all->GetXaxis();
206  massC = 1;
207  par1Axis = all->GetYaxis();
208  par1C = all->GetXaxis()->GetNbins() + 2;
209  } else if (massDimension == 2) {
210  par1Axis = all->GetXaxis();
211  par1C = 1;
212  massAxis = all->GetYaxis();
213  massC = all->GetXaxis()->GetNbins() + 2;
214  } else {
215  return "massDimension > 2 !, skipping...";
216  }
217 
218  //create eff and effChi2 TProfiles
219  if (!par1Axis)
220  return "No par1Axis!";
221  eff =
222  (par1Axis->GetXbins()->GetSize() == 0)
223  ? new TProfile("efficiency", "efficiency", par1Axis->GetNbins(), par1Axis->GetXmin(), par1Axis->GetXmax())
224  : new TProfile("efficiency", "efficiency", par1Axis->GetNbins(), par1Axis->GetXbins()->GetArray());
225  eff->SetTitle("");
226  eff->SetXTitle(par1Axis->GetTitle());
227  eff->SetLineColor(2);
228  eff->SetLineWidth(2);
229  eff->SetMarkerStyle(20);
230  eff->SetMarkerSize(0.8);
231  eff->SetStats(kFALSE);
232  effChi2 = (TProfile*)eff->Clone("efficiencyChi2");
233  eff->SetYTitle("Efficiency");
234  eff->SetOption("PE");
235  eff->GetYaxis()->SetRangeUser(-0.001, 1.001);
236  effChi2->SetYTitle("Chi^2/NDF");
237  effChi2->SetOption("HIST");
238 
239  //create the 1D mass distribution container histograms
240  TH1D* all1D = (massAxis->GetXbins()->GetSize() == 0)
241  ? new TH1D("all1D", "all1D", massAxis->GetNbins(), massAxis->GetXmin(), massAxis->GetXmax())
242  : new TH1D("all1D", "all1D", massAxis->GetNbins(), massAxis->GetXbins()->GetArray());
243  auto* pass1D = (TH1D*)all1D->Clone("pass1D");
244 
245  //for each parameter bin fit the mass distributions
246  for (int par1 = 1; par1 <= par1Axis->GetNbins(); par1++) {
247  for (int mass = 1; mass <= massAxis->GetNbins(); mass++) {
248  int index = par1 * par1C + mass * massC;
249  all1D->SetBinContent(mass, all->GetBinContent(index));
250  pass1D->SetBinContent(mass, pass->GetBinContent(index));
251  }
252  fit(pass1D, all1D);
253  int index = par1;
254  eff->SetBinContent(index, getEfficiency());
255  eff->SetBinEntries(index, 1);
256  eff->SetBinError(index, sqrt(getEfficiency() * getEfficiency() + getEfficiencyError() * getEfficiencyError()));
257  effChi2->SetBinContent(index, getChi2());
258  effChi2->SetBinEntries(index, 1);
259  if (plotName != "") {
260  savePlot(TString::Format("%s_%d", plotName.Data(), par1));
261  }
262  }
263  delete all1D;
264  delete pass1D;
265  return ""; //OK
266  }
267  };
268 
269  //concrete fitter: Gaussian signal plus linear background
270  class GaussianPlusLinearFitter : public AbstractFitter {
271  protected:
272  RooGaussian gaussian;
273  RooRealVar slopeFail;
274  RooChebychev linearFail;
275  RooRealVar slopePass;
276  RooChebychev linearPass;
277  RooAddPdf pdfFail;
278  RooAddPdf pdfPass;
279 
280  public:
281  GaussianPlusLinearFitter(bool verbose = false)
282  : AbstractFitter(verbose),
283  gaussian("gaussian", "gaussian", mass, mean, sigma),
284  slopeFail("slopeFail", "slopeFail", 0., -1., 1.),
285  linearFail("linearFail", "linearFail", mass, slopeFail),
286  slopePass("slopePass", "slopePass", 0., -1., 1.),
287  linearPass("linearPass", "linearPass", mass, slopePass),
288  pdfFail("pdfFail", "pdfFail", RooArgList(gaussian, linearFail), RooArgList(nSignalFail, nBackgroundFail)),
289  pdfPass("pdfPass", "pdfPass", RooArgList(gaussian, linearPass), RooArgList(nSignalPass, nBackgroundPass)) {
290  simPdf.addPdf(pdfFail, "fail");
291  simPdf.addPdf(pdfPass, "pass");
292  };
293  ~GaussianPlusLinearFitter() override = default;
294  ;
295  void fit(TH1* pass, TH1* all) override {
296  using namespace RooFit;
297  all->Add(pass, -1);
298  TH1*& fail = all;
299  if (!data)
300  delete data;
301  data = new RooDataHist("data", "data", mass, Index(category), Import("fail", *fail), Import("pass", *pass));
302  if (pass->Integral() + fail->Integral() < 5) {
303  efficiency.setVal(0.5);
304  efficiency.setError(0.5);
305  chi2 = 0;
306  return;
307  }
308  mean.setVal(expectedMean);
309  sigma.setVal(expectedSigma);
310  efficiency.setVal(pass->Integral() / (pass->Integral() + fail->Integral()));
311  nSignalAll.setVal(0.5 * (fail->Integral() + pass->Integral()));
312  nBackgroundFail.setVal(0.5 * fail->Integral());
313  nBackgroundPass.setVal(0.5 * pass->Integral());
314  slopeFail.setVal(0.);
315  slopePass.setVal(0.);
316  if (verbose) {
317  simPdf.fitTo(*data);
318  } else {
319  simPdf.fitTo(*data, Verbose(kFALSE), PrintLevel(-1), Warnings(kFALSE), PrintEvalErrors(-1));
320  }
321  RooDataHist dataFail("fail", "fail", mass, fail);
322  RooDataHist dataPass("pass", "pass", mass, pass);
323  chi2 = (RooChi2Var("chi2Fail", "chi2Fail", pdfFail, dataFail, DataError(RooAbsData::Poisson)).getVal() +
324  RooChi2Var("chi2Pass", "chi2Pass", pdfPass, dataPass, DataError(RooAbsData::Poisson)).getVal()) /
325  (2 * pass->GetNbinsX() - 8);
326  if (chi2 > 3) {
327  efficiency.setVal(0.5);
328  efficiency.setError(0.5);
329  }
330  }
331  };
332 
333  //concrete fitter: voigtian signal plus exponential background
334  class VoigtianPlusExponentialFitter : public AbstractFitter {
335  protected:
336  RooRealVar width;
337  RooVoigtian voigtian;
338  RooRealVar slopeFail;
339  RooExponential exponentialFail;
340  RooRealVar slopePass;
341  RooExponential exponentialPass;
342  RooAddPdf pdfFail;
343  RooAddPdf pdfPass;
344 
345  public:
346  VoigtianPlusExponentialFitter(bool verbose = false)
347  : AbstractFitter(verbose),
348  width("width", "width", 2.5, "GeV"),
349  voigtian("voigtian", "voigtian", mass, mean, width, sigma),
350  slopeFail("slopeFail", "slopeFail", 0., -1., 0.),
351  exponentialFail("linearFail", "linearFail", mass, slopeFail),
352  slopePass("slopePass", "slopePass", 0., -1., 0.),
353  exponentialPass("linearPass", "linearPass", mass, slopePass),
354  pdfFail(
355  "pdfFail", "pdfFail", RooArgList(voigtian, exponentialFail), RooArgList(nSignalFail, nBackgroundFail)),
356  pdfPass(
357  "pdfPass", "pdfPass", RooArgList(voigtian, exponentialPass), RooArgList(nSignalPass, nBackgroundPass)) {
358  width.setConstant(kTRUE);
359  simPdf.addPdf(pdfFail, "fail");
360  simPdf.addPdf(pdfPass, "pass");
361  };
362  ~VoigtianPlusExponentialFitter() override = default;
363  ;
364  void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_, double width_) {
365  expectedMean = expectedMean_;
366  expectedSigma = expectedSigma_;
367  mass.setRange(massLow, massHigh);
368  mean.setRange(massLow, massHigh);
369  width.setVal(width_);
370  }
371  void fit(TH1* pass, TH1* all) override {
372  using namespace RooFit;
373  all->Add(pass, -1);
374  TH1*& fail = all;
375  if (!data)
376  delete data;
377  data = new RooDataHist("data", "data", mass, Index(category), Import("fail", *fail), Import("pass", *pass));
378  if (pass->Integral() + fail->Integral() < 5) {
379  efficiency.setVal(0.5);
380  efficiency.setError(0.5);
381  chi2 = 0;
382  return;
383  }
384  mean.setVal(expectedMean);
385  sigma.setVal(expectedSigma);
386  efficiency.setVal(pass->Integral() / (pass->Integral() + fail->Integral()));
387  nSignalAll.setVal(0.5 * (fail->Integral() + pass->Integral()));
388  nBackgroundFail.setVal(0.5 * fail->Integral());
389  nBackgroundPass.setVal(0.5 * pass->Integral());
390  slopeFail.setVal(0.);
391  slopePass.setVal(0.);
392  if (verbose) {
393  simPdf.fitTo(*data);
394  } else {
395  simPdf.fitTo(*data, Verbose(kFALSE), PrintLevel(-1), Warnings(kFALSE), PrintEvalErrors(-1));
396  }
397  RooDataHist dataFail("fail", "fail", mass, fail);
398  RooDataHist dataPass("pass", "pass", mass, pass);
399  chi2 = (RooChi2Var("chi2Fail", "chi2Fail", pdfFail, dataFail, DataError(RooAbsData::Poisson)).getVal() +
400  RooChi2Var("chi2Pass", "chi2Pass", pdfPass, dataPass, DataError(RooAbsData::Poisson)).getVal()) /
401  (2 * all->GetNbinsX() - 8);
402  if (chi2 > 3) {
403  efficiency.setVal(0.5);
404  efficiency.setError(0.5);
405  }
406  }
407  };
408 
409 } //namespace dqmTnP
mps_fire.i
i
Definition: mps_fire.py:428
instance
static PFTauRenderPlugin instance
Definition: PFTauRenderPlugin.cc:70
dqmTnP::AbstractFitter::setup
void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_)
Definition: GenericTnPFitter.h:67
dqmTnP::AbstractFitter::savePlot
void savePlot(const TString &name)
Definition: GenericTnPFitter.h:77
L1TEGammaDiff_cfi.plotName
string plotName
Definition: L1TEGammaDiff_cfi.py:20
dqmTnP::VoigtianPlusExponentialFitter::setup
void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_, double width_)
Definition: GenericTnPFitter.h:364
EcalTangentSkim_cfg.Verbose
tuple Verbose
Definition: EcalTangentSkim_cfg.py:30
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
dqmTnP::AbstractFitter::~AbstractFitter
virtual ~AbstractFitter()=default
python.cmstools.all
def all
workaround iterator generators for ROOT classes
Definition: cmstools.py:25
dqmTnP::AbstractFitter::calculateEfficiency
TString calculateEfficiency(TH3 *pass, TH3 *all, int massDimension, TProfile2D *&eff, TProfile2D *&effChi2, const TString &plotName="")
Definition: GenericTnPFitter.h:90
dqmTnP::GaussianPlusLinearFitter::GaussianPlusLinearFitter
GaussianPlusLinearFitter(bool verbose=false)
Definition: GenericTnPFitter.h:281
dqmTnP::VoigtianPlusExponentialFitter::fit
void fit(TH1 *pass, TH1 *all) override
Definition: GenericTnPFitter.h:371
dqmTnP::VoigtianPlusExponentialFitter::VoigtianPlusExponentialFitter
VoigtianPlusExponentialFitter(bool verbose=false)
Definition: GenericTnPFitter.h:346
dqmTnP::GaussianPlusLinearFitter::fit
void fit(TH1 *pass, TH1 *all) override
Definition: GenericTnPFitter.h:295
dqmTnP::AbstractFitter::AbstractFitter
AbstractFitter(bool verbose_=false)
Definition: GenericTnPFitter.h:43
dqmTnP::AbstractFitter::calculateEfficiency
TString calculateEfficiency(TH2 *pass, TH2 *all, int massDimension, TProfile *&eff, TProfile *&effChi2, const TString &plotName="")
Definition: GenericTnPFitter.h:199
dqmTnP::AbstractFitter::fit
virtual void fit(TH1 *num, TH1 *den)=0
dqmTnP::AbstractFitter::simPdf
RooSimultaneous simPdf
Definition: GenericTnPFitter.h:37
Generated for CMSSW Reference Manual by doxygen 1.8.5