JFit2D.cc File Reference

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

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include <chrono>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH2.h"
#include "TGraph2D.h"
#include "TGraph2DErrors.h"
#include "TF2.h"
#include "TFitResult.h"
#include "TString.h"
#include "TRegexp.h"
#include "JROOT/JMinimizer.hh"
#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 102 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>         startErrors;
  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['E'] = make_field(startErrors, "start errors, e.g: \"p0 = 0.01 * GetMaximum;\"")            = JPARSER::initialised();
    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 function(s) to ROOT file " << "(\"<name>" << _F2 << "\")");
    zap['M'] = make_field(minimizer,   "ROOT minimizer type and algorithm [debug]")                 = JMinimizer();
    zap['d'] = make_field(debug)        = 1;

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


  if (option.find('O') == string::npos) { option += "O"; }
  if (option.find("R") == string::npos) { option += "R"; }
  if (option.find("S") == string::npos) { option += "S"; }
  //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(1000);
  fcn->SetNpy(1000);

  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 (Int_t i = 0; i != fcn->GetNpar(); ++i) {
            fcn->SetParError(i, 0.0);
          }

          for (vector<JToken_t>::const_iterator i = startErrors.begin(); i != startErrors.end(); ++i) {
            fcn->SetParError (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);

          
        // execute fit

        const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();

        JFit fit(*object, fcn, option);

        for_each(fit, JType<JTYPELIST<TH2, TGraph2D>::typelist>());
  
        const chrono::steady_clock::time_point t1 = chrono::steady_clock::now();

        if (fit.result != -1) {

          // output fit results
            
          NOTICE("Fit values  " << object->GetName() << endl);
          NOTICE("Fit formula " << formula           << endl);
          NOTICE("chi2/NDF "    << FIXED(7,3) << fit.result->Chi2() << '/' << fit.result->Ndf() << ' ' << (fit.result->IsValid() ? "" : "failed") << endl);
          NOTICE("Number of calls " << fit.result->NCalls() << endl);
          NOTICE("Elapsed time [us] " << setw(8) << chrono::duration_cast<chrono::microseconds>(t1 - t0).count() << 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();
}