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