JOpera3D.cc

Go to the documentation of this file.

 
#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
 
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH3.h"
#include "TF3.h"
#include "TProfile3D.h"
#include "TString.h"
#include "TRegexp.h"
 
#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 * Auxiliary program for operations on 3D histograms.\n
 * The option <tt>-f</tt> corresponds to <tt><file name>:<object name></tt>.
 *
 * Possible operations (option <tt>-u <operation></tt>:
 *     - Add\n
 *       ROOT addition of histograms;
 *     - add\n
 *       Addition of histograms by look up of corresponding bin in second histogram;
 *     - Subtract\n
 *       ROOT subtraction of histograms;
 *     - subtract\n
 *       Subtraction of histograms by look up of corresponding bin in second histogram;
 *     - Multiply\n
 *       ROOT multiplication of histograms;
 *     - multiply\n
 *       Multiplication of histograms by look up of corresponding bin in second histogram;
 *     - Divide\n
 *       ROOT division of histograms;
 *     - divide\n
 *       Division of histograms by look up of corresponding bin in second histogram;
 *     - efficiency\n
 *       Like "Divide", but with correction for errors;
 *     - stdev\n
 *       Standard deviation;
 *     - Replace\n
 *       Replace histogram contents by function value;
 *
 * If only one histogram is specified, the operation applies to the histogram and the first associated function.\n
 * Only the quoted ROOT operations and the additional operation Replace are then allowed.\n
 * The option Replace will replace the histogram contents by the function value.
 * 
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  
  vector<JRootObjectID> inputFile;
  string                outputFile;
  string                opera;
  string                name;
  int                   debug;
 
  try {
 
    JParser<> zap("Auxiliary program for histogram operations.");
 
    zap['f'] = make_field(inputFile,  "<input file>:<object name>");
    zap['u'] = make_field(opera)  = 
      JOpera::Add(), 
      JOpera::add(), 
      JOpera::Subtract(), 
      JOpera::subtract(), 
      JOpera::Multiply(), 
      JOpera::multiply(),  
      JOpera::Divide(), 
      JOpera::divide(), 
      JOpera::efficiency(), 
      JOpera::stdev(),
      JOpera::Replace();
    zap['o'] = make_field(outputFile, "ROOT file with histogram(s)")    = "opera.root";
    zap['O'] = make_field(name,       "histogram name:"
                          << "\n\t\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
                          << "\n\t\"" << JOpera::SAME_AS_INPUT()     << "\" -> same as input; else"
                          << "\n\t as specified")                       = JOpera::SAME_AS_OPERATION();
    zap['d'] = make_field(debug)      = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  vector<TObject*> listOfObjects;
 
  for (vector<JRootObjectID>::const_iterator input = inputFile.begin(); input != inputFile.end(); ++input) {
 
    DEBUG("Input: " << *input << endl);
 
    TDirectory* dir = getDirectory(*input);
 
    if (dir == NULL) {
      ERROR("File: " << input->getFullFilename() << " not opened." << endl);
      continue;
    }
 
    const TRegexp regexp(input->getObjectName());
 
    TIter iter(dir->GetListOfKeys());
 
    for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
 
      const TString tag(key->GetName());
 
      DEBUG("Key: " << tag << " match = " << tag.Contains(regexp) << endl);
 
      // option match
 
      if (tag.Contains(regexp) && isTObject(key)) {
 
        TObject* p = key->ReadObj();
 
        if (dynamic_cast<TH3*>(p) == NULL) {
          FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
        }
 
        listOfObjects.push_back(p);
      }
    }
  }
 
 
  TH3* h3 = NULL;
 
  if        (listOfObjects.size() == 0) {
 
    FATAL("Number of histograms " << listOfObjects.size() << " = 0." << endl);
 
  } else if (listOfObjects.size() == 1 && (opera == JOpera::Add()      ||
                                           opera == JOpera::Subtract() ||
                                           opera == JOpera::Multiply() ||
                                           opera == JOpera::Divide()   ||
                                           opera == JOpera::Replace())) {
 
    // operation between histogram contents and first associated function
 
    TH3* h1 = dynamic_cast<TH3*>(listOfObjects[0]);
    TF2* f1 = dynamic_cast<TF2*>(h1->GetListOfFunctions()->First());
 
    if (f1 == NULL) {
      FATAL(h1->GetName() << " has no associated function." << endl);
    }
 
    NOTICE(h1->GetName() << ' ' << opera << ' ' << f1->GetName() << endl);
 
    if      (name == JOpera::SAME_AS_OPERATION())
      h3 = (TH3*) h1->Clone(opera.c_str());
    else if (name == JOpera::SAME_AS_INPUT())
      h3 = (TH3*) h1->Clone(h1->GetName());
    else
      h3 = (TH3*) h1->Clone(name.c_str());
 
    if        (opera == JOpera::Add()) {
 
      h3->Add     (f1, +1.0);
 
    } else if (opera == JOpera::Subtract()) {
 
      h3->Add     (f1, -1.0);
 
    } else if (opera == JOpera::Multiply()) {
 
      h3->Multiply(f1);
 
    } else if (opera == JOpera::Divide()) {
 
      h3->Divide  (f1);
    
    } else if (opera == JOpera::Replace()) {
 
      h3->Reset();
      h3->Fill(0.0, 0.0, 0.0);     // ensure one entry for drawing after reset
      h3->Add     (f1);
    }
 
  } else if (listOfObjects.size() == 2) {
 
    TH3* h1 = dynamic_cast<TH3*>(listOfObjects[0]);
    TH3* h2 = dynamic_cast<TH3*>(listOfObjects[1]);
 
    NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);
 
    if      (name == JOpera::SAME_AS_OPERATION())
      h3 = (TH3*) h1->Clone(opera.c_str());
    else if (name == JOpera::SAME_AS_INPUT())
      h3 = (TH3*) h1->Clone(h1->GetName());
    else
      h3 = (TH3*) h1->Clone(name.c_str());
 
    if        (opera == JOpera::Add()) {
 
      h3->Add     (h1, h2, +1.0, +1.0);
 
    } else if (opera == JOpera::Subtract()) {
        
      h3->Add     (h1, h2, +1.0, -1.0);
 
    } else if (opera == JOpera::Multiply()) {
        
      h3->Multiply(h1, h2, +1.0, +1.0);
 
    } else if (opera == JOpera::Divide()) {
 
      h3->Divide  (h1, h2, +1.0, +1.0);
 
    } else if (opera == JOpera::add()       ||
               opera == JOpera::subtract()  ||
               opera == JOpera::multiply()  ||
               opera == JOpera::divide()) {
 
      for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
        for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
          for (Int_t i3 = 1; i3 <= h1->GetNbinsZ(); ++i3) {
 
            const Double_t x  = h1->GetXaxis()->GetBinCenter(i1);
            const Double_t y  = h1->GetYaxis()->GetBinCenter(i2);
            const Double_t z  = h1->GetZaxis()->GetBinCenter(i3);
          
            const Int_t    j1 = h2->GetXaxis()->FindBin(x);
            const Int_t    j2 = h2->GetYaxis()->FindBin(y);
            const Int_t    j3 = h2->GetYaxis()->FindBin(z);
            
            const Double_t w1 = h1->GetBinContent(i1,i2,i3);
            const Double_t w2 = h2->GetBinContent(j1,j2,j3);
 
            Double_t w3 = 0.0;
            
            if        (opera == JOpera::add()) {
 
              w3 = w1 + w2;
 
            } else if (opera == JOpera::subtract()) {
 
              w3 = w1 - w2;
 
            } else if (opera == JOpera::multiply()) {
 
              w3 = w1 * w2;
 
            } else if (opera == JOpera::divide()) {
 
              if (w2 == 0.0) {
                ERROR("Bin " << h2->GetName() << "[" << j1 << "," << j2 << "] empty." << endl);
                continue;
              }
 
              w3 = w1 / w2;
            }
 
            h3->SetBinContent(i1, i2, i3, w3);
          }
        }
      }
 
    } else if (opera == JOpera::efficiency()  ||
               opera == JOpera::stdev()       ||
               opera == JOpera::sqrt()) {
 
      for (Int_t i1 = 1; i1 <= h1->GetNbinsX(); ++i1) {
        for (Int_t i2 = 1; i2 <= h1->GetNbinsY(); ++i2) {
          for (Int_t i3 = 1; i3 <= h1->GetNbinsZ(); ++i3) {       
 
            const Double_t y1 = h1->GetBinContent(i1,i2,i3);
            const Double_t y2 = h2->GetBinContent(i1,i2,i3);
 
            Double_t w1 = h1->GetBinError(i1,i2,i3);
            Double_t w2 = h2->GetBinError(i1,i2,i3);
            
            Double_t y3 = 0.0;
            Double_t w3 = 0.0;
            
            if        (opera == JOpera::efficiency()) {
              
              if (y2 == 0.0) {
                ERROR("Bin " << h2->GetName() << "[" << i1 << "," << i2 << "] empty." << endl);
                continue;
              }
              
              y3 = y1 / y2;
              
              if (y1 != 0.0 &&
                  y2 != 0.0) {
                
                w1 /= y1;
                w2 /= y2;
 
                w3  = y3 * fabs(y3 - 1.0) * sqrt(w1*w1 + w2*w2);
              }
              
            } else if (opera == JOpera::stdev()) {
 
              w3 = sqrt(w1*w1 + w2*w2);
 
              if (w3 != 0.0) {
                y3 = (y1 - y2) / w3;
                w3 = 0.0;
              }
 
            } else if (opera == JOpera::sqrt()) {
 
              y3 = (y1+y2) * (y1-y2);
 
              if (y3 < 0.0)
                y3 = -sqrt(-y3);
              else
                y3 = +sqrt(+y3);
 
              w3 = 0.0;
            }
        
            h3->SetBinContent(i1, i2, i3, y3);
            h3->SetBinError  (i1, i2, i3, w3);
          }
        }
      }
    }
 
  } else if (opera == JOpera::Add()      ||
             opera == JOpera::Multiply()) {
 
    for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {
 
      TH3* h1 = dynamic_cast<TH3*>(*i);
 
      NOTICE(h1->GetName() << ' ' << opera << endl);
 
      if (h3 == NULL) {
 
        if      (name == JOpera::SAME_AS_OPERATION())
          h3 = (TH3*) h1->Clone(opera.c_str());
        else if (name == JOpera::SAME_AS_INPUT())
          h3 = (TH3*) h1->Clone(h1->GetName());
        else
          h3 = (TH3*) h1->Clone(name.c_str());
 
      } else {
 
        if        (opera == JOpera::Add()) {
 
          h3->Add     (h3, h1, +1.0, +1.0);
 
        } else if (opera == JOpera::Multiply()) {
 
          h3->Multiply(h3, h1, +1.0, +1.0);
        }
      }
    }
 
  } else {
 
    FATAL("Incompatible number of histograms " << listOfObjects.size() << " with option " << opera << endl);
  }
 
  if (h3 != NULL) {
 
    TFile out(outputFile.c_str(), "recreate");
    
    h3->Write();
    
    out.Write();
    out.Close();
  }
}