/Users/jmonk/Physics/ForIA/src/D3PD/RootHistogrammer.cxx

Go to the documentation of this file.

#include "ForIA/D3PD/RootHistogrammer.hh"

#include "TH1D.h"
#include "TH2D.h"
#include "TH3D.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TProfile.h"

#include <iostream>

namespace ForIA {

  RootHistogrammer::RootHistogrammer(const string &rootFileName) : IHistogrammer(), m_rootFileName(rootFileName), m_rootDir("/"){
    m_directories["/"] = &m_rootDir;
    TH1::AddDirectory(false);
  }
  
  const string &RootHistogrammer::outputName()const{
    return m_rootFileName;
  }
  
  void RootHistogrammer::bookHistogram1D(const string &path, const string &name, 
                                         const string &title, 
                                         int nBins, double xMin, double xMax,
                                         bool divideBinWidths){
    
    if(xMin > xMax){
      throw BinError(xMin, xMax);
    }
    
    TH1D *histo = new TH1D(name.c_str(), title.c_str(), nBins, xMin, xMax);
    histo->Sumw2();
    histo->ResetStats();
    histo->Reset();
    
    try{
      addHisto(path, name, histo, divideBinWidths);
    }catch(DuplicateHistogramException er){
      delete histo;
      throw er;
    }
        
    return; 
  }
  
  void RootHistogrammer::bookHistogram1D(const string &path, const string &name, 
                                         const string &title, const vector<double> &bins,
                                         bool divideBinWidths){

    TH1D *histo = new TH1D(name.c_str(), title.c_str(), bins.size() - 1, &(bins[0]));
    histo->Sumw2();
    histo->ResetStats();
    histo->Reset();
    
    try{
      addHisto(path, name, histo, divideBinWidths);
    }catch(DuplicateHistogramException er){
      delete histo;
      throw er;
    }
    
    return;
  }

  void RootHistogrammer::bookProfile1D(const string &path, const string &name, 
                                       const string &title, int nBins, 
                                       double xMin, double xMax){
    
    TProfile *prof = new TProfile(name.c_str(), title.c_str(), nBins, xMin, xMax);
    prof->Sumw2();    
    try{
      addProfile(path, name, prof);
    }catch(DuplicateHistogramException er){
      delete prof;
      throw er;
    }
    
    return;
  }
  
  void RootHistogrammer::bookHistogram2D(const string &path, const string &name, const string &title, 
                                         int nBinsX, double xMin, double xMax,
                                         int nBinsY, double yMin, double yMax,
                                         bool divideBinWidths){
   
    if(xMin > xMax){
      throw BinError(xMin, xMax);
    }
    
    if(yMin > yMax){
      throw BinError(xMin, xMax);
    }
    
    TH2D *histo = new TH2D(name.c_str(), title.c_str(), nBinsX, xMin, xMax, nBinsY, yMin, yMax);
    histo->Sumw2();
    
    try{
      addHisto(path, name, histo, divideBinWidths);
    }catch(DuplicateHistogramException er){
      delete histo;
      throw er;
    }
        
    return ;
  }
  
  void RootHistogrammer::bookHistogram2D(const string &path, const string &name, const string &title,
                                         const vector<double> &xbins, const vector<double> &ybins,
                                         bool divideBinWidths){
   
    TH2D *histo = new TH2D(name.c_str(), title.c_str(), 
                           xbins.size() - 1, &(xbins[0]), ybins.size() - 1, &(ybins[0]));
    histo->Sumw2();
    
    try{
      addHisto(path, name, histo, divideBinWidths);
    }catch(DuplicateHistogramException er){
    
      delete histo;
      throw er;
    }
    return;
  }
  
  void RootHistogrammer::bookHistogram3D(const string &path, const string &name, const string &title, 
                          int nBinsX, double xMin, double xMax,
                          int nBinsY, double yMin, double yMax,
                          int nBinsZ, double zMin, double zMax,
                          bool divideBinWidths){
    
    if(xMin > xMax){
      throw BinError(xMin, xMax);
    }
    
    if(yMin > yMax){
      throw BinError(xMin, xMax);
    }
  
    if(zMin > zMax){
      throw BinError(zMin, zMax);
    }
    
    TH3D *histo = new TH3D(name.c_str(), title.c_str(), 
                           nBinsX, xMin, xMax, 
                           nBinsY, yMin, yMax,
                           nBinsZ, zMin, zMax);
    histo->Sumw2();
    
    try{
      addHisto(path, name, histo, divideBinWidths);
    }catch(DuplicateHistogramException er){
      
      delete histo;
      throw er;
    }
    return;
  }
  void RootHistogrammer::fillHistogram1D(const string &path, const string &name, 
                                         double value, double weight){
    
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name, this);
    
    TH1D *histo = dir->second->histogram(name);
    
    if(histo == 0) throw NoSuchHistogramException(path, name, this);
    
    histo->Fill(value, weight);
    
    return;
  }

  void RootHistogrammer::fillProfile1D(const string &path, const string &name,
                                       double xVal, double yVal, double weight){
  
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    
    if(dir == m_directories.end())throw NoSuchHistogramException(path, name, this);
    
    TProfile *prof = dir->second->profile(name);
    
    if(prof == 0) throw NoSuchHistogramException(path, name, this);

    prof->Fill(xVal, yVal, weight);
    
    return;
  }
  
  void RootHistogrammer::fillHistogram2D(const string &path, const string &name, 
                                         double xVal, double yVal, double weight){
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name, this);
    
    TH2D *histo = dir->second->histogram2D(name);
    if(histo ==0) throw NoSuchHistogramException(path, name, this);
      
    histo->Fill(xVal, yVal, weight);
    return;
  }
  
  void RootHistogrammer::fillHistogram3D(const string &path, const string &name, 
                                         double xVal, double yVal, double zVal,
                                         double weight){
  
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name, this);
    
    TH3D *histo = dir->second->histogram3D(name);
    if(histo ==0) throw NoSuchHistogramException(path, name, this);
    
    histo->Fill(xVal, yVal, zVal, weight);
    return;
  }
  void RootHistogrammer::finalise(){
   
    
    for(boost::unordered_map<string, Directory*>::iterator dir = m_directories.begin();
        dir != m_directories.end(); ++dir){
            
      for(Directory::iterator histo = dir->second->begin(); histo != dir->second->end(); ++histo){
        
        if(histo->second == 0) throw NoSuchHistogramException(dir->first, histo->first);
        
        if(! dir->second->doDivideBinWidths(histo->second)){
          std::cout<<"Histogram "<< histo->first<<" is not normalised by bin widths"<<std::endl;
          continue;
        }
        
        double nEntries = histo->second->GetEntries();
        
        int nBins = histo->second->GetNbinsX() +1;
        for(int bin=1; bin != nBins; ++bin){
          double area = histo->second->GetBinContent(bin);
          double width = histo->second->GetBinWidth(bin);
          double error = histo->second->GetBinError(bin) / width;
          double height = area / width;
          histo->second->SetBinContent(bin, height);
          histo->second->SetBinError(bin, error);
        }
        
        histo->second->SetEntries(nEntries);
        
      }
      
      for(Directory::iterator2D histo = dir->second->begin2D(); histo != dir->second->end2D(); ++histo){
        
        if(histo->second == 0) throw NoSuchHistogramException(dir->first, histo->first);
        
        if(! dir->second->doDivideBinWidths(histo->second)){
          std::cout<<"Histogram "<< histo->first<<" is not normalised by bin widths"<<std::endl;
          continue;
        }
        
        int nBinsX = histo->second->GetNbinsX() + 1;
        int nBinsY = histo->second->GetNbinsY() + 1;
        
        double nEntries = histo->second->GetEntries();
        
        for(int xbin =1; xbin != nBinsX; ++xbin){
          for(int ybin = 1; ybin != nBinsY; ++ybin){
            int bin = histo->second->GetBin(xbin, ybin);
            double volume = histo->second->GetBinContent(bin);
            double xWidth = histo->second->GetXaxis()->GetBinWidth(xbin);
            double yWidth = histo->second->GetYaxis()->GetBinWidth(ybin);
            double area = xWidth * yWidth;
            double height = volume / area;
            double error = histo->second->GetBinError(bin) / area;
            histo->second->SetBinContent(bin, height);
            histo->second->SetBinError(bin, error);
          }
        }
        histo->second->SetEntries(nEntries);
      }
     
      for(Directory::iterator3D histo = dir->second->begin3D(); histo != dir->second->end3D(); ++histo){
        
        if(histo->second == 0) throw NoSuchHistogramException(dir->first, histo->first);
        
        if(! dir->second->doDivideBinWidths(histo->second)){
          std::cout<<"Histogram "<<histo->first<<" is not normalised by bin widths"<<std::endl;
          continue;
        }
        
        int nBinsX = histo->second->GetNbinsX() + 1;
        int nBinsY = histo->second->GetNbinsY() + 1;
        int nBinsZ = histo->second->GetNbinsZ() + 1;
        
        double nEntries = histo->second->GetEntries();
        
        for(int xbin = 1; xbin != nBinsX; ++xbin){
          for(int ybin = 1; ybin != nBinsY; ++ybin){
            for(int zbin = 1; zbin != nBinsZ; ++zbin){
              int bin = histo->second->GetBin(xbin, ybin, zbin);
              double volume = histo->second->GetBinContent(bin);
              double xWidth = histo->second->GetXaxis()->GetBinWidth(bin);
              double yWidth = histo->second->GetYaxis()->GetBinWidth(bin);
              double zWidth = histo->second->GetZaxis()->GetBinWidth(bin);
              double dxdydz = xWidth * yWidth * zWidth;
              double height = volume / dxdydz;
              double error = histo->second->GetBinError(bin) / dxdydz;
              histo->second->SetBinContent(bin, height);
              histo->second->SetBinError(bin, error);
            }
          }
        }
        histo->second->SetEntries(nEntries);
      }
      
    }
    
    string filename = m_rootFileName + ".root";
    
    TFile *file = new TFile(filename.c_str(), "RECREATE");
    
    m_rootDir.write(file);
    
    file->Close();
    return;
  }
  
  void RootHistogrammer::normaliseHistogram1D(const string &path, const string &name, double norm){
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name);
    
    TH1D *histo = dir->second->histogram(name);
    
    if(histo == 0) throw NoSuchHistogramException(path, name);

    double intgl = histo->Integral();//"width");
    histo->Scale(norm / intgl);
    
    return;
  }
  
  void RootHistogrammer::scaleHistogram1D(const string &path, const string &name, double scale){
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name);
    
    TH1D *histo = dir->second->histogram(name);
    
    if(histo == 0) throw NoSuchHistogramException(path, name);
    
    histo->Scale(scale);
    
    return;
  }
  void RootHistogrammer::normaliseHistogram2D(const string &path, const string &name, double norm){
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name);
    
    TH2D *histo = dir->second->histogram2D(name);
    
    if(histo == 0) throw NoSuchHistogramException(path, name);
    
    double intgl = histo->Integral();//"width");
    histo->Scale(norm / intgl);
    
    return;
  }
  void RootHistogrammer::scaleHistogram2D(const string &path, const string &name, double scale){
    boost::unordered_map<string, Directory*>::iterator dir = m_directories.find(path);
    
    if(dir == m_directories.end()) throw NoSuchHistogramException(path, name);
    
    TH2D *histo = dir->second->histogram2D(name);
    
    if(histo == 0) throw NoSuchHistogramException(path, name);
    
    histo->Scale(scale);
    
    return;
  }
  
  void RootHistogrammer::showHistograms(std::ostream &out)const{
  
    for(boost::unordered_map<string, Directory*>::const_iterator dir = m_directories.begin();
        dir != m_directories.end(); ++dir){
      
      for(Directory::const_iterator histo = dir->second->begin();
          histo != dir->second->end(); ++histo){
        out<<dir->first<<histo->first<<std::endl;
      }
      
      for(Directory::const_iterator2D histo = dir->second->begin2D();
          histo != dir->second->end2D(); ++histo){
        out<<dir->first<<histo->first<<std::endl;
      }
      
    }
        
    return;
  }
  
  Directory::Directory(const string &name): 
  m_name(name){
    stripSlash(m_name);
  }
  
  const string &Directory::name()const{
    return m_name;
  }
  
  pair<string, Directory*> Directory::addDirectory(string path){
   
    stripSlash(path);
    
    size_t pos = path.find("/");
    
    string dirName;
    
    bool addMore = false;
    string nextDir = "";
    
    if(pos == string::npos){
      dirName = path;
    }else{
      dirName = path.substr(0, pos);
      nextDir = path.substr(pos);
      addMore = true;
    }
    
    map<string, Directory>::iterator dir = m_dirs.find(dirName);
    if(dir==m_dirs.end()){
      Directory newDir(dirName);
      dir = m_dirs.insert(std::pair<string, Directory>(dirName, newDir)).first;
    }
    
    string resultName = "/" + dirName;
    Directory *resultDir = &(dir->second);
    
    if(addMore){
      pair<string, Directory*> added = dir->second.addDirectory(nextDir);
      resultName = resultName + added.first;
      resultDir =  added.second;
    }
    
    return pair<string, Directory*>(resultName, resultDir);
  }
  
  string &Directory::stripSlash(string &path){
    string::iterator c = path.begin();
    
    if(strncmp(&(*c), "/", 1) == 0){
      path.erase(c);
    }
    
    return path;
  }
  
  void Directory::addHistogram(const string &name, TH1D *histo, bool divideBinWidths){
   
    if(m_histos[name]!=0) throw DuplicateHistogramException(this->name(), name);
    
    m_histos[name] = histo;
    
    if(!divideBinWidths){
      m_noDivideBinWidths.insert(histo);
    }
    
    return;
  }
  
  void Directory::addProfile(const string &name, TProfile *prof){
    if(m_profile1Ds[name]!=0) throw DuplicateHistogramException(this->name(), name);
    
    m_profile1Ds[name] = prof;
    
    return;
  }
  
  void Directory::addHistogram(const string &name, TH2D *histo, bool divideBinWidths){
    
    if(m_2dHistos[name] !=0) throw DuplicateHistogramException(this->name(), name);

    m_2dHistos[name] = histo;
    
    if(!divideBinWidths){
      m_noDivideBinWidths.insert(histo);
    }
    
    return;
  }
  
  void Directory::addHistogram(const string &name, TH3D *histo, bool divideBinWidths){
    if(m_3dHistos[name] != 0) throw DuplicateHistogramException(this->name(), name);
    
    m_3dHistos[name] = histo;
    
    if(!divideBinWidths){
      m_noDivideBinWidths.insert(histo);
    }
    
    return;
  }
  
  TH1D *Directory::histogram(const string &name){
    return m_histos[name];
  }
  
  TProfile *Directory::profile(const string &name){
    return m_profile1Ds[name];
  }
  
  TH2D *Directory::histogram2D(const string &name){
    return m_2dHistos[name];
  }
  
  TH3D *Directory::histogram3D(const string &name){
    return m_3dHistos[name];
  }
  
  bool Directory::doDivideBinWidths(const TH1 *histo)const{
    
    if(m_noDivideBinWidths.count(histo) != 0) return false;
      
    return true;
  }
  
  bool Directory::write(TDirectory *dir) const{

    if(dir==0) return false;
    
    dir->cd();
        
    for(boost::unordered_map<string, TH1D*>::const_iterator histo = m_histos.begin();
        histo != m_histos.end(); ++histo){
      histo->second->Write();
    }
    
    for(boost::unordered_map<string, TProfile*>::const_iterator prof = m_profile1Ds.begin();
        prof != m_profile1Ds.end(); ++prof){
      prof->second->Write();
    }
    
    for(boost::unordered_map<string, TH2D*>::const_iterator histo = m_2dHistos.begin();
        histo != m_2dHistos.end(); ++histo){
      histo->second->Write();
    }
    
    for(boost::unordered_map<string, TH3D*>::const_iterator histo = m_3dHistos.begin();
        histo != m_3dHistos.end(); ++histo){
      histo->second->Write();
    }
    
    for(map<string, Directory>::const_iterator dirIt = m_dirs.begin();
        dirIt != m_dirs.end(); ++dirIt){
      TDirectory *newDir = dir->mkdir(dirIt->second.m_name.c_str());
      dirIt->second.write(newDir);
    }
    
    return true;
  }
 
}