JOpera2D.cc File Reference

Auxiliary program for operations on 2D histograms. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TF2.h"
#include "TProfile2D.h"
#include "TString.h"
#include "TRegexp.h"
#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Auxiliary program for operations on 2D histograms.

The option -f corresponds to <file name>:<object name>.

Possible operations (option -u <operation>:

• Add
  ROOT addition of histograms;
• add
  Addition of histograms by look up of corresponding bin in second histogram;
• Subtract
  ROOT subtraction of histograms;
• subtract
  Subtraction of histograms by look up of corresponding bin in second histogram;
• Multiply
  ROOT multiplication of histograms;
• multiply
  Multiplication of histograms by look up of corresponding bin in second histogram;
• Divide
  ROOT division of histograms;
• divide
  Division of histograms by look up of corresponding bin in second histogram;
• efficiency
  Like "Divide", but with correction for errors;
• stdev
  Standard deviation;

If only one histogram is specified, the operation applies to the histogram and the first associated function.
Only the quoted ROOT operations and the additional operation Replace are then allowed.
The option Replace will replace the histogram contents by the function value.

Author

mdejong

Definition in file JOpera2D.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 57 of file JOpera2D.cc.

{
  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,       MAKE_STRING("histogram name;"
                                                  << "\n\"" << JOpera::SAME_AS_OPERATION() << "\" -> same as operation; or"
                                                  << "\n\"" << JOpera::SAME_AS_INPUT()     << "\" -> same as input; else"
                                                  << "\nas 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)) {

        TObject* p = key->ReadObj();

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

        listOfObjects.push_back(p);
      }
    }
  }


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

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

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

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

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

          const Double_t x  = h1->GetXaxis()->GetBinCenter(i1);
          const Double_t y  = h1->GetYaxis()->GetBinCenter(i2);
          
          const Int_t    j1 = h2->GetXaxis()->FindBin(x);
          const Int_t    j2 = h2->GetYaxis()->FindBin(y);
            
          const Double_t w1 = h1->GetBinContent(i1,i2);
          const Double_t w2 = h2->GetBinContent(j1,j2);

          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, 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) {

          const Double_t y1 = h1->GetBinContent(i1,i2);
          const Double_t y2 = h2->GetBinContent(i1,i2);

          Double_t w1 = h1->GetBinError(i1,i2);
          Double_t w2 = h2->GetBinError(i1,i2);
            
          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, y3);
          h3->SetBinError  (i1, i2, w3);
        }
      }
    }

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

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

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

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

      if (h3 == NULL) {

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