CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros 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(0),
56  verbose(verbose_)
57  {
58  //turn on/off default messaging of roofit
59  RooMsgService::instance().setSilentMode( !verbose?kTRUE:kFALSE );
60  for(int i=0; i<RooMsgService::instance().numStreams(); i++){
61  RooMsgService::instance().setStreamStatus( i, verbose?kTRUE:kFALSE );
62  }
63  category.defineType("pass");
64  category.defineType("fail");
65  };
66  virtual ~AbstractFitter(){};
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(TH3 *pass, TH3 *all, int massDimension, TProfile2D* &eff, TProfile2D* &effChi2, TString plotName=""){
91  //sort out the TAxis
92  TAxis *par1Axis, *par2Axis, *massAxis;
93  int par1C, par2C, massC;
94  if(massDimension==1){
95  massAxis = all->GetXaxis();
96  massC = 1;
97  par1Axis = all->GetYaxis();
98  par1C = all->GetXaxis()->GetNbins()+2;
99  par2Axis = all->GetZaxis();
100  par2C = (all->GetXaxis()->GetNbins()+2)*(all->GetYaxis()->GetNbins()+2);
101  }else if(massDimension==2){
102  par1Axis = all->GetXaxis();
103  par1C = 1;
104  massAxis = all->GetYaxis();
105  massC = all->GetXaxis()->GetNbins()+2;
106  par2Axis = all->GetZaxis();
107  par2C = (all->GetXaxis()->GetNbins()+2)*(all->GetYaxis()->GetNbins()+2);
108  }else if(massDimension==3){
109  par1Axis = all->GetXaxis();
110  par1C = 1;
111  par2Axis = all->GetYaxis();
112  par2C = all->GetXaxis()->GetNbins()+2;
113  massAxis = all->GetZaxis();
114  massC = (all->GetXaxis()->GetNbins()+2)*(all->GetYaxis()->GetNbins()+2);
115  }else{
116  return "massDimension > 3 !, skipping...";
117  }
118 
119  //create eff and effChi2 TProfiles
120  if(!par1Axis || !par2Axis) return "No par1Axis or par2Axis!";
121  if( par1Axis->GetXbins()->GetSize()==0 && par2Axis->GetXbins()->GetSize()==0 ){
122  eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXmin(),par1Axis->GetXmax(),par2Axis->GetNbins(),par2Axis->GetXmin(),par2Axis->GetXmax());
123  }else if( par1Axis->GetXbins()->GetSize()==0 ){
124  eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXmin(),par1Axis->GetXmax(),par2Axis->GetNbins(),par2Axis->GetXbins()->GetArray());
125  }else if( par2Axis->GetXbins()->GetSize()==0 ){
126  eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXbins()->GetArray(),par2Axis->GetNbins(),par2Axis->GetXmin(),par2Axis->GetXmax());
127  }else{
128  eff = new TProfile2D("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXbins()->GetArray(),par2Axis->GetNbins(),par2Axis->GetXbins()->GetArray());
129  }
130  eff->SetTitle("");
131  eff->SetXTitle( par1Axis->GetTitle() );
132  eff->SetYTitle( par2Axis->GetTitle() );
133  eff->SetStats(kFALSE);
134  effChi2 = (TProfile2D*)eff->Clone("efficiencyChi2");
135  eff->SetZTitle("Efficiency");
136  eff->SetOption("colztexte");
137  eff->GetZaxis()->SetRangeUser(-0.001,1.001);
138  effChi2->SetZTitle("Chi^2/NDF");
139  effChi2->SetOption("colztext");
140 
141  //create the 1D mass distribution container histograms
142  TH1D* all1D = (massAxis->GetXbins()->GetSize()==0) ?
143  new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXmin(),massAxis->GetXmax()):
144  new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXbins()->GetArray());
145  TH1D* pass1D = (TH1D *)all1D->Clone("pass1D");
146 
147  //for each parameter bin fit the mass distributions
148  for(int par1=1; par1<=par1Axis->GetNbins(); par1++){
149  for(int par2=1; par2<=par2Axis->GetNbins(); par2++){
150  for(int mass=1; mass<=massAxis->GetNbins(); mass++){
151  int index = par1*par1C + par2*par2C + mass*massC;
152  all1D->SetBinContent(mass,all->GetBinContent(index));
153  pass1D->SetBinContent(mass,pass->GetBinContent(index));
154  }
155  fit( pass1D, all1D );
156  int index = par1 + par2*(par1Axis->GetNbins()+2);
157  eff->SetBinContent( index, getEfficiency() );
158  eff->SetBinEntries( index, 1 );
159  eff->SetBinError( index, sqrt( getEfficiency()*getEfficiency() + getEfficiencyError()*getEfficiencyError() ) );
160  effChi2->SetBinContent( index, getChi2() );
161  effChi2->SetBinEntries( index, 1 );
162  if(plotName!=""){
163  savePlot( TString::Format("%s_%d_%d",plotName.Data(),par1,par2) );
164  }
165  }
166  }
167  delete all1D;
168  delete pass1D;
169  return "";//OK
170  }
171 
172  TString calculateEfficiency(TH2 *pass, TH2 *all, int massDimension, TProfile* &eff, TProfile* &effChi2, TString plotName=""){
173  //sort out the TAxis
174  TAxis *par1Axis, *massAxis;
175  int par1C, massC;
176  if(massDimension==1){
177  massAxis = all->GetXaxis();
178  massC = 1;
179  par1Axis = all->GetYaxis();
180  par1C = all->GetXaxis()->GetNbins()+2;
181  }else if(massDimension==2){
182  par1Axis = all->GetXaxis();
183  par1C = 1;
184  massAxis = all->GetYaxis();
185  massC = all->GetXaxis()->GetNbins()+2;
186  }else{
187  return "massDimension > 2 !, skipping...";
188  }
189 
190  //create eff and effChi2 TProfiles
191  if(!par1Axis) return "No par1Axis!";
192  eff = (par1Axis->GetXbins()->GetSize()==0)?
193  new TProfile("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXmin(),par1Axis->GetXmax()):
194  new TProfile("efficiency","efficiency",par1Axis->GetNbins(),par1Axis->GetXbins()->GetArray());
195  eff->SetTitle("");
196  eff->SetXTitle( par1Axis->GetTitle() );
197  eff->SetLineColor(2);
198  eff->SetLineWidth(2);
199  eff->SetMarkerStyle(20);
200  eff->SetMarkerSize(0.8);
201  eff->SetStats(kFALSE);
202  effChi2 = (TProfile*)eff->Clone("efficiencyChi2");
203  eff->SetYTitle("Efficiency");
204  eff->SetOption("PE");
205  eff->GetYaxis()->SetRangeUser(-0.001,1.001);
206  effChi2->SetYTitle("Chi^2/NDF");
207  effChi2->SetOption("HIST");
208 
209  //create the 1D mass distribution container histograms
210  TH1D * all1D = (massAxis->GetXbins()->GetSize()==0) ?
211  new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXmin(),massAxis->GetXmax()):
212  new TH1D("all1D","all1D",massAxis->GetNbins(),massAxis->GetXbins()->GetArray());
213  TH1D * pass1D = (TH1D *)all1D->Clone("pass1D");
214 
215  //for each parameter bin fit the mass distributions
216  for(int par1=1; par1<=par1Axis->GetNbins(); par1++){
217  for(int mass=1; mass<=massAxis->GetNbins(); mass++){
218  int index = par1*par1C + mass*massC;
219  all1D->SetBinContent(mass,all->GetBinContent(index));
220  pass1D->SetBinContent(mass,pass->GetBinContent(index));
221  }
222  fit( pass1D, all1D );
223  int index = par1;
224  eff->SetBinContent( index, getEfficiency() );
225  eff->SetBinEntries( index, 1 );
226  eff->SetBinError( index, sqrt( getEfficiency()*getEfficiency() + getEfficiencyError()*getEfficiencyError() ) );
227  effChi2->SetBinContent( index, getChi2() );
228  effChi2->SetBinEntries( index, 1 );
229  if(plotName!=""){
230  savePlot( TString::Format("%s_%d",plotName.Data(),par1) );
231  }
232  }
233  delete all1D;
234  delete pass1D;
235  return "";//OK
236  }
237 };
238 
239 //concrete fitter: Gaussian signal plus linear background
240 class GaussianPlusLinearFitter: public AbstractFitter{
241  protected:
242  RooGaussian gaussian;
243  RooRealVar slopeFail;
244  RooChebychev linearFail;
245  RooRealVar slopePass;
246  RooChebychev linearPass;
247  RooAddPdf pdfFail;
248  RooAddPdf pdfPass;
249  public:
250  GaussianPlusLinearFitter(bool verbose = false):
251  AbstractFitter(verbose),
252  gaussian("gaussian","gaussian",mass,mean,sigma),
253  slopeFail("slopeFail","slopeFail",0.,-1.,1.),
254  linearFail("linearFail","linearFail",mass,slopeFail),
255  slopePass("slopePass","slopePass",0.,-1.,1.),
256  linearPass("linearPass","linearPass",mass,slopePass),
257  pdfFail("pdfFail","pdfFail", RooArgList(gaussian,linearFail), RooArgList(nSignalFail,nBackgroundFail)),
258  pdfPass("pdfPass","pdfPass", RooArgList(gaussian,linearPass), RooArgList(nSignalPass,nBackgroundPass))
259  {
260  simPdf.addPdf(pdfFail,"fail");
261  simPdf.addPdf(pdfPass,"pass");
262  };
263  ~GaussianPlusLinearFitter(){};
264  void fit(TH1* pass, TH1* all){
265  using namespace RooFit;
266  all->Add(pass,-1);
267  TH1* &fail = all;
268  if(!data) delete data;
269  data = new RooDataHist("data", "data", mass, Index(category), Import("fail",*fail), Import("pass",*pass) );
270  if(pass->Integral()+fail->Integral() < 5){
271  efficiency.setVal(0.5);
272  efficiency.setError(0.5);
273  chi2 = 0;
274  return;
275  }
276  mean.setVal(expectedMean);
277  sigma.setVal(expectedSigma);
278  efficiency.setVal(pass->Integral()/(pass->Integral()+fail->Integral()));
279  nSignalAll.setVal(0.5*(fail->Integral()+pass->Integral()));
280  nBackgroundFail.setVal(0.5*fail->Integral());
281  nBackgroundPass.setVal(0.5*pass->Integral());
282  slopeFail.setVal(0.);
283  slopePass.setVal(0.);
284  if(verbose){
285  simPdf.fitTo( *data );
286  }else{
287  simPdf.fitTo( *data, Verbose(kFALSE), PrintLevel(-1), Warnings(kFALSE), PrintEvalErrors(-1) );
288  }
289  RooDataHist dataFail("fail", "fail", mass, fail );
290  RooDataHist dataPass("pass", "pass", mass, pass );
291  chi2 = ( RooChi2Var("chi2Fail","chi2Fail",pdfFail,dataFail,DataError(RooAbsData::Poisson)).getVal()
292  +RooChi2Var("chi2Pass","chi2Pass",pdfPass,dataPass,DataError(RooAbsData::Poisson)).getVal() )/(2*pass->GetNbinsX()-8);
293  if(chi2>3){
294  efficiency.setVal(0.5);
295  efficiency.setError(0.5);
296  }
297  }
298 };
299 
300 //concrete fitter: voigtian signal plus exponential background
301 class VoigtianPlusExponentialFitter: public AbstractFitter{
302  protected:
303  RooRealVar width;
304  RooVoigtian voigtian;
305  RooRealVar slopeFail;
306  RooExponential exponentialFail;
307  RooRealVar slopePass;
308  RooExponential exponentialPass;
309  RooAddPdf pdfFail;
310  RooAddPdf pdfPass;
311  public:
312  VoigtianPlusExponentialFitter(bool verbose = false):
313  AbstractFitter(verbose),
314  width("width","width",2.5,"GeV"),
315  voigtian("voigtian","voigtian",mass,mean,width,sigma),
316  slopeFail("slopeFail","slopeFail",0.,-1.,0.),
317  exponentialFail("linearFail","linearFail",mass,slopeFail),
318  slopePass("slopePass","slopePass",0.,-1.,0.),
319  exponentialPass("linearPass","linearPass",mass,slopePass),
320  pdfFail("pdfFail","pdfFail", RooArgList(voigtian,exponentialFail), RooArgList(nSignalFail,nBackgroundFail)),
321  pdfPass("pdfPass","pdfPass", RooArgList(voigtian,exponentialPass), RooArgList(nSignalPass,nBackgroundPass))
322  {
323  width.setConstant(kTRUE);
324  simPdf.addPdf(pdfFail,"fail");
325  simPdf.addPdf(pdfPass,"pass");
326  };
327  ~VoigtianPlusExponentialFitter(){};
328  void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_, double width_){
329  expectedMean = expectedMean_;
330  expectedSigma = expectedSigma_;
331  mass.setRange(massLow,massHigh);
332  mean.setRange(massLow,massHigh);
333  width.setVal(width_);
334  }
335  void fit(TH1* pass, TH1* all){
336  using namespace RooFit;
337  all->Add(pass,-1);
338  TH1* &fail = all;
339  if(!data) delete data;
340  data = new RooDataHist("data", "data", mass, Index(category), Import("fail",*fail), Import("pass",*pass) );
341  if(pass->Integral()+fail->Integral() < 5){
342  efficiency.setVal(0.5);
343  efficiency.setError(0.5);
344  chi2 = 0;
345  return;
346  }
347  mean.setVal(expectedMean);
348  sigma.setVal(expectedSigma);
349  efficiency.setVal(pass->Integral()/(pass->Integral()+fail->Integral()));
350  nSignalAll.setVal(0.5*(fail->Integral()+pass->Integral()));
351  nBackgroundFail.setVal(0.5*fail->Integral());
352  nBackgroundPass.setVal(0.5*pass->Integral());
353  slopeFail.setVal(0.);
354  slopePass.setVal(0.);
355  if(verbose){
356  simPdf.fitTo( *data );
357  }else{
358  simPdf.fitTo( *data, Verbose(kFALSE), PrintLevel(-1), Warnings(kFALSE), PrintEvalErrors(-1) );
359  }
360  RooDataHist dataFail("fail", "fail", mass, fail );
361  RooDataHist dataPass("pass", "pass", mass, pass );
362  chi2 = ( RooChi2Var("chi2Fail","chi2Fail",pdfFail,dataFail,DataError(RooAbsData::Poisson)).getVal()
363  +RooChi2Var("chi2Pass","chi2Pass",pdfPass,dataPass,DataError(RooAbsData::Poisson)).getVal() )/(2*all->GetNbinsX()-8);
364  if(chi2>3){
365  efficiency.setVal(0.5);
366  efficiency.setError(0.5);
367  }
368  }
369 };
370 
371 }//namespace dqmTnP
i
int i
Definition: DBlmapReader.cc:9
create_public_lumi_plots.width
int width
Definition: create_public_lumi_plots.py:1009
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
dqmTnP::VoigtianPlusExponentialFitter::setup
void setup(double expectedMean_, double massLow, double massHigh, double expectedSigma_, double width_)
Definition: GenericTnPFitter.h:328
dqmTnP::AbstractFitter::calculateEfficiency
TString calculateEfficiency(TH3 *pass, TH3 *all, int massDimension, TProfile2D *&eff, TProfile2D *&effChi2, TString plotName="")
Definition: GenericTnPFitter.h:90
EcalTangentSkim_cfg.Verbose
tuple Verbose
Definition: EcalTangentSkim_cfg.py:24
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
dqmTnP::GaussianPlusLinearFitter::GaussianPlusLinearFitter
GaussianPlusLinearFitter(bool verbose=false)
Definition: GenericTnPFitter.h:250
dqmTnP::VoigtianPlusExponentialFitter::VoigtianPlusExponentialFitter
VoigtianPlusExponentialFitter(bool verbose=false)
Definition: GenericTnPFitter.h:312
dqmTnP::AbstractFitter::AbstractFitter
AbstractFitter(bool verbose_=false)
Definition: GenericTnPFitter.h:43
dqmTnP::AbstractFitter::fit
virtual void fit(TH1 *num, TH1 *den)=0
dqmTnP::AbstractFitter::calculateEfficiency
TString calculateEfficiency(TH2 *pass, TH2 *all, int massDimension, TProfile *&eff, TProfile *&effChi2, TString plotName="")
Definition: GenericTnPFitter.h:172
dqmTnP::AbstractFitter::simPdf
RooSimultaneous simPdf
Definition: GenericTnPFitter.h:37
dqmTnP::GaussianPlusLinearFitter::fit
void fit(TH1 *pass, TH1 *all)
Definition: GenericTnPFitter.h:264
Generated for CMSSW Reference Manual by doxygen 1.8.5