JFit.cc File Reference

General purpose fit program using ROOT. More...

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <algorithm>
#include "TROOT.h"
#include "TFile.h"
#include "TKey.h"
#include "TH1.h"
#include "TGraph.h"
#include "TProfile.h"
#include "TF1.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 using ROOT.

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 JFit.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 96 of file 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>         fixedValues;
  JRange_t                 X;
  JRange_t                 Y;
  char                     project;
  string                   option;
  int                      debug;

  try {

    JParser<> zap("General purpose fit program using ROOT.");

    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;\"");
    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['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"; }

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

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

          
        // execute fit

        JFit fit(*object, fcn, option);
        
        for_each(fit, JType<JTYPELIST<TH1, TGraph>::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();
      }
    }

    dir->Close();
  }

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