Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-01-27 02:45:32

0001 #include "PhysicsTools/TagAndProbe/interface/TagProbeFitter.h"
0002 #include <memory>
0003 
0004 #include <stdexcept>
0005 
0006 #include "TROOT.h"
0007 #include "TFile.h"
0008 #include "TPad.h"
0009 #include "TText.h"
0010 #include "TCanvas.h"
0011 #include "TGraphAsymmErrors.h"
0012 #include "TH2F.h"
0013 #include "TStyle.h"
0014 #include "Math/QuantFuncMathCore.h"
0015 #include "RooAbsReal.h"
0016 #include "RooBinning.h"
0017 #include "RooBinningCategory.h"
0018 #include "RooCategory.h"
0019 #include "RooDataHist.h"
0020 #include "RooDataSet.h"
0021 #include "RooFitLegacy/RooCatTypeLegacy.h"
0022 #include "RooFitResult.h"
0023 #include "RooFormulaVar.h"
0024 #include "RooGlobalFunc.h"
0025 #include "RooLinkedListIter.h"
0026 #include "RooMappedCategory.h"
0027 #include "RooMinimizer.h"
0028 #include "RooMsgService.h"
0029 #include "RooMultiCategory.h"
0030 #include "RooNumIntConfig.h"
0031 #include "RooPlot.h"
0032 #include "RooProdPdf.h"
0033 #include "RooProfileLL.h"
0034 #include "RooRealVar.h"
0035 #include "RooThresholdCategory.h"
0036 #include "RooTrace.h"
0037 #include "RooWorkspace.h"
0038 #include "RooTreeDataStore.h"
0039 
0040 using namespace std;
0041 using namespace RooFit;
0042 
0043 TagProbeFitter::TagProbeFitter(const std::vector<std::string>& inputFileNames,
0044                                string inputDirectoryName,
0045                                string inputTreeName,
0046                                string outputFileName,
0047                                int numCPU_,
0048                                bool saveWorkspace_,
0049                                bool floatShapeParameters_,
0050                                const std::vector<std::string>& fixVars_) {
0051   inputTree = new TChain((inputDirectoryName + "/" + inputTreeName).c_str());
0052   for (size_t i = 0; i < inputFileNames.size(); i++) {
0053     inputTree->Add(inputFileNames[i].c_str());
0054   }
0055   outputFile = new TFile(outputFileName.c_str(), "recreate");
0056   outputDirectory = outputFile->mkdir(inputDirectoryName.c_str());
0057   numCPU = numCPU_;
0058   saveWorkspace = saveWorkspace_;
0059   massBins = 0;  // automatic default
0060   floatShapeParameters = floatShapeParameters_;
0061   fixVars = fixVars_;
0062   weightVar = "";
0063   if (!floatShapeParameters && fixVars.empty())
0064     std::cout << "TagProbeFitter: "
0065               << "You wnat to fix some variables but do not specify them!";
0066 
0067   gROOT->SetStyle("Plain");
0068   gStyle->SetTitleFillColor(0);
0069   gStyle->SetPalette(1);
0070   gStyle->SetOptStat(0);
0071   gStyle->SetPaintTextFormat(".2f");
0072 
0073   quiet = false;
0074 
0075   binnedFit = false;
0076 
0077   doSaveDistributionsPlot = true;
0078 
0079   // make integration very precise
0080   RooAbsReal::defaultIntegratorConfig()->setEpsAbs(1e-13);
0081   RooAbsReal::defaultIntegratorConfig()->setEpsRel(1e-13);
0082 
0083   split_mode = 0;
0084 }
0085 
0086 TagProbeFitter::~TagProbeFitter() {
0087   if (inputTree)
0088     delete inputTree;
0089   if (outputFile)
0090     outputFile->Close();
0091 }
0092 
0093 void TagProbeFitter::setQuiet(bool quiet_) {
0094   quiet = quiet_;
0095   if (quiet) {
0096     RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
0097   } else {
0098     RooMsgService::instance().setGlobalKillBelow(RooFit::WARNING);
0099   }
0100 }
0101 
0102 void TagProbeFitter::setSplitMode(unsigned int nevents) { split_mode = nevents; }
0103 
0104 bool TagProbeFitter::addVariable(string name, string title, double low, double hi, string units) {
0105   RooRealVar temp(name.c_str(), title.c_str(), low, hi, units.c_str());
0106   temp.setBins(5000, "cache");
0107   variables.addClone(temp);
0108   return true;
0109 }
0110 
0111 bool TagProbeFitter::addCategory(string name, string title, string expression) {
0112   RooCategory* c = (RooCategory*)parameterParser.factory(expression.c_str());
0113   if (!c)
0114     return false;
0115   //set the name of the category to the passed name instead of the one in the expression
0116   c->SetName(name.c_str());
0117   c->SetTitle(title.c_str());
0118   variables.addClone(*c);
0119   return true;
0120 }
0121 
0122 bool TagProbeFitter::addExpression(string expressionName,
0123                                    string title,
0124                                    string expression,
0125                                    const std::vector<string>& arguments) {
0126   expressionVars.push_back(make_pair(make_pair(expressionName, title), make_pair(expression, arguments)));
0127   return true;
0128 }
0129 
0130 bool TagProbeFitter::addThresholdCategory(string categoryName, string title, string varName, double cutValue) {
0131   thresholdCategories.push_back(make_pair(make_pair(categoryName, title), make_pair(varName, cutValue)));
0132   return true;
0133 }
0134 
0135 void TagProbeFitter::addPdf(string name, vector<string>& pdfCommands) { pdfs[name] = pdfCommands; }
0136 
0137 void TagProbeFitter::setBinsForMassPlots(int bins) { massBins = bins; }
0138 
0139 void TagProbeFitter::setWeightVar(const std::string& var) { weightVar = var; }
0140 
0141 string TagProbeFitter::calculateEfficiency(string dirName,
0142                                            const std::vector<string>& effCats,
0143                                            const std::vector<string>& effStates,
0144                                            vector<string>& unbinnedVariables,
0145                                            map<string, vector<double> >& binnedReals,
0146                                            map<string, std::vector<string> >& binnedCategories,
0147                                            vector<string>& binToPDFmap) {
0148   //go to home directory
0149   outputDirectory->cd();
0150   //make a directory corresponding to this efficiency binning
0151   gDirectory->mkdir(dirName.c_str())->cd();
0152 
0153   RooArgSet dataVars;
0154 
0155   //collect unbinned variables
0156   for (vector<string>::iterator v = unbinnedVariables.begin(); v != unbinnedVariables.end(); v++) {
0157     dataVars.addClone(variables[v->c_str()], true);
0158     if (binnedFit && (v == unbinnedVariables.begin())) {
0159       ((RooRealVar&)dataVars[v->c_str()]).setBins(massBins);
0160     }
0161   }
0162   //collect the binned variables and the corresponding bin categories
0163   RooArgSet binnedVariables;
0164   RooArgSet binCategories;
0165   for (map<string, vector<double> >::iterator v = binnedReals.begin(); v != binnedReals.end(); v++) {
0166     TString name = v->first;
0167     if (variables.find(name) == nullptr) {
0168       cerr << "Binned variable '" << name << "' not found." << endl;
0169       return "Error";
0170     }
0171     binnedVariables.add(*dataVars.addClone(variables[name]));
0172     ((RooRealVar&)binnedVariables[name]).setBinning(RooBinning(v->second.size() - 1, &v->second[0]));
0173     binCategories.addClone(RooBinningCategory(name + "_bins", name + "_bins", (RooRealVar&)binnedVariables[name]));
0174   }
0175 
0176   //collect the category variables and the corresponding mapped categories
0177   RooArgSet categories;
0178   RooArgSet mappedCategories;
0179   for (map<string, vector<string> >::iterator v = binnedCategories.begin(); v != binnedCategories.end(); v++) {
0180     TString name = v->first;
0181     if (variables.find(name) == nullptr) {
0182       cerr << "Binned category '" << name << "' not found." << endl;
0183       return "Error";
0184     }
0185     categories.add(*dataVars.addClone(variables[name]));
0186     mappedCategories.addClone(RooMappedCategory(name + "_bins", name + "_bins", (RooCategory&)categories[name]));
0187     for (unsigned int i = 0; i < v->second.size(); i++) {
0188       ((RooMappedCategory&)mappedCategories[name + "_bins"])
0189           .map(v->second[i].c_str(), name + "_" + TString(v->second[i].c_str()).ReplaceAll(",", "_"));
0190     }
0191   }
0192 
0193   // add the efficiency category if it's not a dynamic one
0194   for (vector<string>::const_iterator effCat = effCats.begin(); effCat != effCats.end(); ++effCat) {
0195     if (variables.find(effCat->c_str()) != nullptr) {
0196       dataVars.addClone(variables[effCat->c_str()], true);
0197     }
0198   }
0199 
0200   //  add all variables used in expressions
0201   for (vector<pair<pair<string, string>, pair<string, vector<string> > > >::const_iterator ev = expressionVars.begin(),
0202                                                                                            eve = expressionVars.end();
0203        ev != eve;
0204        ++ev) {
0205     for (vector<string>::const_iterator it = ev->second.second.begin(), ed = ev->second.second.end(); it != ed; ++it) {
0206       // provided that they are real variables themselves
0207       if (variables.find(it->c_str()))
0208         dataVars.addClone(variables[it->c_str()], true);
0209     }
0210   }
0211   // add all real variables used in cuts
0212   for (vector<pair<pair<string, string>, pair<string, double> > >::const_iterator tc = thresholdCategories.begin(),
0213                                                                                   tce = thresholdCategories.end();
0214        tc != tce;
0215        ++tc) {
0216     if (variables.find(tc->second.first.c_str()))
0217       dataVars.addClone(variables[tc->second.first.c_str()], true);
0218   }
0219 
0220   //now add the necessary mass and passing variables to make the unbinned RooDataSet
0221   RooDataSet* data(nullptr);
0222   if (not split_mode) {
0223     data = new RooDataSet("data",
0224                           "data",
0225                           dataVars,
0226                           Import(*inputTree),
0227                           /*selExpr=*/Cut(""),
0228                           /*wgtVarName=*/WeightVar(weightVar.empty() ? nullptr : weightVar.c_str()));
0229 
0230     // Now add all expressions that are computed dynamically
0231     for (vector<pair<pair<string, string>, pair<string, vector<string> > > >::const_iterator
0232              ev = expressionVars.begin(),
0233              eve = expressionVars.end();
0234          ev != eve;
0235          ++ev) {
0236       RooArgList args;
0237       for (vector<string>::const_iterator it = ev->second.second.begin(), ed = ev->second.second.end(); it != ed;
0238            ++it) {
0239         args.add(dataVars[it->c_str()]);
0240       }
0241       RooFormulaVar expr(ev->first.first.c_str(), ev->first.second.c_str(), ev->second.first.c_str(), args);
0242       RooRealVar* col = (RooRealVar*)data->addColumn(expr);
0243       dataVars.addClone(*col);
0244     }
0245 
0246     // And add all dynamic categories from thresholds
0247     for (vector<pair<pair<string, string>, pair<string, double> > >::const_iterator tc = thresholdCategories.begin(),
0248                                                                                     tce = thresholdCategories.end();
0249          tc != tce;
0250          ++tc) {
0251       RooThresholdCategory tmp(tc->first.first.c_str(),
0252                                tc->first.second.c_str(),
0253                                (RooAbsReal&)dataVars[tc->second.first.c_str()],
0254                                "above",
0255                                1);
0256       tmp.addThreshold(tc->second.second, "below", 0);
0257       RooCategory* cat = (RooCategory*)data->addColumn(tmp);
0258       dataVars.addClone(*cat);
0259     }
0260   }
0261 
0262   //merge the bin categories to a MultiCategory for convenience
0263   RooMultiCategory allCats("allCats", "allCats", RooArgSet(binCategories, mappedCategories));
0264   string effName;
0265 
0266   if (not split_mode) {
0267     data->addColumn(allCats);
0268     //setup the efficiency category
0269     if (effCats.size() == 1) {
0270       effName = effCats.front() + "::" + effStates.front();
0271       RooMappedCategory efficiencyCategory(
0272           "_efficiencyCategory_", "_efficiencyCategory_", (RooCategory&)dataVars[effCats.front().c_str()], "Failed");
0273       efficiencyCategory.map(effStates.front().c_str(), "Passed");
0274       data->addColumn(efficiencyCategory);
0275     } else {
0276       RooArgSet rooEffCats;
0277       string multiState = "{";
0278       for (size_t i = 0; i < effCats.size(); ++i) {
0279         if (i) {
0280           multiState += ";";
0281           effName += " && ";
0282         }
0283         rooEffCats.add((RooCategory&)dataVars[effCats[i].c_str()]);
0284         multiState += effStates[i];
0285         effName = effCats[i] + "::" + effStates[i];
0286       }
0287       multiState += "}";
0288       RooMultiCategory efficiencyMultiCategory("_efficiencyMultiCategory_", "_efficiencyMultiCategory_", rooEffCats);
0289       RooMappedCategory efficiencyCategory(
0290           "_efficiencyCategory_", "_efficiencyCategory_", efficiencyMultiCategory, "Failed");
0291       efficiencyCategory.map(multiState.c_str(), "Passed");
0292       data->addColumn(efficiencyCategory);
0293     }
0294   }
0295 
0296   //setup the pdf category
0297   RooMappedCategory pdfCategory(
0298       "_pdfCategory_", "_pdfCategory_", allCats, (!binToPDFmap.empty()) ? binToPDFmap[0].c_str() : "all");
0299   for (unsigned int i = 1; i < binToPDFmap.size(); i += 2) {
0300     pdfCategory.map(binToPDFmap[i].c_str(), binToPDFmap[i + 1].c_str());
0301   }
0302   if (not split_mode)
0303     data->addColumn(pdfCategory);
0304 
0305   //create the empty efficiency datasets from the binned variables
0306   RooRealVar efficiency("efficiency", "Efficiency", 0, 1);
0307 
0308   RooDataSet fitEfficiency("fit_eff",
0309                            "Efficiency from unbinned ML fit",
0310                            RooArgSet(RooArgSet(binnedVariables, categories), efficiency),
0311                            StoreAsymError(RooArgSet(binnedVariables, efficiency)));
0312 
0313   RooDataSet cntEfficiency("cnt_eff",
0314                            "Efficiency from counting",
0315                            RooArgSet(RooArgSet(binnedVariables, categories), efficiency),
0316                            StoreAsymError(RooArgSet(binnedVariables, efficiency)));
0317 
0318   if (not split_mode) {
0319     if (!floatShapeParameters) {
0320       //fitting whole dataset to get initial values for some parameters
0321       RooWorkspace* w = new RooWorkspace();
0322       w->import(*data);
0323       efficiency.setVal(0);  //reset
0324       efficiency.setAsymError(0, 0);
0325       std::cout << "ALL dataset: calling doFitEfficiency with pdf: " << pdfCategory.getLabel() << std::endl;
0326       doFitEfficiency(w, pdfCategory.getLabel(), efficiency);
0327       delete w;
0328     }
0329   } else {
0330     // disactive not needed branches
0331     inputTree->SetBranchStatus("*", false);
0332     for (TObject* obj : dataVars)
0333       inputTree->SetBranchStatus(obj->GetName(), true);
0334   }
0335 
0336   // loop over all bins with the help of allCats
0337   // store index values in a separate vector to avoid issues
0338   // with iteration over changing allCats object
0339   std::vector<Int_t> allCatValues;
0340   TIterator* it = allCats.typeIterator();
0341   for (RooCatType* t = (RooCatType*)it->Next(); t != nullptr; t = (RooCatType*)it->Next())
0342     allCatValues.push_back(t->getVal());
0343   for (auto iCat : allCatValues) {
0344     const RooCatType* t = allCats.lookupType(iCat);
0345     //name of the multicategory
0346     TString catName = t->GetName();
0347     //skip unmapped states
0348     if (catName.Contains("NotMapped"))
0349       continue;
0350 
0351     RooDataSet* data_bin(nullptr);
0352     RooArgSet tmpVars;
0353 
0354     if (not split_mode) {
0355       //create the dataset
0356       data_bin = (RooDataSet*)data->reduce(Cut(TString::Format("allCats==%d", iCat)));
0357     } else {
0358       data_bin = new RooDataSet("data", "data", dataVars, WeightVar(weightVar.empty() ? nullptr : weightVar.c_str()));
0359 
0360       TDirectory* tmp = gDirectory;
0361       gROOT->cd();
0362 
0363       // loop over input data and fill the dataset with events for
0364       // current category
0365       unsigned int n_entries = inputTree->GetEntries();
0366       printf("Input number of events: %u\n", n_entries);
0367       unsigned int first_entry = 0;
0368       while (first_entry < n_entries) {
0369         TTree* copyTree = inputTree->CopyTree("", "", split_mode, first_entry);
0370         RooTreeDataStore store("reader",
0371                                "reader",
0372                                dataVars,
0373                                *copyTree,
0374                                /*selExpr=*/"",
0375                                /*wgtVarName=*/(weightVar.empty() ? nullptr : weightVar.c_str()));
0376         for (unsigned int i = 0; i < store.GetEntries(); ++i) {
0377           store.get(i);
0378           if (allCats.getIndex() == iCat) {
0379             data_bin->add(dataVars, weightVar.empty() ? 1.0 : dataVars.getRealValue(weightVar.c_str()));
0380           }
0381         }
0382         delete copyTree;
0383         first_entry += split_mode;
0384         data_bin->Print("V");
0385       }
0386       data_bin->Print("V");
0387       tmp->cd();
0388 
0389       // Now add all expressions that are computed dynamically
0390       for (vector<pair<pair<string, string>, pair<string, vector<string> > > >::const_iterator
0391                ev = expressionVars.begin(),
0392                eve = expressionVars.end();
0393            ev != eve;
0394            ++ev) {
0395         RooArgList args;
0396         for (vector<string>::const_iterator it = ev->second.second.begin(), ed = ev->second.second.end(); it != ed;
0397              ++it) {
0398           args.add(dataVars[it->c_str()]);
0399         }
0400         RooFormulaVar expr(ev->first.first.c_str(), ev->first.second.c_str(), ev->second.first.c_str(), args);
0401         RooRealVar* col = (RooRealVar*)data_bin->addColumn(expr);
0402         tmpVars.add(*dataVars.addClone(*col));
0403       }
0404 
0405       // And add all dynamic categories from thresholds
0406       for (vector<pair<pair<string, string>, pair<string, double> > >::const_iterator tc = thresholdCategories.begin(),
0407                                                                                       tce = thresholdCategories.end();
0408            tc != tce;
0409            ++tc) {
0410         RooThresholdCategory tmp(tc->first.first.c_str(),
0411                                  tc->first.second.c_str(),
0412                                  (RooAbsReal&)dataVars[tc->second.first.c_str()],
0413                                  "above",
0414                                  1);
0415         tmp.addThreshold(tc->second.second, "below", 0);
0416         RooCategory* cat = (RooCategory*)data_bin->addColumn(tmp);
0417         tmpVars.add(*dataVars.addClone(*cat));
0418       }
0419 
0420       //setup the efficiency category
0421       if (effCats.size() == 1) {
0422         effName = effCats.front() + "::" + effStates.front();
0423         RooMappedCategory efficiencyCategory(
0424             "_efficiencyCategory_", "_efficiencyCategory_", (RooCategory&)dataVars[effCats.front().c_str()], "Failed");
0425         efficiencyCategory.map(effStates.front().c_str(), "Passed");
0426         data_bin->addColumn(efficiencyCategory);
0427       } else {
0428         RooArgSet rooEffCats;
0429         string multiState = "{";
0430         for (size_t i = 0; i < effCats.size(); ++i) {
0431           if (i) {
0432             multiState += ";";
0433             effName += " && ";
0434           }
0435           rooEffCats.add((RooCategory&)dataVars[effCats[i].c_str()]);
0436           multiState += effStates[i];
0437           effName = effCats[i] + "::" + effStates[i];
0438         }
0439         multiState += "}";
0440         RooMultiCategory efficiencyMultiCategory("_efficiencyMultiCategory_", "_efficiencyMultiCategory_", rooEffCats);
0441         RooMappedCategory efficiencyCategory(
0442             "_efficiencyCategory_", "_efficiencyCategory_", efficiencyMultiCategory, "Failed");
0443         efficiencyCategory.map(multiState.c_str(), "Passed");
0444         data_bin->addColumn(efficiencyCategory);
0445       }
0446       data_bin->addColumn(pdfCategory);
0447     }
0448 
0449     //set the category variables by reading the first event
0450     const RooArgSet* row = data_bin->get();
0451 
0452     //get PDF name
0453     TString pdfName(((RooCategory*)row->find("_pdfCategory_"))->getLabel());
0454 
0455     //make directory name
0456     TString dirName = catName;
0457     dirName.ReplaceAll("{", "").ReplaceAll("}", "").ReplaceAll(";", "__");
0458     if (pdfName.Length() > 0) {
0459       dirName.Append("__").Append(pdfName);
0460     }
0461 
0462     cout << "Fitting bin:  " << dirName << endl;
0463     //make a directory for each bin
0464     gDirectory->mkdir(dirName)->cd();
0465 
0466     //create a workspace
0467     RooWorkspace* w = new RooWorkspace();
0468     //import the data
0469     w->import(*data_bin);
0470     delete data_bin;  // clean up earlier
0471     data_bin =
0472         dynamic_cast<RooDataSet*>(w->data("data"));  // point to the data that's in the workspace now (saves memory)
0473 
0474     //save the distribution of variables
0475     if (doSaveDistributionsPlot)
0476       saveDistributionsPlot(w);
0477     //do the fitting only if there is sufficient number of events
0478     if (data_bin->numEntries() > 0) {
0479       //set the values of binnedVariables to the mean value in this data bin
0480       RooArgSet meanOfVariables;
0481       for (const RooAbsArg* vv : binnedVariables) {
0482         const RooRealVar* v = dynamic_cast<const RooRealVar*>(vv);
0483         meanOfVariables.addClone(*v);
0484         double mean = w->data("data")->mean(*v);
0485         RooBinning binning((RooBinning&)v->getBinning());
0486         int ind = binning.binNumber(mean);
0487         RooRealVar& newVar = (RooRealVar&)meanOfVariables[v->GetName()];
0488         newVar.setVal(mean);
0489         newVar.setAsymError(binning.binLow(ind) - mean, binning.binHigh(ind) - mean);
0490       }
0491 
0492       //put an entry in the efficiency dataset
0493       //note that the category values are coming from data_bin->get(0)
0494       meanOfVariables.addClone(*data_bin->get(0), true);
0495 
0496       efficiency.setVal(0);  //reset
0497       efficiency.setAsymError(0, 0);
0498       doFitEfficiency(w, pdfName.Data(), efficiency);
0499       fitEfficiency.add(RooArgSet(meanOfVariables, efficiency));
0500 
0501       /*      efficiency.setVal(0);//reset
0502       doSBSEfficiency(w, efficiency);
0503       sbsEfficiency.add( RooArgSet(meanOfVariables, efficiency) );*/
0504 
0505       efficiency.setVal(0);  //reset
0506       doCntEfficiency(w, efficiency);
0507       cntEfficiency.add(RooArgSet(meanOfVariables, efficiency));
0508     }
0509     //save the workspace if requested
0510     if (saveWorkspace) {
0511       w->Write("w");
0512     }
0513     //clean up
0514     delete w;
0515     if (split_mode)
0516       dataVars.remove(tmpVars);
0517     //get back to the initial directory
0518     gDirectory->cd("..");
0519   }
0520 
0521   //save the efficiency data
0522   fitEfficiency.Write();
0523   gDirectory->mkdir("fit_eff_plots")->cd();
0524   saveEfficiencyPlots(fitEfficiency, effName, binnedVariables, mappedCategories);
0525   gDirectory->cd("..");
0526 
0527   /*  sbsEfficiency.Write();
0528   gDirectory->mkdir("sbs_eff_plots")->cd();
0529   saveEfficiencyPlots(sbsEfficiency, effCat+"::"+effState, binnedVariables, mappedCategories);
0530   gDirectory->cd("..");*/
0531 
0532   cntEfficiency.Write();
0533   gDirectory->mkdir("cnt_eff_plots")->cd();
0534   saveEfficiencyPlots(cntEfficiency, effName, binnedVariables, mappedCategories);
0535   gDirectory->cd("..");
0536 
0537   if (not split_mode)
0538     delete data;
0539 
0540   //empty string means no error
0541   return "";
0542 }
0543 
0544 void TagProbeFitter::doFitEfficiency(RooWorkspace* w, string pdfName, RooRealVar& efficiency) {
0545   //if pdfName is empty skip the fit
0546   if (pdfName == "all") {
0547     return;
0548   }
0549   //create the simultaneous pdf of name pdfName
0550   createPdf(w, pdfs[pdfName]);
0551   //set the initial values for the yields of signal and background
0552   setInitialValues(w);
0553   std::unique_ptr<RooFitResult> res;
0554 
0555   RooAbsData* data = w->data("data");
0556   std::unique_ptr<RooDataHist> bdata;
0557   if (binnedFit) {
0558     // get variables from data, which contain also other binning or expression variables
0559     const RooArgSet* dataObs = data->get(0);
0560     // remove everything which is not a dependency of the pdf
0561     RooArgSet* obs = w->pdf("simPdf")->getObservables(dataObs);
0562     bdata = std::make_unique<RooDataHist>("data_binned", "data_binned", *obs, *data);
0563     w->import(*bdata);
0564     data = w->data("data_binned");
0565     delete obs;
0566   }
0567 
0568   double totPassing = data->sumEntries("_efficiencyCategory_==_efficiencyCategory_::Passed");
0569   double totFailing = data->sumEntries("_efficiencyCategory_==_efficiencyCategory_::Failed");
0570 
0571   RooAbsReal* simNLL = w->pdf("simPdf")->createNLL(*data, Extended(true), NumCPU(numCPU));
0572 
0573   RooMinimizer minimizer(*simNLL);  // we are going to use this for 'scan'
0574   minimizer.setStrategy(1);
0575   minimizer.setProfile(true);
0576   RooProfileLL profileLL("simPdfNLL", "", *simNLL, *w->var("efficiency"));
0577 
0578   //******* The block of code below is to make the fit converge faster.
0579   // ****** This part is OPTIONAL, i.e., off be default. User can activate this
0580   // ****** by setting the following parameters: "fixVars" and "floatShapeParameters"
0581   // ****** Here is the full logic:
0582   /////   ---> if "floatShapeParameters==true" && "fixVars is empty" :
0583   ////             Just perform the fit without any of these optimizations (this is the default logic)
0584   /////   ---> if "floatShapeParameters==true" && "fixVars is NOT empty" :
0585   ////             In each bin, fix the desired parameters and perform the fit to obtain a good starting point.
0586   /////            Then float these parameters and fit again.
0587   ////    ---> if "floatShapeParameters==false" && "fixVars is empty" :
0588   ////             Do not perform any fit at all. Just print error message.
0589   ////    ---> if "floatShapeParameters==false" && "fixVars is NOT empty" :
0590   ///              Perform a global fit to the whole sample, save the fitted values of the
0591   ///              user specified parameters, and fix them for bin-by-bin fit.
0592 
0593   if (!fixVars.empty()) {
0594     // calculate initial values for parameters user want to fix
0595     if (!floatShapeParameters && fixVarValues.empty()) {
0596       // we want to fix these parameters for each bin.
0597       // the following sequence fixes them, fits, releases and fits again
0598       // to get reasonable values.
0599       // ----------------------------------------------------------------------
0600       // This procedure works only once with a whole dataset (without binning)
0601       // ----------------------------------------------------------------------
0602 
0603       // fix them
0604       varFixer(w, true);
0605       //do fit
0606       minimizer.minimize("Minuit2", "Scan");
0607       minimizer.migrad();
0608       minimizer.hesse();
0609       //minuit.minos();
0610       //w->pdf("simPdf")->fitTo(*data, Save(true), Extended(true), NumCPU(numCPU), Strategy(2),
0611       //PrintLevel(quiet?-1:1), PrintEvalErrors(quiet?-1:1), Warnings(!quiet));
0612       //release vars
0613       varFixer(w, false);
0614       //do fit
0615       minimizer.minimize("Minuit2", "Scan");
0616       minimizer.migrad();
0617       minimizer.hesse();
0618       //minuit.minos();
0619       //w->pdf("simPdf")->fitTo(*data, Save(true), Extended(true), NumCPU(numCPU), Strategy(2),
0620       //PrintLevel(quiet?-1:1), PrintEvalErrors(quiet?-1:1), Warnings(!quiet));
0621       //save vars
0622       varSaver(w);
0623       // now we have a starting point. Fit will converge faster.
0624     }
0625 
0626     // here we can use initial values if we want (this works for each bin)
0627     if (!floatShapeParameters)
0628       varRestorer(w);  //restore vars
0629 
0630     //do fit
0631     minimizer.minimize("Minuit2", "Scan");
0632     minimizer.migrad();
0633     minimizer.hesse();
0634     // initialize the profile likelihood
0635     profileLL.getVal();
0636     RooMinimizer* profMinuit = profileLL.minimizer();
0637     profMinuit->setProfile(true);
0638     profMinuit->setStrategy(2);
0639     profMinuit->setPrintLevel(1);
0640     profMinuit->minos(*w->var("efficiency"));
0641     res.reset(profMinuit->save());
0642     //res = w->pdf("simPdf")->fitTo(*data, Save(true), Extended(true), NumCPU(numCPU), Strategy(2),
0643     //Minos(*w->var("efficiency")), PrintLevel(quiet?-1:1),
0644     //PrintEvalErrors(quiet?-1:1), Warnings(!quiet));
0645   }  //if(!fixVars.empty())
0646 
0647   // (default = true) if we don't want to fix any parameters or want to fit each bin with all parameters floating
0648   if (floatShapeParameters) {
0649     //release vars
0650     varFixer(w, false);
0651 
0652     //do fit
0653     minimizer.minimize("Minuit2", "Scan");
0654     minimizer.migrad();
0655     minimizer.hesse();
0656     res.reset(w->pdf("simPdf")->fitTo(*data,
0657                                       Save(true),
0658                                       Extended(true),
0659                                       NumCPU(numCPU),
0660                                       Strategy(2),
0661                                       Minos(*w->var("efficiency")),
0662                                       PrintLevel(quiet ? -1 : 1),
0663                                       PrintEvalErrors(quiet ? -1 : 1),
0664                                       Warnings(!quiet)));
0665   }
0666 
0667   // save everything
0668   res->Write("fitresults");
0669   w->saveSnapshot("finalState", w->components());
0670   saveFitPlot(w);
0671   //extract the efficiency parameter from the results
0672   RooRealVar* e = (RooRealVar*)res->floatParsFinal().find("efficiency");
0673   //What's wrong with this? and why don't they copy the errors!
0674   //efficiency = *e;
0675   efficiency.setVal(e->getVal());
0676   Double_t errLo = e->getErrorLo(), errHi = e->getErrorHi();
0677   if (errLo == 0 && e->getVal() < 0.5)
0678     errLo = e->getMin() - e->getVal();
0679   if (errHi == 0 && e->getVal() > 0.5)
0680     errHi = e->getMax() - e->getVal();
0681   efficiency.setAsymError(errLo, errHi);
0682 
0683   if (totPassing * totFailing == 0) {
0684     RooRealVar* nTot = (RooRealVar*)res->floatParsFinal().find("numTot");
0685     RooRealVar* fSig = (RooRealVar*)res->floatParsFinal().find("fSigAll");
0686     double cerr = ROOT::Math::beta_quantile(1 - (1.0 - .68540158589942957) / 2, 1, nTot->getVal() * fSig->getVal());
0687     /*
0688     std::cout << "======================================================================================" << std::endl;
0689     std::cout << "======= totPassing "  << totPassing << ", totFailing " << totFailing << std::endl;
0690     std::cout << "======= FIT: e  "  <<  e->getVal() << ",  e Lo " << e->getErrorLo()  << ",  e Hi " <<  e->getErrorHi() << std::endl;
0691     std::cout << "======= FIT:nS  "  << nS->getVal() << ", nS Lo " << nS->getErrorLo() << ", nS Hi " << nS->getErrorHi() << std::endl;
0692     std::cout << "======= FIT:nB  "  << nB->getVal() << ", nB Lo " << nB->getErrorLo() << ", nB Hi " << nB->getErrorHi() << std::endl;
0693     std::cout << "======= CNT:    "  << cerr << std::endl;
0694     std::cout << "======================================================================================" << std::endl;
0695     */
0696     if (totPassing == 0) {
0697       efficiency.setVal(0);
0698       efficiency.setAsymError(0, cerr);
0699     } else {
0700       efficiency.setVal(1);
0701       efficiency.setAsymError(-cerr, 0);
0702     }
0703   }
0704 
0705   delete simNLL;
0706 }
0707 
0708 void TagProbeFitter::createPdf(RooWorkspace* w, vector<string>& pdfCommands) {
0709   // create the signal and background pdfs defined by the user
0710   for (unsigned int i = 0; i < pdfCommands.size(); i++) {
0711     const std::string& command = pdfCommands[i];
0712     if (command.find("#import ") == 0) {
0713       TDirectory* here = gDirectory;
0714       w->import(command.substr(8).c_str());
0715       here->cd();
0716     } else {
0717       TDirectory* here = gDirectory;
0718       w->factory(command.c_str());
0719       here->cd();
0720     }
0721   }
0722   // setup the simultaneous extended pdf
0723 
0724   w->factory("expr::nSignalPass('efficiency*fSigAll*numTot', efficiency, fSigAll[.9,0,1],numTot[1,0,1e10])");
0725   w->factory("expr::nSignalFail('(1-efficiency)*fSigAll*numTot', efficiency, fSigAll,numTot)");
0726   w->factory("expr::nBkgPass('effBkg*(1-fSigAll)*numTot', effBkg[.9,0,1],fSigAll,numTot)");
0727   w->factory("expr::nBkgFail('(1-effBkg)*(1-fSigAll)*numTot', effBkg,fSigAll,numTot)");
0728   TString sPass = "signal", sFail = "signal";
0729   if (w->pdf("signalPass") != nullptr && w->pdf("signalFail") != nullptr) {
0730     if (w->pdf("signal") != nullptr)
0731       throw std::logic_error(
0732           "You must either define one 'signal' PDF or two PDFs ('signalPass', 'signalFail'), not both!");
0733     sPass = "signalPass";
0734     sFail = "signalFail";
0735   } else if (w->pdf("signal") != nullptr) {
0736     if (w->pdf("signalPass") != nullptr || w->pdf("signalFail") != nullptr) {
0737       throw std::logic_error(
0738           "You must either define one 'signal' PDF or two PDFs ('signalPass', 'signalFail'), not both!");
0739     }
0740   } else {
0741     throw std::logic_error("You must either define one 'signal' PDF or two PDFs ('signalPass', 'signalFail')");
0742   }
0743   w->factory("SUM::pdfPass(nSignalPass*" + sPass + ", nBkgPass*backgroundPass)");  //fBkgPass*
0744   w->factory("SUM::pdfFail(nSignalFail*" + sFail + ", nBkgFail*backgroundFail)");  //fBkgFail*
0745 
0746   w->factory("SIMUL::simPdf(_efficiencyCategory_, Passed=pdfPass, Failed=pdfFail)");
0747   // signalFractionInPassing is not used in the fit just to set the initial values
0748   if (w->pdf("simPdf") == nullptr)
0749     throw std::runtime_error("Could not create simultaneous fit pdf.");
0750   if (w->var("signalFractionInPassing") == nullptr)
0751     w->factory("signalFractionInPassing[0.9]");
0752 }
0753 
0754 void TagProbeFitter::setInitialValues(RooWorkspace* w) {
0755   // calculate initial values
0756   double totPassing = w->data("data")->sumEntries("_efficiencyCategory_==_efficiencyCategory_::Passed");
0757   double totFailing = w->data("data")->sumEntries("_efficiencyCategory_==_efficiencyCategory_::Failed");
0758   //std::cout << "Number of probes: " << totPassing+totFailing << std::endl;
0759 
0760   // now set the values
0761   w->var("numTot")->setVal(totPassing + totFailing);
0762   w->var("numTot")->setMax(2.0 * (totPassing + totFailing) + 10);  //wiggle room in case of 0 events in bin
0763 
0764   if (totPassing == 0) {
0765     w->var("efficiency")->setVal(0.0);
0766     w->var("efficiency")->setAsymError(0, 1);
0767     w->var("efficiency")->setConstant(false);
0768   } else if (totFailing == 0) {
0769     w->var("efficiency")->setVal(1.0);
0770     w->var("efficiency")->setAsymError(-1, 0);
0771     w->var("efficiency")->setConstant(false);
0772   } else {
0773     w->var("efficiency")->setConstant(false);
0774   }
0775 
0776   // save initial state for reference
0777   w->saveSnapshot("initialState", w->components());
0778 }
0779 
0780 void TagProbeFitter::saveFitPlot(RooWorkspace* w) {
0781   // save refferences for convenience
0782   RooCategory& efficiencyCategory = *w->cat("_efficiencyCategory_");
0783   RooAbsData* dataAll = (binnedFit ? w->data("data_binned") : w->data("data"));
0784   RooAbsData* dataPass = dataAll->reduce(Cut("_efficiencyCategory_==_efficiencyCategory_::Passed"));
0785   RooAbsData* dataFail = dataAll->reduce(Cut("_efficiencyCategory_==_efficiencyCategory_::Failed"));
0786   RooAbsPdf& pdf = *w->pdf("simPdf");
0787   std::unique_ptr<RooArgSet> obs(pdf.getObservables(*dataAll));
0788   RooRealVar* mass = nullptr;
0789   for (RooAbsArg* v : *obs) {
0790     if (!v->InheritsFrom("RooRealVar"))
0791       continue;
0792     mass = (RooRealVar*)v;
0793     break;
0794   }
0795   if (!mass)
0796     return;
0797   // make a 2x2 canvas
0798   TCanvas canvas("fit_canvas");
0799   canvas.Divide(2, 2);
0800   vector<RooPlot*> frames;
0801   // plot the Passing Probes
0802   canvas.cd(1);
0803   if (massBins == 0) {
0804     frames.push_back(mass->frame(Name("Passing"), Title("Passing Probes")));
0805     frames.push_back(mass->frame(Name("Failing"), Title("Failing Probes")));
0806     frames.push_back(mass->frame(Name("All"), Title("All Probes")));
0807   } else {
0808     frames.push_back(mass->frame(Name("Passing"), Title("Passing Probes"), Bins(massBins)));
0809     frames.push_back(mass->frame(Name("Failing"), Title("Failing Probes"), Bins(massBins)));
0810     frames.push_back(mass->frame(Name("All"), Title("All Probes"), Bins(massBins)));
0811   }
0812   dataPass->plotOn(frames[0]);
0813   pdf.plotOn(frames[0], Slice(efficiencyCategory, "Passed"), ProjWData(*dataPass), LineColor(kGreen));
0814   pdf.plotOn(frames[0],
0815              Slice(efficiencyCategory, "Passed"),
0816              ProjWData(*dataPass),
0817              LineColor(kGreen),
0818              Components("backgroundPass"),
0819              LineStyle(kDashed));
0820   frames[0]->Draw();
0821   // plot the Failing Probes
0822   canvas.cd(2);
0823   dataFail->plotOn(frames[1]);
0824   pdf.plotOn(frames[1], Slice(efficiencyCategory, "Failed"), ProjWData(*dataFail), LineColor(kRed));
0825   pdf.plotOn(frames[1],
0826              Slice(efficiencyCategory, "Failed"),
0827              ProjWData(*dataFail),
0828              LineColor(kRed),
0829              Components("backgroundFail"),
0830              LineStyle(kDashed));
0831   frames[1]->Draw();
0832   // plot the All Probes
0833   canvas.cd(3);
0834   dataAll->plotOn(frames.back());
0835   pdf.plotOn(frames.back(), ProjWData(*dataAll), LineColor(kBlue));
0836   pdf.plotOn(frames.back(),
0837              ProjWData(*dataAll),
0838              LineColor(kBlue),
0839              Components("backgroundPass,backgroundFail"),
0840              LineStyle(kDashed));
0841   frames.back()->Draw();
0842   // plot the Fit Results
0843   canvas.cd(4);
0844   frames.push_back(mass->frame(Name("Fit Results"), Title("Fit Results")));
0845   pdf.paramOn(
0846       frames.back(), RooFit::Label(""), RooFit::Format("NELU", AutoPrecision(0)), RooFit::Layout(0.1, 0.9, 0.9));
0847   // draw only the parameter box not the whole frame
0848   frames.back()->findObject(Form("%s_paramBox", pdf.GetName()))->Draw();
0849   //save and clean up
0850   canvas.Draw();
0851   canvas.Write();
0852   for (size_t i = 0; i < frames.size(); i++) {
0853     delete frames[i];
0854   }
0855   delete dataPass;
0856   delete dataFail;
0857 }
0858 
0859 void TagProbeFitter::saveDistributionsPlot(RooWorkspace* w) {
0860   // save pointers to datasets to manage memory
0861   RooAbsData* dataAll = w->data("data");
0862   RooAbsData* dataPass = dataAll->reduce(Cut("_efficiencyCategory_==_efficiencyCategory_::Passed"));
0863   RooAbsData* dataFail = dataAll->reduce(Cut("_efficiencyCategory_==_efficiencyCategory_::Failed"));
0864 
0865   const RooArgSet* vars = dataAll->get();
0866   vector<RooRealVar*> reals;
0867   for (RooAbsArg* v : *vars) {
0868     if (!v->InheritsFrom("RooRealVar"))
0869       continue;
0870     reals.push_back((RooRealVar*)v);
0871   }
0872   TCanvas canvas("distributions_canvas");
0873   canvas.Divide(3, reals.size());
0874   vector<RooPlot*> frames;
0875   for (unsigned int i = 0; i < reals.size(); i++) {
0876     // plot the Passing Probes
0877     canvas.cd(3 * i + 1);
0878     frames.push_back(reals[i]->frame(Name("Passing"), Title("Passing Probes"), Bins(100)));
0879     dataPass->plotOn(frames.back(), LineColor(kGreen));
0880     dataPass->statOn(frames.back());
0881     frames.back()->Draw();
0882     // plot the Failing Probes
0883     canvas.cd(3 * i + 2);
0884     frames.push_back(reals[i]->frame(Name("Failing"), Title("Failing Probes"), Bins(100)));
0885     dataFail->plotOn(frames.back(), LineColor(kRed));
0886     dataFail->statOn(frames.back());
0887     frames.back()->Draw();
0888     // plot the All Probes
0889     canvas.cd(3 * i + 3);
0890     frames.push_back(reals[i]->frame(Name("All"), Title("All Probes"), Bins(100)));
0891     dataAll->plotOn(frames.back(), LineColor(kBlue));
0892     dataAll->statOn(frames.back());
0893     frames.back()->Draw();
0894   }
0895   canvas.Draw();
0896   canvas.Write();
0897   for (size_t i = 0; i < frames.size(); i++) {
0898     delete frames[i];
0899   }
0900   delete dataPass;
0901   delete dataFail;
0902 }
0903 
0904 void TagProbeFitter::saveEfficiencyPlots(RooDataSet& eff,
0905                                          const TString& effName,
0906                                          RooArgSet& binnedVariables,
0907                                          RooArgSet& mappedCategories) {
0908   bool isOnePoint =
0909       (eff.numEntries() == 1);  // for datasets with > 1 entry, we don't make plots for variables with just one bin
0910   for (auto it1 = binnedVariables.begin(); it1 != binnedVariables.end(); it1++) {
0911     RooRealVar* v1 = dynamic_cast<RooRealVar*>(*it1);
0912     RooArgSet binCategories1D;
0913     if (v1->numBins() == 1 && !isOnePoint)
0914       continue;
0915     for (auto it2 = binnedVariables.begin(); it2 != binnedVariables.end(); it2++) {
0916       RooRealVar* v2 = dynamic_cast<RooRealVar*>(*it2);
0917       if (v2 == v1)
0918         continue;
0919       if (v2->numBins() == 1 && !isOnePoint)
0920         continue;
0921       binCategories1D.addClone(
0922           RooBinningCategory(TString(v2->GetName()) + "_bins", TString(v2->GetName()) + "_bins", *v2));
0923 
0924       RooArgSet binCategories2D;
0925       for (auto it3 = binnedVariables.begin(); it3 != binnedVariables.end(); it3++) {
0926         RooRealVar* v3 = dynamic_cast<RooRealVar*>(*it3);
0927         if (v3 == v1 || v3 == v2)
0928           continue;
0929         binCategories2D.addClone(
0930             RooBinningCategory(TString(v3->GetName()) + "_bins", TString(v3->GetName()) + "_bins", *v3));
0931       }
0932       RooMultiCategory allCats2D("allCats2D", "allCats2D", RooArgSet(binCategories2D, mappedCategories));
0933       if (allCats2D.numTypes() == 0) {
0934         makeEfficiencyPlot2D(eff, *v1, *v2, TString::Format("%s_%s_PLOT", v1->GetName(), v2->GetName()), "", effName);
0935       } else {
0936         RooDataSet myEff(eff);
0937         myEff.addColumn(allCats2D);
0938         std::unique_ptr<TIterator> catIt(allCats2D.typeIterator());
0939         for (RooCatType* t = (RooCatType*)catIt->Next(); t != nullptr; t = (RooCatType*)catIt->Next()) {
0940           TString catName = t->GetName();
0941           if (catName.Contains("NotMapped"))
0942             continue;
0943           catName.ReplaceAll("{", "").ReplaceAll("}", "").ReplaceAll(";", "_&_");
0944           makeEfficiencyPlot2D(myEff,
0945                                *v1,
0946                                *v2,
0947                                TString::Format("%s_%s_PLOT_%s", v1->GetName(), v2->GetName(), catName.Data()),
0948                                catName,
0949                                effName,
0950                                "allCats2D",
0951                                t->getVal());
0952         }
0953       }
0954     }
0955     RooMultiCategory allCats1D("allCats1D", "allCats1D", RooArgSet(binCategories1D, mappedCategories));
0956     if (allCats1D.numTypes() == 0) {
0957       makeEfficiencyPlot1D(eff, *v1, TString::Format("%s_PLOT", v1->GetName()), "", effName);
0958     } else {
0959       RooDataSet myEff(eff);
0960       myEff.addColumn(allCats1D);
0961       std::unique_ptr<TIterator> catIt(allCats1D.typeIterator());
0962       for (RooCatType* t = (RooCatType*)catIt->Next(); t != nullptr; t = (RooCatType*)catIt->Next()) {
0963         TString catName = t->GetName();
0964         if (catName.Contains("NotMapped"))
0965           continue;
0966         catName.ReplaceAll("{", "").ReplaceAll("}", "").ReplaceAll(";", "_&_");
0967         makeEfficiencyPlot1D(myEff,
0968                              *v1,
0969                              TString::Format("%s_PLOT_%s", v1->GetName(), catName.Data()),
0970                              catName,
0971                              effName,
0972                              "allCats1D",
0973                              t->getVal());
0974       }
0975     }
0976   }
0977 }
0978 
0979 void TagProbeFitter::makeEfficiencyPlot1D(RooDataSet& eff,
0980                                           RooRealVar& v,
0981                                           const TString& plotName,
0982                                           const TString& plotTitle,
0983                                           const TString& effName,
0984                                           const char* catName,
0985                                           int catIndex) {
0986   TGraphAsymmErrors* p = new TGraphAsymmErrors();
0987   const RooArgSet* entry = eff.get();
0988   const RooRealVar& vi = dynamic_cast<const RooRealVar&>(*entry->find(v.GetName()));
0989   const RooRealVar& ei = dynamic_cast<const RooRealVar&>(*entry->find("efficiency"));
0990   for (unsigned int i = 0, n = eff.numEntries(); i < n; ++i) {
0991     entry = eff.get(i);
0992     if (catName != nullptr && entry->getCatIndex(catName) != catIndex)
0993       continue;
0994     int j = p->GetN();
0995     p->Set(j + 1);
0996     p->SetPoint(j, vi.getVal(), ei.getVal());
0997     p->SetPointError(j, -vi.getAsymErrorLo(), vi.getAsymErrorHi(), -ei.getAsymErrorLo(), ei.getAsymErrorHi());
0998   }
0999   TCanvas canvas(plotName);
1000   TH1F* frame = new TH1F("frame", "Efficiency of " + effName, 1, v.getMin(), v.getMax());
1001   frame->SetDirectory(nullptr);
1002   p->SetNameTitle(Form("hxy_%s", eff.GetName()), "Efficiency of " + effName);
1003   p->GetXaxis()->SetTitle(strlen(v.getUnit()) ? Form("%s (%s)", v.GetName(), v.getUnit()) : v.GetName());
1004   p->GetYaxis()->SetTitle("Efficiency of " + effName);
1005   frame->GetXaxis()->SetTitle(strlen(v.getUnit()) ? Form("%s (%s)", v.GetName(), v.getUnit()) : v.GetName());
1006   frame->GetYaxis()->SetTitle("Efficiency of " + effName);
1007   frame->GetYaxis()->SetRangeUser(0, 1);
1008   frame->Draw();
1009   p->SetLineWidth(2);
1010   p->SetMarkerStyle(kFullCircle);
1011   p->SetMarkerSize(1.2);
1012   p->Draw("P SAME");
1013   canvas.Write();
1014   delete frame;
1015   delete p;
1016 }
1017 
1018 void TagProbeFitter::makeEfficiencyPlot2D(RooDataSet& eff,
1019                                           RooRealVar& v1,
1020                                           RooRealVar& v2,
1021                                           const TString& plotName,
1022                                           const TString& plotTitle,
1023                                           const TString& effName,
1024                                           const char* catName,
1025                                           int catIndex) {
1026   TCanvas canvas(plotName);
1027   canvas.SetRightMargin(0.15);
1028   TH2F* h = new TH2F(plotName,
1029                      plotName,
1030                      v1.getBinning().numBins(),
1031                      v1.getBinning().array(),
1032                      v2.getBinning().numBins(),
1033                      v2.getBinning().array());
1034   const RooArgSet* set = eff.get();
1035   RooRealVar* e = (RooRealVar*)set->find("efficiency");
1036   RooRealVar* v1_ = (RooRealVar*)set->find(v1.GetName());
1037   RooRealVar* v2_ = (RooRealVar*)set->find(v2.GetName());
1038   h->SetTitle(TString::Format("%s;%s%s;%s%s;Efficiency of %s",
1039                               plotTitle.Data(),
1040                               v1.GetTitle(),
1041                               TString(v1.getUnit()).Length() == 0 ? "" : TString::Format(" (%s)", v1.getUnit()).Data(),
1042                               v2.GetTitle(),
1043                               TString(v2.getUnit()).Length() == 0 ? "" : TString::Format(" (%s)", v2.getUnit()).Data(),
1044                               effName.Data()));
1045   h->SetOption("colztexte");
1046   h->GetZaxis()->SetRangeUser(-0.001, 1.001);
1047   h->SetStats(kFALSE);
1048   for (int i = 0; i < eff.numEntries(); i++) {
1049     const RooArgSet* entry = eff.get(i);
1050     if (catName != nullptr && entry->getCatIndex(catName) != catIndex)
1051       continue;
1052     h->SetBinContent(h->FindBin(v1_->getVal(), v2_->getVal()), e->getVal());
1053     h->SetBinError(h->FindBin(v1_->getVal(), v2_->getVal()), (e->getErrorHi() - e->getErrorLo()) / 2.);
1054   }
1055   h->Draw();
1056   canvas.Draw();
1057   canvas.Write();
1058   delete h;
1059 }
1060 
1061 void TagProbeFitter::doSBSEfficiency(RooWorkspace* w, RooRealVar& efficiency) {}
1062 
1063 void TagProbeFitter::doCntEfficiency(RooWorkspace* w, RooRealVar& efficiency) {
1064   int pass = w->data("data")->sumEntries("_efficiencyCategory_==_efficiencyCategory_::Passed");
1065   int fail = w->data("data")->sumEntries("_efficiencyCategory_==_efficiencyCategory_::Failed");
1066   double e = (pass + fail == 0) ? 0 : pass / double(pass + fail);
1067   // Use Clopper-Pearson
1068   double alpha = (1.0 - .68540158589942957) / 2;
1069   double lo = (pass == 0) ? 0.0 : ROOT::Math::beta_quantile(alpha, pass, fail + 1);
1070   double hi = (fail == 0) ? 1.0 : ROOT::Math::beta_quantile(1 - alpha, pass + 1, fail);
1071   ////from TGraphAsymmErrors
1072   //double lob, hib;
1073   //Efficiency( pass, pass+fail, .68540158589942957, e, lob, hib );
1074   //std::cerr << "CNT " << pass << "/" << fail << ":  Clopper Pearson [" << lo << ", "  << hi << "], Bayes [" << lob << ", " << hib << "]" << std::endl;
1075   efficiency.setVal(e);
1076   efficiency.setAsymError(lo - e, hi - e);
1077 }
1078 
1079 void TagProbeFitter::varFixer(RooWorkspace* w, bool fix) {
1080   std::vector<std::string>::const_iterator it;
1081   for (it = fixVars.begin(); it < fixVars.end(); it++) {
1082     if (w->var((*it).c_str()))
1083       w->var((*it).c_str())->setAttribute("Constant", fix);
1084     else {
1085       std::cout << "TagProbeFitter: "
1086                 << "Can't find a variable to fix: " << *it;
1087     }
1088   }
1089 }
1090 
1091 void TagProbeFitter::varSaver(RooWorkspace* w) {
1092   if (!fixVarValues.empty()) {
1093     std::cout << "attempt to save variables more than once!" << std::endl;
1094     return;
1095   }
1096   std::vector<std::string>::const_iterator it;
1097   for (it = fixVars.begin(); it < fixVars.end(); it++) {
1098     fixVarValues.push_back(w->var((*it).c_str())->getVal());
1099   }
1100 }
1101 
1102 void TagProbeFitter::varRestorer(RooWorkspace* w) {
1103   if (fixVarValues.size() == fixVars.size())
1104     for (unsigned int i = 0; i < fixVars.size(); i++) {
1105       std::cout << "setting variable " << fixVars[i].c_str() << std::endl;
1106       w->var(fixVars[i].c_str())->setVal(fixVarValues[i]);
1107     }
1108   else {
1109     std::cout << "fixVars and fixVarValues are not of the same size!" << std::endl;
1110   }
1111 }