testLikelihoodDiffPlots.C

Go to the documentation of this file.

//
//    Copyright 2010
//
//    This file is part of rootpwa
//
//    rootpwa is free software: you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//
//    rootpwa is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
//
//-------------------------------------------------------------------------
// File and Version Information:
// $Rev::                             $: revision of last commit
// $Author::                          $: author of last commit
// $Date::                            $: date of last commit
//
// Description:
//      generates difference plots for two trees with log likelihood
//      values and corresponding derivatives generated by
//      testLikelihood
//
//
// Author List:
//      Boris Grube          TUM            (original author)
//
//
//-------------------------------------------------------------------------


#include <iostream>
#include <sstream>
#include <string>
#include <map>
#include <algorithm>
#include <cmath>

#include "TFile.h"
#include "TTree.h"
#include "TH1.h"


using namespace std;


void testLikelihoodDiffPlots(const string& inFileNameA = "testLikelihood.root",
                             const string& inFileNameB = "testLikelihood.root",
                             const string& outFileName = "testLikelihoodDiff.root",
                             const string& inTreeName  = "testLikelihoodTree")
{
        // open input files, get trees, and read number of parameters
        const string inFileNames[2] = {inFileNameA, inFileNameB};
        TFile*       inFiles[2]     = {0, 0};
        TTree*       inTrees[2]     = {0, 0};
        // leaf variables
        Double_t             logLikelihood[2] = {0, 0};
        map<string, double>* derivatives  [2] = {0, 0};
        for (unsigned int i = 0; i < 2; ++i) {
                inFiles[i] = TFile::Open(inFileNames[i].c_str(), "READ");
                if (!inFiles[i] || inFiles[i]->IsZombie()) {
                        cout << "cannot open file '" << inFileNames[i] << "'. aborting." << endl;
                        throw;
                }
                inFiles[i]->GetObject(inTreeName.c_str(), inTrees[i]);
                if (!inTrees[i]) {
                        cout << "cannot find tree '" << inTreeName << "'. aborting." << endl;
                        throw;
                }
                // connect leafs
                inTrees[i]->SetBranchAddress("logLikelihood", &logLikelihood[i]);
                inTrees[i]->SetBranchAddress("derivatives",   &derivatives  [i]);
                inTrees[i]->GetEntry(0);
        }
        if (derivatives[0]->size() != derivatives[1]->size()) {
                cout << "derivative arrays have different size: "
                     << derivatives[0]->size() << " vs. " << derivatives[1]->size() << ". aborting." << endl;
                throw;
        }

        // create output file
        TFile* outFile = TFile::Open(outFileName.c_str(), "RECREATE");
        if (!outFile || outFile->IsZombie()) {
                cout << "cannot open file '" << outFileName << "'. aborting." << endl;
                throw;
        }

        // define histograms
        TH1D* hLikelihoodDiffAbs = new TH1D
                ("hLikelihoodDiffAbs", "hLikelihoodDiffAbs;L_{CPU} - L_{GPU};Count", 10000, -1e-7, 1e-7);
        TH1D* hLikelihoodDiffRel = new TH1D
                ("hLikelihoodDiffRel", "hLikelihoodDiffRel;1 - L_{GPU} / L_{CPU};Count", 10000, -1e-11, 1e-11);
        TH1D* hDerivDiffTotAbs = new TH1D
                ("hDerivDiffTotAbs", "hDerivDiffTotAbs;#nabla_{CPU} - #nabla_{GPU};Count",
                 10000, -1e-10, 1e-10);
        TH1D* hDerivDiffTotRel = new TH1D
                ("hDerivDiffTotRel", "hDerivDiffTotRel;1 - #nabla_{GPU} / #nabla_{CPU};Count",
                 10000, -1e-10, 1e-10);
        map<string, TH1D*> hDerivDiffAbs;
        map<string, TH1D*> hDerivDiffRel;
        for (map<string, double>::const_iterator i = derivatives[0]->begin();
             i != derivatives[0]->end(); ++i) {
                stringstream suf;
                suf << "_" << i->first;
                hDerivDiffAbs[i->first] =
                        new TH1D(("hDerivDiffAbs" + suf.str()).c_str(),
                                 ("hDerivDiffAbs" + suf.str() + ";#nabla_{CPU} - #nabla_{GPU};Count").c_str(),
                                 10000, -1e-11, 1e-11);
                hDerivDiffRel[i->first] =
                        new TH1D(("hDerivDiffRel" + suf.str()).c_str(),
                                 ("hDerivDiffRel" + suf.str() + ";1 - #nabla_{GPU} / #nabla_{CPU};Count").c_str(),
                                 10000, -1e-11, 1e-11);
        }

        // loop over tree and fill histograms
        // long int nmbEntries = min(inTrees[0]->GetEntries(), inTrees[1]->GetEntries());
        long int nmbEntries = 100;
        for (long int entry = 0; entry < nmbEntries; ++entry) {
                Long64_t treeEntries[2] = {inTrees[0]->LoadTree(entry), inTrees[1]->LoadTree(entry)};
                inTrees[0]->GetEntry(treeEntries[0]);
                inTrees[1]->GetEntry(treeEntries[1]);
                
                hLikelihoodDiffAbs->Fill(logLikelihood[0] - logLikelihood[1]);
                hLikelihoodDiffRel->Fill(1 - logLikelihood[1] / logLikelihood[0]);
                // cout << "[" << entry << "] = " << logLikelihood[0] << " vs. " << logLikelihood[1] << endl;
                for (map<string, double>::const_iterator i = derivatives[0]->begin();
                     i != derivatives[0]->end(); ++i) {
                        map<string, double>::const_iterator j  = derivatives[1]->find(i->first);
                        if (j == derivatives[1]->end()) {
                                cout << "cannot find derivative for '" << i->first << "' in second tree. skipping." << endl;
                                continue;
                        }
                        double diffAbs = i->second - j->second;
                        double diffRel = 1 - j->second / i->second;
                        hDerivDiffAbs[i->first]->Fill(diffAbs);
                        hDerivDiffRel[i->first]->Fill(diffRel);
                        hDerivDiffTotAbs->Fill(diffAbs);
                        hDerivDiffTotRel->Fill(diffRel);
                }
        }

        // find gradients with maximum mean and RMS of relative difference
        double maxMean = 0;
        string maxMeanParName;
        double maxRms = 0;
        string maxRmsParName;
        for (map<string, TH1D*>::const_iterator i = hDerivDiffRel.begin(); i != hDerivDiffRel.end(); ++i) {
                const double mean = abs((double)i->second->GetMean());
                if (mean > maxMean) {
                        maxMean        = mean;
                        maxMeanParName = i->first;
                }
                const double rms = i->second->GetRMS();
                if (rms > maxRms) {
                        maxRms        = rms;
                        maxRmsParName = i->first;
                }
        }
        cout << "max. mean of relative diff. = " << maxMean << " @ [" << maxMeanParName << "]" << endl;
        cout << "max. RMS  of relative diff. = " << maxRms  << " @ [" << maxRmsParName  << "]" << endl;

        // write outout and cleanup
        inFiles[0]->Close();
        inFiles[1]->Close();
        outFile->Write();
        outFile->Close();
}