10 #include "TGraphAsymmErrors.h"
13 #include "Math/QuantFuncMathCore.h"
14 #include "Roo1DTable.h"
15 #include "RooAbsDataStore.h"
16 #include "RooAbsReal.h"
17 #include "RooAddPdf.h"
18 #include "RooBinning.h"
19 #include "RooBinningCategory.h"
20 #include "RooCategory.h"
21 #include "RooChebychev.h"
22 #include "RooDataHist.h"
23 #include "RooDataSet.h"
24 #include "RooEfficiency.h"
25 #include "RooExtendPdf.h"
26 #include "RooFitResult.h"
27 #include "RooFormulaVar.h"
28 #include "RooGaussian.h"
29 #include "RooGenericPdf.h"
30 #include "RooGlobalFunc.h"
31 #include "RooLinkedListIter.h"
32 #include "RooMappedCategory.h"
33 #include "RooMinimizer.h"
34 #include "RooMinuit.h"
35 #include "RooMsgService.h"
36 #include "RooMultiCategory.h"
37 #include "RooNLLVar.h"
38 #include "RooNumIntConfig.h"
40 #include "RooProdPdf.h"
41 #include "RooProfileLL.h"
42 #include "RooRealVar.h"
43 #include "RooSimultaneous.h"
44 #include "RooThresholdCategory.h"
46 #include "RooWorkspace.h"
49 using namespace RooFit;
52 inputTree =
new TChain((inputDirectoryName+
"/"+inputTreeName).c_str());
53 for(
size_t i=0;
i<inputFileNames.size();
i++){
54 inputTree->Add(inputFileNames[
i].c_str());
56 outputFile =
new TFile(outputFileName.c_str(),
"recreate");
57 outputDirectory =
outputFile->mkdir(inputDirectoryName.c_str());
59 saveWorkspace = saveWorkspace_;
61 floatShapeParameters = floatShapeParameters_;
64 if(!floatShapeParameters && fixVars.empty())
std::cout <<
"TagProbeFitter: " <<
"You wnat to fix some variables but do not specify them!";
66 gROOT->SetStyle(
"Plain");
67 gStyle->SetTitleFillColor(0);
68 gStyle->SetPalette(1);
69 gStyle->SetOptStat(0);
70 gStyle->SetPaintTextFormat(
".2f");
76 doSaveDistributionsPlot =
true;
79 RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1
e-13);
80 RooAbsReal::defaultIntegratorConfig()->setEpsRel(1
e-13);
99 RooRealVar
temp(name.c_str(), title.c_str(), low, hi, units.c_str());
100 temp.setBins(5000,
"cache");
106 RooCategory*
c = (RooCategory*) parameterParser.factory(expression.c_str());
110 c->SetName(name.c_str());
111 c->SetTitle(title.c_str());
117 expressionVars.push_back(make_pair(make_pair(expressionName,title), make_pair(expression,arguments)));
123 thresholdCategories.push_back(make_pair(make_pair(categoryName,title), make_pair(varName,cutValue)));
129 pdfs[
name] = pdfCommands;
140 string TagProbeFitter::calculateEfficiency(
string dirName,
const std::vector<string>& effCats,
const std::vector<string>& effStates, vector<string>& unbinnedVariables,
map<
string, vector<double> >& binnedReals,
map<
string, std::vector<string> >& binnedCategories, vector<string>& binToPDFmap){
142 outputDirectory->cd();
144 gDirectory->mkdir(dirName.c_str())->cd();
149 for(vector<string>::iterator
v=unbinnedVariables.begin();
v!=unbinnedVariables.end();
v++){
150 dataVars.addClone(
variables[
v->c_str()],
true);
151 if (binnedFit && (
v == unbinnedVariables.begin())) {
152 ((RooRealVar&)dataVars[
v->c_str()]).setBins(massBins);
156 RooArgSet binnedVariables;
157 RooArgSet binCategories;
158 for(
map<
string, vector<double> >::iterator
v=binnedReals.begin();
v!=binnedReals.end();
v++){
159 TString
name =
v->first;
160 if (
variables.find(name) == 0) {
cerr <<
"Binned variable '"<<name<<
"' not found." << endl;
return "Error"; }
161 binnedVariables.addClone(
variables[name]);
162 ((RooRealVar&)binnedVariables[name]).setBinning( RooBinning(
v->second.size()-1, &
v->second[0]) );
163 binCategories.addClone( RooBinningCategory(name+
"_bins", name+
"_bins", (RooRealVar&)binnedVariables[name]) );
165 dataVars.addClone(binnedVariables,
true);
169 RooArgSet mappedCategories;
170 for(
map<
string, vector<string> >::iterator
v=binnedCategories.begin();
v!=binnedCategories.end();
v++){
171 TString name =
v->first;
172 if (
variables.find(name) == 0) {
cerr <<
"Binned category '"<<name<<
"' not found." << endl;
return "Error"; }
174 mappedCategories.addClone(RooMappedCategory(name+
"_bins", name+
"_bins", (RooCategory&)categories[name]));
175 for(
unsigned int i = 0;
i<
v->second.size();
i++){
176 ((RooMappedCategory&)mappedCategories[name+
"_bins"]).map(
v->second[
i].c_str(), name+
"_"+TString(
v->second[
i].c_str()).ReplaceAll(
",",
"_"));
179 dataVars.addClone(categories,
true);
182 for (vector<string>::const_iterator effCat = effCats.begin(); effCat != effCats.end(); ++effCat) {
183 if (
variables.find(effCat->c_str()) != 0) {
184 dataVars.addClone(
variables[effCat->c_str()],
true);
189 for(vector<pair<pair<string,string>, pair<
string, vector<string> > > >::const_iterator
ev = expressionVars.begin(), eve = expressionVars.end();
ev != eve; ++
ev){
190 for (vector<string>::const_iterator it =
ev->second.second.begin(), ed =
ev->second.second.end(); it != ed; ++it) {
196 for(vector<pair<pair<string,string>, pair<string, double> > >::const_iterator tc = thresholdCategories.begin(), tce = thresholdCategories.end(); tc != tce; ++tc){
197 if (
variables.find(tc->second.first.c_str())) dataVars.addClone(
variables[tc->second.first.c_str()],
true);
202 RooDataSet
data(
"data",
"data", inputTree,
204 "", (weightVar.empty() ? 0 : weightVar.c_str()));
207 for(vector<pair<pair<string,string>, pair<
string, vector<string> > > >::const_iterator
ev = expressionVars.begin(), eve = expressionVars.end();
ev != eve; ++
ev){
209 for (vector<string>::const_iterator it =
ev->second.second.begin(), ed =
ev->second.second.end(); it != ed; ++it) {
210 args.add(dataVars[it->c_str()]);
212 RooFormulaVar expr(
ev->first.first.c_str(),
ev->first.second.c_str(),
ev->second.first.c_str(),
args);
213 RooRealVar *
col = (RooRealVar *)
data.addColumn(expr);
214 dataVars.addClone(*col);
218 for(vector<pair<pair<string,string>, pair<string, double> > >::const_iterator tc = thresholdCategories.begin(), tce = thresholdCategories.end(); tc != tce; ++tc){
219 RooThresholdCategory
tmp(tc->first.first.c_str(), tc->first.second.c_str(), (RooAbsReal &)dataVars[tc->second.first.c_str()],
"above", 1);
220 tmp.addThreshold(tc->second.second,
"below",0);
221 RooCategory *
cat = (RooCategory *)
data.addColumn(
tmp);
222 dataVars.addClone(*cat);
227 RooMultiCategory allCats(
"allCats",
"allCats", RooArgSet(binCategories, mappedCategories));
228 data.addColumn(allCats);
231 if (effCats.size() == 1) {
232 effName = effCats.front() +
"::" + effStates.front();
233 RooMappedCategory efficiencyCategory(
"_efficiencyCategory_",
"_efficiencyCategory_", (RooCategory&)dataVars[effCats.front().c_str()],
"Failed");
234 efficiencyCategory.map(effStates.front().c_str(),
"Passed");
235 data.addColumn( efficiencyCategory );
237 RooArgSet rooEffCats;
239 for (
size_t i = 0;
i < effCats.size(); ++
i) {
240 if (
i) { multiState +=
";"; effName +=
" && "; }
241 rooEffCats.add((RooCategory &) dataVars[effCats[
i].c_str()]);
242 multiState += effStates[
i];
243 effName = effCats[
i] +
"::" + effStates[
i];
246 RooMultiCategory efficiencyMultiCategory(
"_efficiencyMultiCategory_",
"_efficiencyMultiCategory_", rooEffCats);
247 RooMappedCategory efficiencyCategory(
"_efficiencyCategory_",
"_efficiencyCategory_", efficiencyMultiCategory,
"Failed");
248 efficiencyCategory.map(multiState.c_str(),
"Passed");
249 data.addColumn( efficiencyCategory );
252 RooMappedCategory pdfCategory(
"_pdfCategory_",
"_pdfCategory_", allCats, (binToPDFmap.size()>0)?binToPDFmap[0].c_str():
"all");
253 for(
unsigned int i = 1;
i<binToPDFmap.size();
i+=2){
254 pdfCategory.map(binToPDFmap[
i].c_str(), binToPDFmap[
i+1].c_str());
256 data.addColumn( pdfCategory );
259 RooRealVar
efficiency(
"efficiency",
"Efficiency", 0, 1);
261 RooDataSet fitEfficiency(
"fit_eff",
"Efficiency from unbinned ML fit", RooArgSet(RooArgSet(binnedVariables, categories),
efficiency), StoreAsymError(RooArgSet(binnedVariables,
efficiency)));
263 RooDataSet cntEfficiency(
"cnt_eff",
"Efficiency from counting", RooArgSet(RooArgSet(binnedVariables, categories),
efficiency), StoreAsymError(RooArgSet(binnedVariables,
efficiency)));
266 if(!floatShapeParameters){
268 RooWorkspace*
w =
new RooWorkspace();
272 std::cout <<
"ALL dataset: calling doFitEfficiency with pdf: " << pdfCategory.getLabel() << std::endl;
273 doFitEfficiency(w, pdfCategory.getLabel(),
efficiency);
278 TIterator* it = allCats.typeIterator();
279 for(RooCatType*
t = (RooCatType*)it->Next();
t!=0;
t = (RooCatType*)it->Next() ){
281 TString catName =
t->GetName();
283 if(catName.Contains(
"NotMapped"))
continue;
285 RooAbsData* data_bin = (RooDataSet*)
data.reduce(Cut(TString::Format(
"allCats==%d",
t->getVal())));
287 const RooArgSet* row = data_bin->get();
290 TString pdfName(((RooCategory*)row->find(
"_pdfCategory_"))->getLabel());
294 TString dirName = catName;
295 dirName.ReplaceAll(
"{",
"").ReplaceAll(
"}",
"").ReplaceAll(
";",
"__");
296 if(pdfName.Length() > 0){
297 dirName.Append(
"__").Append(pdfName);
300 cout<<
"Fitting bin: "<<dirName<<endl;
302 gDirectory->mkdir(dirName)->cd();
306 RooWorkspace* w =
new RooWorkspace();
308 w->import(*data_bin);
310 data_bin = w->data(
"data");
313 if (doSaveDistributionsPlot) saveDistributionsPlot(w);
315 if(data_bin->numEntries()>0){
317 RooArgSet meanOfVariables;
318 RooLinkedListIter vit = binnedVariables.iterator();
319 for(RooRealVar*
v = (RooRealVar*)vit.Next();
v!=0;
v = (RooRealVar*)vit.Next() ){
320 meanOfVariables.addClone(*
v);
321 double mean = w->data(
"data")->mean(*
v);
322 RooBinning
binning((RooBinning&)
v->getBinning());
323 int ind =
binning.binNumber(mean);
324 RooRealVar& newVar = (RooRealVar&)meanOfVariables[
v->GetName()];
331 meanOfVariables.addClone(*data_bin->get(0),
true);
335 doFitEfficiency(w, pdfName.Data(),
efficiency);
336 fitEfficiency.add( RooArgSet(meanOfVariables,
efficiency) );
344 cntEfficiency.add( RooArgSet(meanOfVariables,
efficiency) );
353 gDirectory->cd(
"..");
357 fitEfficiency.Write();
358 gDirectory->mkdir(
"fit_eff_plots")->cd();
359 saveEfficiencyPlots(fitEfficiency, effName, binnedVariables, mappedCategories);
360 gDirectory->cd(
"..");
367 cntEfficiency.Write();
368 gDirectory->mkdir(
"cnt_eff_plots")->cd();
369 saveEfficiencyPlots(cntEfficiency, effName, binnedVariables, mappedCategories);
370 gDirectory->cd(
"..");
377 if(pdfName ==
"all"){
381 createPdf(w, pdfs[pdfName]);
384 std::auto_ptr<RooFitResult> res(0);
386 RooAbsData *
data = w->data(
"data");
387 std::auto_ptr<RooDataHist> bdata;
390 const RooArgSet *dataObs = data->get(0);
392 RooArgSet *obs = w->pdf(
"simPdf")->getObservables(dataObs);
393 bdata.reset(
new RooDataHist(
"data_binned",
"data_binned", *obs, *data));
395 data = w->data(
"data_binned");
399 double totPassing = data->sumEntries(
"_efficiencyCategory_==_efficiencyCategory_::Passed");
400 double totFailing = data->sumEntries(
"_efficiencyCategory_==_efficiencyCategory_::Failed");
402 RooAbsReal* simNLL = w->pdf(
"simPdf")->createNLL(*data,Extended(
true),NumCPU(numCPU));
404 RooMinimizer minimizer(*simNLL);
405 RooMinuit minuit(*simNLL);
406 minuit.setStrategy(1);
407 minuit.setProfile(
true);
408 RooProfileLL profileLL(
"simPdfNLL",
"",*simNLL,*w->var(
"efficiency"));
427 if(!fixVars.empty()){
429 if(!floatShapeParameters && fixVarValues.empty()){
440 minimizer.minimize(
"Minuit2",
"Scan");
449 minimizer.minimize(
"Minuit2",
"Scan");
461 if(!floatShapeParameters) varRestorer(w);
464 minimizer.minimize(
"Minuit2",
"Scan");
469 RooMinimizer* profMinuit = profileLL.minimizer();
470 profMinuit->setProfile(
true);
471 profMinuit->setStrategy(2);
472 profMinuit->setPrintLevel(1);
473 profMinuit->minos(*w->var(
"efficiency"));
474 res.reset( profMinuit->save() );
481 if(floatShapeParameters){
486 minimizer.minimize(
"Minuit2",
"Scan");
489 res.reset( w->pdf(
"simPdf")->fitTo(*data, Save(
true), Extended(
true), NumCPU(numCPU), Strategy(2),
490 Minos(*w->var(
"efficiency")), PrintLevel(quiet?-1:1),
491 PrintEvalErrors(quiet?-1:1), Warnings(!quiet)) );
497 res->Write(
"fitresults");
498 w->saveSnapshot(
"finalState",w->components());
501 RooRealVar*
e = (RooRealVar*) res->floatParsFinal().find(
"efficiency");
504 efficiency.setVal(e->getVal());
505 Double_t errLo = e->getErrorLo(), errHi = e->getErrorHi();
506 if (errLo == 0 && e->getVal() < 0.5) errLo = e->getMin()-e->getVal();
507 if (errHi == 0 && e->getVal() > 0.5) errHi = e->getMax()-e->getVal();
508 efficiency.setAsymError(errLo, errHi);
510 if (totPassing * totFailing == 0) {
511 RooRealVar* nTot = (RooRealVar*) res->floatParsFinal().find(
"numTot");
512 RooRealVar* fSig = (RooRealVar*) res->floatParsFinal().find(
"fSigAll");
513 double cerr = ROOT::Math::beta_quantile( 1-(1.0-.68540158589942957)/2, 1, nTot->getVal() * fSig->getVal() );
523 if (totPassing == 0) {
524 efficiency.setVal(0);
525 efficiency.setAsymError(0,cerr);
527 efficiency.setVal(1);
528 efficiency.setAsymError(-cerr,0);
537 for(
unsigned int i=0;
i<pdfCommands.size();
i++){
539 if (command.find(
"#import ") == 0) {
540 TDirectory *here = gDirectory;
541 w->import(command.substr(8).c_str());
544 TDirectory *here = gDirectory;
545 w->factory(command.c_str());
551 w->factory(
"expr::nSignalPass('efficiency*fSigAll*numTot', efficiency, fSigAll[.9,0,1],numTot[1,0,1e10])");
552 w->factory(
"expr::nSignalFail('(1-efficiency)*fSigAll*numTot', efficiency, fSigAll,numTot)");
553 w->factory(
"expr::nBkgPass('effBkg*(1-fSigAll)*numTot', effBkg[.9,0,1],fSigAll,numTot)");
554 w->factory(
"expr::nBkgFail('(1-effBkg)*(1-fSigAll)*numTot', effBkg,fSigAll,numTot)");
555 TString sPass =
"signal", sFail =
"signal";
556 if (w->pdf(
"signalPass") != 0 && w->pdf(
"signalFail") != 0) {
557 if (w->pdf(
"signal") != 0)
throw std::logic_error(
"You must either define one 'signal' PDF or two PDFs ('signalPass', 'signalFail'), not both!");
558 sPass =
"signalPass"; sFail =
"signalFail";
559 }
else if (w->pdf(
"signal") != 0) {
560 if (w->pdf(
"signalPass") != 0 || w->pdf(
"signalFail") != 0) {
561 throw std::logic_error(
"You must either define one 'signal' PDF or two PDFs ('signalPass', 'signalFail'), not both!");
564 throw std::logic_error(
"You must either define one 'signal' PDF or two PDFs ('signalPass', 'signalFail')");
566 w->factory(
"SUM::pdfPass(nSignalPass*"+sPass+
", nBkgPass*backgroundPass)");
567 w->factory(
"SUM::pdfFail(nSignalFail*"+sFail+
", nBkgFail*backgroundFail)");
569 w->factory(
"SIMUL::simPdf(_efficiencyCategory_, Passed=pdfPass, Failed=pdfFail)");
571 if (w->pdf(
"simPdf") == 0)
throw std::runtime_error(
"Could not create simultaneous fit pdf.");
572 if(w->var(
"signalFractionInPassing") == 0)
573 w->factory(
"signalFractionInPassing[0.9]");
578 double signalEfficiency = w->var(
"efficiency")->getVal();
579 double signalFractionInPassing = w->var(
"signalFractionInPassing")->getVal();
580 double totPassing = w->data(
"data")->sumEntries(
"_efficiencyCategory_==_efficiencyCategory_::Passed");
581 double totFailinging = w->data(
"data")->sumEntries(
"_efficiencyCategory_==_efficiencyCategory_::Failed");
582 double numSignalAll = totPassing*signalFractionInPassing/signalEfficiency;
587 if(numSignalAll*(1-signalEfficiency) > totFailinging)
588 numSignalAll = totFailinging;
590 w->var(
"numTot")->setVal(totPassing+totFailinging);
591 w->var(
"numTot")->setMax(2.0*(totPassing+totFailinging)+10);
593 if (totPassing == 0) {
594 w->var(
"efficiency")->setVal(0.0);
595 w->var(
"efficiency")->setAsymError(0,1);
596 w->var(
"efficiency")->setConstant(
false);
597 }
else if (totFailinging == 0) {
598 w->var(
"efficiency")->setVal(1.0);
599 w->var(
"efficiency")->setAsymError(-1,0);
600 w->var(
"efficiency")->setConstant(
false);
602 w->var(
"efficiency")->setConstant(
false);
610 w->saveSnapshot(
"initialState",w->components());
615 RooCategory& efficiencyCategory = *w->cat(
"_efficiencyCategory_");
616 RooAbsData* dataAll = (binnedFit ? w->data(
"data_binned") : w->data(
"data") );
617 RooAbsData* dataPass = dataAll->reduce(Cut(
"_efficiencyCategory_==_efficiencyCategory_::Passed"));
618 RooAbsData* dataFail = dataAll->reduce(Cut(
"_efficiencyCategory_==_efficiencyCategory_::Failed"));
619 RooAbsPdf& pdf = *w->pdf(
"simPdf");
620 std::auto_ptr<RooArgSet> obs(pdf.getObservables(*dataAll));
621 RooRealVar* mass = 0;
622 RooLinkedListIter it = obs->iterator();
623 for(RooAbsArg*
v = (RooAbsArg*)it.Next();
v!=0;
v = (RooAbsArg*)it.Next() ){
624 if(!
v->InheritsFrom(
"RooRealVar"))
continue;
625 mass = (RooRealVar*)
v;
630 TCanvas
canvas(
"fit_canvas");
632 vector<RooPlot*> frames;
636 frames.push_back(mass->frame(Name(
"Passing"), Title(
"Passing Probes")));
637 frames.push_back(mass->frame(Name(
"Failing"), Title(
"Failing Probes")));
638 frames.push_back(mass->frame(Name(
"All"), Title(
"All Probes")));
640 frames.push_back(mass->frame(Name(
"Passing"), Title(
"Passing Probes"), Bins(massBins)));
641 frames.push_back(mass->frame(Name(
"Failing"), Title(
"Failing Probes"), Bins(massBins)));
642 frames.push_back(mass->frame(Name(
"All"), Title(
"All Probes"), Bins(massBins)));
644 dataPass->plotOn(frames[0]);
645 pdf.plotOn(frames[0], Slice(efficiencyCategory,
"Passed"), ProjWData(*dataPass), LineColor(kGreen));
646 pdf.plotOn(frames[0], Slice(efficiencyCategory,
"Passed"), ProjWData(*dataPass), LineColor(kGreen), Components(
"backgroundPass"), LineStyle(kDashed));
650 dataFail->plotOn(frames[1]);
651 pdf.plotOn(frames[1], Slice(efficiencyCategory,
"Failed"), ProjWData(*dataFail), LineColor(kRed));
652 pdf.plotOn(frames[1], Slice(efficiencyCategory,
"Failed"), ProjWData(*dataFail), LineColor(kRed), Components(
"backgroundFail"), LineStyle(kDashed));
656 dataAll->plotOn(frames.back());
657 pdf.plotOn(frames.back(), ProjWData(*dataAll), LineColor(kBlue));
658 pdf.plotOn(frames.back(), ProjWData(*dataAll), LineColor(kBlue), Components(
"backgroundPass,backgroundFail"), LineStyle(kDashed));
659 frames.back()->Draw();
662 frames.push_back(mass->frame(Name(
"Fit Results"), Title(
"Fit Results")));
663 pdf.paramOn(frames.back(), dataAll,
"", 0,
"NELU", 0.1, 0.9, 0.9);
665 frames.back()->findObject(Form(
"%s_paramBox",pdf.GetName()))->Draw();
669 for (
size_t i=0;
i<frames.size();
i++) {
678 RooAbsData* dataAll = w->data(
"data");
679 RooAbsData* dataPass = dataAll->reduce(Cut(
"_efficiencyCategory_==_efficiencyCategory_::Passed"));
680 RooAbsData* dataFail = dataAll->reduce(Cut(
"_efficiencyCategory_==_efficiencyCategory_::Failed"));
682 const RooArgSet* vars = dataAll->get();
683 vector<RooRealVar*> reals;
684 RooLinkedListIter it = vars->iterator();
685 for(RooAbsArg*
v = (RooAbsArg*)it.Next();
v!=0;
v = (RooAbsArg*)it.Next() ){
686 if(!
v->InheritsFrom(
"RooRealVar"))
continue;
687 reals.push_back((RooRealVar*)
v);
689 TCanvas
canvas(
"distributions_canvas");
690 canvas.Divide(3,reals.size());
691 vector<RooPlot*> frames;
692 for(
unsigned int i=0;
i<reals.size();
i++){
695 frames.push_back(reals[
i]->frame(Name(
"Passing"), Title(
"Passing Probes"), Bins(100)));
696 dataPass->plotOn(frames.back(), LineColor(kGreen));
697 dataPass->statOn(frames.back());
698 frames.back()->Draw();
701 frames.push_back(reals[
i]->frame(Name(
"Failing"), Title(
"Failing Probes"), Bins(100)));
702 dataFail->plotOn(frames.back(), LineColor(kRed));
703 dataFail->statOn(frames.back());
704 frames.back()->Draw();
707 frames.push_back(reals[
i]->frame(Name(
"All"), Title(
"All Probes"), Bins(100)));
708 dataAll->plotOn(frames.back(), LineColor(kBlue));
709 dataAll->statOn(frames.back());
710 frames.back()->Draw();
714 for (
size_t i=0;
i<frames.size();
i++) {
722 RooLinkedListIter v1it = binnedVariables.iterator();
723 bool isOnePoint = (eff.numEntries() == 1);
724 for(RooRealVar* v1 = (RooRealVar*)v1it.Next(); v1!=0; v1 = (RooRealVar*)v1it.Next() ){
725 RooArgSet binCategories1D;
726 if (v1->numBins() == 1 && !isOnePoint)
continue;
727 RooLinkedListIter v2it = binnedVariables.iterator();
728 for(RooRealVar* v2 = (RooRealVar*)v2it.Next(); v2!=0; v2 = (RooRealVar*)v2it.Next() ){
729 if(v2 == v1)
continue;
730 if (v2->numBins() == 1 && !isOnePoint)
continue;
731 binCategories1D.addClone( RooBinningCategory(TString(v2->GetName())+
"_bins", TString(v2->GetName())+
"_bins", *v2) );
733 RooArgSet binCategories2D;
734 RooLinkedListIter v3it = binnedVariables.iterator();
735 for(RooRealVar* v3 = (RooRealVar*)v3it.Next(); v3!=0; v3 = (RooRealVar*)v3it.Next() ){
736 if(v3 == v1 || v3 == v2)
continue;
737 binCategories2D.addClone( RooBinningCategory(TString(v3->GetName())+
"_bins", TString(v3->GetName())+
"_bins", *v3) );
739 RooMultiCategory allCats2D(
"allCats2D",
"allCats2D", RooArgSet(binCategories2D, mappedCategories));
740 if(allCats2D.numTypes()==0){
741 makeEfficiencyPlot2D(eff, *v1, *v2, TString::Format(
"%s_%s_PLOT", v1->GetName(), v2->GetName()),
"", effName);
743 RooDataSet myEff(eff);
744 myEff.addColumn(allCats2D);
745 std::auto_ptr<TIterator> catIt(allCats2D.typeIterator());
746 for(RooCatType*
t = (RooCatType*)catIt->Next();
t!=0;
t = (RooCatType*)catIt->Next() ){
747 TString catName =
t->GetName();
748 if(catName.Contains(
"NotMapped"))
continue;
749 catName.ReplaceAll(
"{",
"").ReplaceAll(
"}",
"").ReplaceAll(
";",
"_&_");
750 makeEfficiencyPlot2D(myEff, *v1, *v2, TString::Format(
"%s_%s_PLOT_%s",v1->GetName(), v2->GetName(), catName.Data()), catName, effName,
"allCats1D",
t->getVal());
754 RooMultiCategory allCats1D(
"allCats1D",
"allCats1D", RooArgSet(binCategories1D, mappedCategories));
755 if(allCats1D.numTypes()==0){
756 makeEfficiencyPlot1D(eff, *v1, TString::Format(
"%s_PLOT", v1->GetName()),
"", effName);
758 RooDataSet myEff(eff);
759 myEff.addColumn(allCats1D);
760 std::auto_ptr<TIterator> catIt(allCats1D.typeIterator());
761 for(RooCatType*
t = (RooCatType*)catIt->Next();
t!=0;
t = (RooCatType*)catIt->Next() ){
762 TString catName =
t->GetName();
763 if(catName.Contains(
"NotMapped"))
continue;
764 catName.ReplaceAll(
"{",
"").ReplaceAll(
"}",
"").ReplaceAll(
";",
"_&_");
765 makeEfficiencyPlot1D(myEff, *v1, TString::Format(
"%s_PLOT_%s", v1->GetName(), catName.Data()), catName, effName,
"allCats1D",
t->getVal());
772 TGraphAsymmErrors *
p =
new TGraphAsymmErrors();
773 const RooArgSet *entry = eff.get();
774 const RooRealVar &vi =
dynamic_cast<const RooRealVar &
>(*entry->find(v.GetName()));
775 const RooRealVar &ei =
dynamic_cast<const RooRealVar &
>(*entry->find(
"efficiency"));
776 for (
unsigned int i = 0,
n = eff.numEntries();
i <
n; ++
i) {
778 if (catName != 0 && entry->getCatIndex(catName) != catIndex)
continue;
779 int j = p->GetN(); p->Set(j+1);
780 p->SetPoint(j, vi.getVal(), ei.getVal() );
781 p->SetPointError(j, -vi.getAsymErrorLo(), vi.getAsymErrorHi(), -ei.getAsymErrorLo(), ei.getAsymErrorHi() );
784 TH1F *frame =
new TH1F(
"frame",
"Efficiency of "+effName, 1, v.getMin(), v.getMax()); frame->SetDirectory(0);
785 p->SetNameTitle(Form(
"hxy_%s", eff.GetName()),
"Efficiency of "+effName);
786 p->GetXaxis()->SetTitle(strlen(v.getUnit()) ? Form(
"%s (%s)", v.GetName(), v.getUnit()) : v.GetName());
787 p->GetYaxis()->SetTitle(
"Efficiency of "+effName);
788 frame->GetXaxis()->SetTitle(strlen(v.getUnit()) ? Form(
"%s (%s)", v.GetName(), v.getUnit()) : v.GetName());
789 frame->GetYaxis()->SetTitle(
"Efficiency of "+effName);
790 frame->GetYaxis()->SetRangeUser(0,1);
792 p->SetLineWidth(2); p->SetMarkerStyle(kFullCircle); p->SetMarkerSize(1.2);
799 void TagProbeFitter::makeEfficiencyPlot2D(RooDataSet& eff, RooRealVar& v1, RooRealVar& v2,
const TString& plotName,
const TString& plotTitle,
const TString& effName,
const char *catName,
int catIndex){
801 canvas.SetRightMargin(0.15);
802 TH2F*
h =
new TH2F(plotName, plotName, v1.getBinning().numBins(), v1.getBinning().array(), v2.getBinning().numBins(), v2.getBinning().array());
803 const RooArgSet* set = eff.get();
804 RooRealVar*
e = (RooRealVar*) set->find(
"efficiency");
805 RooRealVar* v1_ = (RooRealVar*) set->find(v1.GetName());
806 RooRealVar* v2_ = (RooRealVar*) set->find(v2.GetName());
807 h->SetTitle(TString::Format(
"%s;%s%s;%s%s;Efficiency of %s", plotTitle.Data(),
808 v1.GetTitle(), TString(v1.getUnit()).Length()==0?
"":TString::Format(
" (%s)", v1.getUnit()).
Data(),
809 v2.GetTitle(), TString(v2.getUnit()).Length()==0?
"":TString::Format(
" (%s)", v2.getUnit()).
Data(), effName.
Data()));
810 h->SetOption(
"colztexte");
811 h->GetZaxis()->SetRangeUser(-0.001,1.001);
813 for(
int i=0;
i<eff.numEntries();
i++){
814 const RooArgSet *entry = eff.get(
i);
815 if (catName != 0 && entry->getCatIndex(catName) != catIndex)
continue;
816 h->SetBinContent(h->FindBin(v1_->getVal(), v2_->getVal()), e->getVal());
817 h->SetBinError(h->FindBin(v1_->getVal(), v2_->getVal()), (e->getErrorHi()-e->getErrorLo())/2.);
829 int pass = w->data(
"data")->sumEntries(
"_efficiencyCategory_==_efficiencyCategory_::Passed");
830 int fail = w->data(
"data")->sumEntries(
"_efficiencyCategory_==_efficiencyCategory_::Failed");
831 double e = (pass+fail == 0) ? 0 : pass/
double(pass+fail);
833 double alpha = (1.0 - .68540158589942957)/2;
834 double lo = (pass == 0) ? 0.0 : ROOT::Math::beta_quantile( alpha, pass, fail+1 );
835 double hi = (fail == 0) ? 1.0 : ROOT::Math::beta_quantile( 1-alpha, pass+1, fail );
840 efficiency.setVal(e);
841 efficiency.setAsymError(lo-e, hi-e);
845 std::vector<std::string>::const_iterator it;
846 for(it=fixVars.begin(); it<fixVars.end(); it++){
847 if(w->var((*it).c_str()))
848 w->var((*it).c_str())->setAttribute(
"Constant",fix);
850 std::cout <<
"TagProbeFitter: " <<
"Can't find a variable to fix: " << *it;
856 if(!fixVarValues.empty()){
857 std::cout <<
"attempt to save variables more than once!" << std::endl;
860 std::vector<std::string>::const_iterator it;
861 for(it=fixVars.begin(); it<fixVars.end(); it++){
862 fixVarValues.push_back(w->var((*it).c_str())->getVal());
868 if(fixVarValues.size()==fixVars.size())
869 for(
unsigned int i=0;
i< fixVars.size();
i++){
870 std::cout <<
"setting variable " << fixVars[
i].c_str() << std::endl;
871 w->var(fixVars[
i].c_str())->setVal(fixVarValues[
i]);
874 std::cout <<
"fixVars and fixVarValues are not of the same size!" << std::endl;
void doCntEfficiency(RooWorkspace *w, RooRealVar &efficiency)
calculate the efficiecny by counting in the dataset found in the workspace
void addPdf(std::string pdfName, std::vector< std::string > &pdfCommands)
static PFTauRenderPlugin instance
void setWeightVar(const std::string &weight)
set a variable to be used as weight for a dataset. empty string means no weights. ...
std::string calculateEfficiency(std::string dirName, std::string efficiencyCategory, std::string efficiencyState, std::vector< std::string > &unbinnedVariables, std::map< std::string, std::vector< double > > &binnedReals, std::map< std::string, std::vector< std::string > > &binnedCategories, std::vector< std::string > &binToPDFmap)
calculate the efficiency for a particular binning of the data; it saves everything in the directory "...
bool addThresholdCategory(std::string categoryName, std::string title, std::string varName, double cutValue)
adds a new category based on a cut
void setQuiet(bool quiet_=true)
suppress most of the output from RooFit and Minuit
TagProbeFitter(const std::vector< std::string > &inputFileNames, std::string inputDirectoryName, std::string inputTreeName, std::string outputFileName, int numCPU=1, bool saveWorkspace_=false, bool floatShapeParameters=true, const std::vector< std::string > &fixVars_=std::vector< std::string >())
construct the fitter with the inputFileName, inputDirectoryName, inputTreeName, outputFileName and sp...
bool addExpression(std::string expressionName, std::string title, std::string expression, const std::vector< std::string > &arguments)
adds a new category based on a cut
void saveEfficiencyPlots(RooDataSet &eff, const TString &effName, RooArgSet &binnedVariables, RooArgSet &mappedCategories)
saves the efficiency plots
void varSaver(RooWorkspace *w)
store values in the vector
void makeEfficiencyPlot2D(RooDataSet &eff, RooRealVar &v1, RooRealVar &v2, const TString &plotName, const TString &plotTitle, const TString &effName, const char *catName=0, int catIndex=-1)
makes the 2D plot
void varRestorer(RooWorkspace *w)
restore variables's values for fit starting point
void saveFitPlot(RooWorkspace *w)
saves the fit canvas
~TagProbeFitter()
destructor closes the files
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
void setInitialValues(RooWorkspace *w)
sets initial values of the PDF parameters based on the data available in the workspace ...
void saveDistributionsPlot(RooWorkspace *w)
saves the distributions canvas
void varFixer(RooWorkspace *w, bool fix)
fix or release variables selected by user
std::vector< std::vector< double > > tmp
TString units(TString variable, Char_t axis)
char data[epos_bytes_allocation]
void setBinsForMassPlots(int bins)
set number of bins to use when making the plots; 0 = automatic
void doSBSEfficiency(RooWorkspace *w, RooRealVar &efficiency)
calculate the efficiecny with side band substraction in the dataset found in the workspace ...
void createPdf(RooWorkspace *w, std::vector< std::string > &pdfCommands)
creates the simultaneous PDF in the workspace according to the "pdfCommands"
void doFitEfficiency(RooWorkspace *w, std::string pdfName, RooRealVar &efficiency)
calculate the efficiecny with a simulataneous maximum likelihood fit in the dataset found in the work...
bool addCategory(std::string categoryName, std::string title, std::string expression)
adds a new category variable to the set of variables describing the data in the tree; "expression" is...
bool addVariable(std::string variableName, std::string title, double low, double hi, std::string units)
adds a new real variable to the set of variables describing the data in the tree
void makeEfficiencyPlot1D(RooDataSet &eff, RooRealVar &v, const TString &plotName, const TString &plotTitle, const TString &effName, const char *catName=0, int catIndex=-1)
makes the 1D plot