JFit.cc File Reference

General purpose fit program for 1D 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 "TH1.h"
#include "TGraph.h"
#include "TProfile.h"
#include "TF1.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 1D 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(x/[1])+[2]"
     -@ "p0 = GetMaximum; p1 = 2*GetRMS"
     -= "p2 = 0"

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

Author

mdejong

Definition in file software/JGizmo/JFit.cc.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Definition at line 100 of file software/JGizmo/JFit.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;
  char                     project;
  string                   option;
  int                      debug;
 
  try {
 
    JParser<> zap("General purpose fit program for 1D 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,           "ordinate range")                          = JRange_t();
    zap['P'] = make_field(project,     "projection")                              = '\0', 'x', 'X', 'y', 'Y';
    zap['O'] = make_field(option,      "Fit option")                              = "";
    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"; }
 
  bool px = (project == 'x' || project == 'X');  // projection on x-axis
  bool py = (project == 'y' || project == 'Y');  // projection on y-axis
 
  if (py) {
    swap(Y, X);   // Y becomes x-axis range and X becomes y-axis range
  }
 
 
  TFile out(outputFile.c_str(), "recreate");
 
 
  TF1* fcn = new TF1("user", formula.c_str());
 
  fcn->SetNpx(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();
      
        try {
 
          TProfile& h1 = dynamic_cast<TProfile&>(*object);
 
          object = h1.ProjectionX();
        }
        catch(exception&) {}
 
        try {
            
          TH2& h2 = dynamic_cast<TH2&>(*object);
 
          if        (px) {
 
            object = h2.ProjectionX("_px", 
                                    h2.GetYaxis()->FindBin(Y.getLowerLimit()),
                                    h2.GetYaxis()->FindBin(Y.getUpperLimit()));
            
          } else if (py) {
 
            object = h2.ProjectionY("_py", 
                                    h2.GetXaxis()->FindBin(Y.getLowerLimit()),
                                    h2.GetXaxis()->FindBin(Y.getUpperLimit()));
 
          } else {
 
            ERROR("For 2D histograms, use option option -P for projections." << endl);
 
            continue;
          }
        }           
        catch(exception&) {}
 
 
        // 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 j = 0; j != fcn->GetNpar(); ++j) {
          DEBUG(left << setw(12) << fcn->GetParName  (j) << ' ' << 
                SCIENTIFIC(12,5) << fcn->GetParameter(j) << endl);
        }
 
        Double_t xmin = numeric_limits<Double_t>::max();
        Double_t xmax = numeric_limits<Double_t>::lowest();
 
        {
          TH1* h1 = dynamic_cast<TH1*>(object);
 
          if (h1 != NULL) {
            xmin = min(xmin, h1->GetXaxis()->GetXmin());
            xmax = max(xmax, h1->GetXaxis()->GetXmax());
          }
        }
        
        {
          TGraph* g1 = dynamic_cast<TGraph*>(object);
 
          if (g1 != NULL) {
            for (Int_t i = 0; i != g1->GetN(); ++i) {
              xmin = min(xmin, g1->GetX()[i]);
              xmax = max(xmax, g1->GetX()[i]);
            }
          }
        }
 
        if (X != JRange_t()) {
          xmin = X.getLowerLimit();
          xmax = X.getUpperLimit();
        }
 
        fcn->SetRange(xmin, xmax);
 
          
        // execute fit
 
        const chrono::steady_clock::time_point t0 = chrono::steady_clock::now();
 
        JFit fit(*object, fcn, option);
        
        for_each(fit, JType<JTYPELIST<TH1, TGraph>::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();
      }
    }
 
    dir->Close();
  }
 
  out.Write();
  out.Close();
}