JOpera1D.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 "TH1.h"
#include "TF1.h"
#include "TProfile.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 1D histograms.
 * 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::sqrt(),
      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();
        TProfile* q = dynamic_cast<TProfile*>(p);

        if (q != NULL) {
          p = q->ProjectionX();
        }

        if (dynamic_cast<TH1*>(p) == NULL) {
          FATAL("Object " << p->GetName() << " not compatible with histogram operations." << endl);
        }
        
        listOfObjects.push_back(p);
      }
    }
  }


  TH1* 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

    TH1* h1 = dynamic_cast<TH1*>(listOfObjects[0]);
    TF1* f1 = (TF1*) 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 = (TH1*) h1->Clone(opera.c_str());
    else if (name == JOpera::SAME_AS_INPUT())
      h3 = (TH1*) h1->Clone(h1->GetName());
    else
      h3 = (TH1*) 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);          // ensure one entry for drawing after reset
      h3->Add     (f1);
    }
    
  } else if (listOfObjects.size() == 2) {

    // operation between two histograms

    TH1* h1 = dynamic_cast<TH1*>(listOfObjects[0]);
    TH1* h2 = dynamic_cast<TH1*>(listOfObjects[1]);

    NOTICE(h1->GetName() << ' ' << opera << ' ' << h2->GetName() << endl);

    if      (name == JOpera::SAME_AS_OPERATION())
      h3 = (TH1*) h1->Clone(opera.c_str());
    else if (name == JOpera::SAME_AS_INPUT())
      h3 = (TH1*) h1->Clone(h1->GetName());
    else
      h3 = (TH1*) 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 i = 1; i <= h1->GetNbinsX(); ++i) {

        const Double_t x  = h1->GetBinCenter (i);
        const Double_t w1 = h1->GetBinContent(i);
          
        const Int_t    j  = h2->FindBin(x);
        const Double_t w2 = h2->GetBinContent(j);

        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() << "[" << j << "] empty." << endl);
            continue;
          }
              
          w3 = w1 / w2;
        }

        h3->SetBinContent(i, w3);
      }

    } else if (opera == JOpera::efficiency() ||
               opera == JOpera::stdev()      ||
               opera == JOpera::sqrt()) {

      for (Int_t i = 1; i <= h1->GetNbinsX(); ++i) {

        const Double_t y1 = h1->GetBinContent(i);
        const Double_t y2 = h2->GetBinContent(i);

        Double_t w1 = h1->GetBinError(i);
        Double_t w2 = h2->GetBinError(i);

        Double_t y3 = 0.0;
        Double_t w3 = 0.0;

        if        (opera == JOpera::efficiency()) {
          
          if (y2 == 0.0) {
            ERROR("Bin " << h2->GetName() << "[" << i << "] 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(i, y3);
        h3->SetBinError  (i, w3);
      }
    }

  } else if (opera == JOpera::Add()      ||
             opera == JOpera::Multiply()) {

    // operation between more than two histograms

    for (vector<TObject*>::iterator i = listOfObjects.begin(); i != listOfObjects.end(); ++i) {

      TH1* h1 = dynamic_cast<TH1*>(*i);

      NOTICE(h1->GetName() << ' ' << opera << endl);

      if (h3 == NULL) {

        if      (name == JOpera::SAME_AS_OPERATION())
          h3 = (TH1*) h1->Clone(opera.c_str());
        else if (name == JOpera::SAME_AS_INPUT())
          h3 = (TH1*) h1->Clone(h1->GetName());
        else
          h3 = (TH1*) 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();
  }
}