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();
  }
}