JFit2D.cc File Reference

General purpose fit program for 2D ROOT objects. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TGraph2D.h"
#include "TGraph2DErrors.h"
#include "TF2.h"
#include "TString.h"
#include "TRegexp.h"
#include "JLang/JType.hh"
#include "JLang/JTypeList.hh"
#include "JLang/JToken.hh"
#include "JTools/JRange.hh"
#include "JGizmo/JRootObjectID.hh"
#include "JGizmo/JGizmoToolkit.hh"
#include "Jeep/JPrint.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

General purpose fit program for 2D ROOT objects.

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

The expressions for the fit function, start and fixed values should comply with ROOT class TFormula.

In the expressions of the start and fixed values, names of member methods of corresponding class of the fit object may appear, such as TH1::GetMaximum, TH1::GetRMS, etc., e.g:

     -F "[0]*exp(-0.5*x*x/([1]*[1])*exp(-0.5*y*y/([2]*[2])"
     -@ "p0 = GetMaximum; p1 = 2*GetRMS(1); .."
     -= "p2 = 1"

The result of the formulas for the start and fixed values will be evaluated for each histogram separately.

Author

mdejong

Definition in file JFit2D.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 96 of file JFit2D.cc.

{
  using namespace std;
  using namespace JPP;

  typedef JToken<';'>      JToken_t;
  typedef JRange<double>   JRange_t;

  vector<JRootObjectID>    inputFile;
  string                   outputFile;
  string                   formula;
  vector<JToken_t>         startValues;
  vector<JToken_t>         fixedValues;
  vector<JToken_t>         limitValues;
  JRange_t                 X;
  JRange_t                 Y;
  string                   option;
  bool                     writeFits;
  int                      debug;

  try {

    JParser<> zap("General purpose fit program for 2D ROOT objects.");

    zap['f'] = make_field(inputFile,   "<input file>:<object name>");
    zap['o'] = make_field(outputFile,  "ROOT file with fit results")         = "fit.root";
    zap['F'] = make_field(formula,     "fit formula, e.g: \"[0]+[1]*x\"");
    zap['@'] = make_field(startValues, "start values, e.g: \"p0 = GetMaximum;\"");
    zap['='] = make_field(fixedValues, "fixed values, e.g: \"p0 = GetMaximum;\"")                   = JPARSER::initialised();
    zap['R'] = make_field(limitValues, "limit values, e.g: \"p0 = <lower limit> <upper limit>;\"")  = JPARSER::initialised();
    zap['x'] = make_field(X,           "abscissa range")                     = JRange_t();
    zap['y'] = make_field(Y,           "abscissa range")                     = JRange_t();
    zap['O'] = make_field(option,      "Fit option")                         = "";
    zap['w'] = make_field(writeFits,   "write fit results to ROOT file");
    zap['d'] = make_field(debug)        = 1;
    
    zap['='] = JPARSER::initialised();

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  if (option.find('O') == string::npos) { option += "O"; }
  //if (option.find('N') == string::npos) { option += "N"; }
  if (debug ==  0 && option.find('Q') == string::npos) { option += "Q"; }


  TFile out(outputFile.c_str(), "recreate");


  TF2* fcn = new TF2("user", formula.c_str());

  fcn->SetNpx(10000);

  if (fcn->IsZombie()) {
    FATAL("Function: " << formula << " is zombie." << endl);
  }

  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* object = key->ReadObj();
      

        // set fit parameters

        try {

          for (vector<JToken_t>::const_iterator i = startValues.begin(); i != startValues.end(); ++i) {
            fcn->SetParameter(getParameter(*i), getValue(*i,object));
          }
          
          for (vector<JToken_t>::const_iterator i = fixedValues.begin(); i != fixedValues.end(); ++i) {
            fcn->FixParameter(getParameter(*i), getValue(*i,object));
          }

          for (vector<JToken_t>::const_iterator i = limitValues.begin(); i != limitValues.end(); ++i) {
            fcn->SetParLimits(getParameter(*i), getValue(*i,0), getValue(*i,1));
          }
        }
        catch(JLANG::JParseError& error) {
          FATAL(error << endl);
        }

        DEBUG("Start values " << object->GetName() << endl);

        for (int i = 0; i != fcn->GetNpar(); ++i) {
          DEBUG(left << setw(12) << fcn->GetParName  (i) << ' ' << 
                SCIENTIFIC(12,5) << fcn->GetParameter(i) << endl);
        }

        Double_t xmin = numeric_limits<Double_t>::max();
        Double_t xmax = numeric_limits<Double_t>::lowest();
        Double_t ymin = numeric_limits<Double_t>::max();
        Double_t ymax = numeric_limits<Double_t>::lowest();

        {
          TH2* h2 = dynamic_cast<TH2*>(object);

          if (h2 != NULL) {

            xmin = min(xmin, h2->GetXaxis()->GetXmin());
            xmax = max(xmax, h2->GetXaxis()->GetXmax());
            ymin = min(ymin, h2->GetYaxis()->GetXmin());
            ymax = max(ymax, h2->GetYaxis()->GetXmax());

            if (X != JRange_t()) { h2->GetXaxis()->SetRangeUser(X.getLowerLimit(), X.getUpperLimit()); }
            if (Y != JRange_t()) { h2->GetYaxis()->SetRangeUser(Y.getLowerLimit(), Y.getUpperLimit()); }
          }
        }
        
        {
          TGraph2D* g2 = dynamic_cast<TGraph2D*>(object);

          if (g2 != NULL) {
            for (Int_t i = 0; i != g2->GetN(); ++i) {
              xmin = min(xmin, g2->GetX()[i]);
              xmax = max(xmax, g2->GetX()[i]);
              ymin = min(ymin, g2->GetY()[i]);
              ymax = max(ymax, g2->GetY()[i]);
            }
          }
        }

        if (X != JRange_t()) {
          xmin = X.getLowerLimit();
          xmax = X.getUpperLimit();
        }

        if (Y != JRange_t()) {
          ymin = Y.getLowerLimit();
          ymax = Y.getUpperLimit();
        }

        fcn->SetRange(xmin, ymin, xmax, ymax);

        if (option.find("R") == string::npos) {
          option += "R";
        }

          
        // execute fit

        JFit fit(*object, fcn, option);

        for_each(fit, JType<JTYPELIST<TH2, TGraph2D>::typelist>());
  

        if (fit.result != -1) {

          // output fit results
            
          NOTICE("Fit values  " << object->GetName() << endl);
          NOTICE("Fit formula " << formula           << endl);
            
          for (int j = 0; j != fcn->GetNpar(); ++j) {
            NOTICE(left << setw(12) << fcn->GetParName  (j) << ' '     << 
                   SCIENTIFIC(12,5) << fcn->GetParameter(j) << " +/- " << 
                   SCIENTIFIC(12,5) << fcn->GetParError (j) << endl);
          }

        } else {

          WARNING("Object: not compatible with ROOT Fit." << endl);
        }

        out.cd();

        object->Write();
        fcn   ->Write();

        if (writeFits) {

          TObject* p = object->Clone(MAKE_CSTRING(object->GetName() << ".f2"));

          {
            TH2* h2 = dynamic_cast<TH2*>(p);

            if (h2 != NULL) {
              for (Int_t ix = 1; ix <= h2->GetXaxis()->GetNbins(); ++ix) {
                for (Int_t iy = 1; iy <= h2->GetYaxis()->GetNbins(); ++iy) {

                  const Double_t x = h2->GetXaxis()->GetBinCenter(ix);
                  const Double_t y = h2->GetYaxis()->GetBinCenter(iy);

                  h2->SetBinContent(ix, iy, fcn->Eval(x,y));
                  h2->SetBinError  (ix, iy, 0.0);
                }
              }
            }
          }

          {
            TGraph2D* g2 = dynamic_cast<TGraph2D*>(p);

            if (g2 != NULL) {
              for (Int_t i = 0; i != g2->GetN(); ++i) {

                const Double_t x = g2->GetX()[i];
                const Double_t y = g2->GetY()[i];

                g2->GetZ()[i] = fcn->Eval(x,y);
              }
            }
          }

          {
            TGraph2DErrors* g2 = dynamic_cast<TGraph2DErrors*>(p);

            if (g2 != NULL) {
              for (Int_t i = 0; i != g2->GetN(); ++i) {
                g2->SetPointError(i, 0.0, 0.0, 0.0);
              }
            }
          }
        }
      }
    }

    dir->Close();
  }

  out.Write();
  out.Close();
}