MVAComputer/bin/mvaExtractPDFs.cpp

0001 #include <iostream>
0002 #include <sstream>
0003 #include <string>
0004 #include <vector>
0005 #include <map>
0006
0007 #include <TFile.h>
0008 #include <TH1.h>
0009
0010 #include "PhysicsTools/MVAComputer/interface/Calibration.h"
0011 #include "PhysicsTools/MVAComputer/interface/MVAComputer.h"
0012 #include "PhysicsTools/MVAComputer/interface/Spline.h"
0013
0014 using namespace PhysicsTools;
0015
0016 static const int precision = 100;
0017
0018 int main(int argc, char **argv) {
0019   using Calibration::HistogramF;
0020   using Calibration::VarProcessor;
0021
0022   if (argc != 3) {
0023     std::cerr << "Syntax: " << argv[0] << " <MVA File> "
0024               << "<output ROOT file>" << std::endl;
0025     return 1;
0026   }
0027
0028   Calibration::MVAComputer *calib = MVAComputer::readCalibration(argv[1]);
0029   if (!calib)
0030     return 1;
0031
0032   std::map<std::string, HistogramF *> histos;
0033
0034   std::vector<VarProcessor *> procs = calib->getProcessors();
0035   for (unsigned int z = 0; z < procs.size(); ++z) {
0036     VarProcessor *proc = procs[z];
0037     if (!proc)
0038       continue;
0039
0040     std::ostringstream ss3;
0041     ss3 << (z + 1);
0042
0043     Calibration::ProcLikelihood *lkh = dynamic_cast<Calibration::ProcLikelihood *>(proc);
0044     Calibration::ProcNormalize *norm = dynamic_cast<Calibration::ProcNormalize *>(proc);
0045
0046     if (lkh) {
0047       for (unsigned int i = 0; i < lkh->pdfs.size(); i++) {
0048         std::ostringstream ss2;
0049         ss2 << (i + 1);
0050         histos["proc" + ss3.str() + "_sig" + ss2.str()] = &lkh->pdfs[i].signal;
0051         histos["proc" + ss3.str() + "_bkg" + ss2.str()] = &lkh->pdfs[i].background;
0052       }
0053     } else if (norm) {
0054       for (unsigned int i = 0; i < norm->distr.size(); i++) {
0055         std::ostringstream ss2;
0056         ss2 << (i + 1);
0057         histos["proc" + ss3.str() + "_norm" + ss2.str()] = &norm->distr[i];
0058       }
0059     }
0060   }
0061
0062   TFile *f = TFile::Open(argv[2], "RECREATE");
0063   if (!f)
0064     return 2;
0065
0066   for (std::map<std::string, HistogramF *>::const_iterator iter = histos.begin(); iter != histos.end(); ++iter) {
0067     std::string name = iter->first;
0068     HistogramF *histo = iter->second;
0069
0070     unsigned int size = histo->values().size() - 2;
0071     std::vector<double> values(histo->values().begin() + 1, histo->values().end() - 1);
0072     Spline spline;
0073     spline.set(values.size(), &values.front());
0074
0075     double min = histo->range().min;
0076     double max = histo->range().max;
0077
0078     TH1F *h = new TH1F((name + "_histo").c_str(),
0079                        (name + "_histo").c_str(),
0080                        size,
0081                        min - 0.5 * (max - min) / size,
0082                        max + 0.5 * (max - min) / size);
0083     TH1F *s = new TH1F((name + "_spline").c_str(), (name + "_spline").c_str(), size * precision, min, max);
0084
0085     for (unsigned int i = 0; i < size; i++) {
0086       h->SetBinContent(i + 1, histo->values()[i + 1]);
0087       for (int j = 0; j < precision; j++) {
0088         unsigned int k = i * precision + j;
0089         double x = (k + 0.5) / (size * precision);
0090         double v = spline.eval(x);
0091         s->SetBinContent(k, v);
0092       }
0093     }
0094   }
0095
0096   f->Write();
0097   delete f;
0098
0099   return 0;
0100 }