JRootFit.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <limits>
 
#include "TROOT.h"
#include "TF1.h"
#include "TH1D.h"
#include "TFitResult.h"
 
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JMinimizer.hh"
 
#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Program to test ROOT fit.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  double   precision;
  int      debug;
 
  try {
 
    JParser<> zap("Program to test ROOT fit.");
 
    zap['e'] = make_field(precision)           = numeric_limits<double>::min();
    zap['d'] = make_field(debug)               = 3;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  string option = "";
 
  if (option.find('S') == string::npos) {
    option += 'S';
  }
 
  if (debug < debug_t && option.find('Q') == string::npos) {
    option += "Q";
  }
 
 
  TH1D h1("h1", NULL, 1, -0.5, +0.5);
 
  
  struct {
    double y;
    double z;
  } data[] = { 
    { 100.0, 10.0 }, 
    {   0.0, 10.0 }
  };
 
 
  for (int i = 0; i != sizeof(data)/sizeof(data[0]); ++i) {
 
    const double y = data[i].y;
    const double z = data[i].z;
 
    h1.SetBinContent(1, y);
    h1.SetBinError  (1, z);
 
    {
      NOTICE("ROOT fit " << FIXED(7,3) << y << " +/- " << FIXED(7,3) << z << endl);
 
      TF1 f1("f1", "[0]");
 
      double ymin, ymax;
 
      f1.SetParLimits(0, 0.0,  0.0);
      f1.GetParLimits(0, ymin, ymax);
 
      const TFitResultPtr result = h1.Fit(&f1, option.c_str());
 
      ASSERT(result.Get() != NULL, "TFitResultPtr");
 
      ERROR("Result: " 
            << FIXED(7,3) << f1.GetParameter(0) << " +/- " 
            << FIXED(7,3) << f1.GetParError (0) << " "
            << "[" << FIXED(7,3) << ymin << "," << FIXED(7,3) << ymax << "]" << " "
            << (result->IsParameterBound(0) ? " *** bound *** " : "")  
            << (result->IsParameterFixed(0) ? " *** fixed *** " : "")  
            << endl);
    }
    {
      NOTICE("Normal fit " << FIXED(7,3) << y << " +/- " << FIXED(7,3) << z << endl);
    
      TF1 f1("f1", "[0]");
 
      ASSERT(setParameter(f1, JFitParameter_t(0, y)) == true);
      ASSERT(setParameter(f1, JFitParameter_t(1, y)) == false);
 
      const TFitResultPtr result = h1.Fit(&f1, option.c_str());
 
      ASSERT(result.Get() != NULL, "TFitResultPtr");
 
      DEBUG("Result: " 
            << FIXED(7,3) << f1.GetParameter(0) << " +/- " 
            << FIXED(7,3) << f1.GetParError (0) << " "
            << (result->IsParameterBound(0) ? " *** bound *** " : "")  
            << (result->IsParameterFixed(0) ? " *** fixed *** " : "")  
            << endl);
 
      ASSERT(result->IsValid());
      ASSERT(!result->IsParameterBound(0));
      ASSERT(!result->IsParameterFixed(0));
      ASSERT(fabs(f1.GetParameter(0) - y) < f1.GetParError(0));
    }
    {
      NOTICE("Fixed parameters fit " << FIXED(7,3) << y << " +/- " << FIXED(7,3) << z << endl);
 
      TF1 f1("f1", "[0]");
 
      ASSERT(fixParameter(f1, JFitParameter_t(0, y)) == true);
      ASSERT(isParameterFixed(f1, 0) == true);
      ASSERT(isParameterFixed(f1, 1) == false);
 
      const JFitParameter_t parameter(0, 0.0);
 
      ASSERT(fixParameter(f1, parameter));
      ASSERT(fixParameter(f1, parameter));
 
      const TFitResultPtr result = h1.Fit(&f1, option.c_str());
 
      ASSERT(result.Get() != NULL, "TFitResultPtr");
 
      DEBUG("Result: " 
            << FIXED(7,3) << f1.GetParameter(0) << " +/- " 
            << FIXED(7,3) << f1.GetParError (0) << " "
            << (result->IsParameterBound(0) ? " *** bound *** " : "")  
            << (result->IsParameterFixed(0) ? " *** fixed *** " : "")  
            << endl);
 
      ASSERT(result->IsValid());
      ASSERT(!result->IsParameterBound(0));
      ASSERT(result->IsParameterFixed(0));
      ASSERT(fabs(f1.GetParameter(0) - parameter.value) < precision);
      ASSERT(fabs(f1.GetParError (0))                   < precision);
    }
    {
      NOTICE("Fixed parameters fit " << FIXED(7,3) << y << " +/- " << FIXED(7,3) << z << endl);
 
      TF1 f1("f1", "[0]");
 
      ASSERT(fixParameter(f1, JFitParameter_t(0, y)) == true);
      ASSERT(isParameterFixed(f1, 0) == true);
      ASSERT(isParameterFixed(f1, 1) == false);
 
      const JFitParameter_t parameter(0, 0.0);
 
      ASSERT(fixParameter(f1, parameter));
      ASSERT(fixParameter(f1, parameter));
 
      const TFitResultPtr result = h1.Fit(&f1, option.c_str());
 
      ASSERT(result.Get() != NULL, "TFitResultPtr");
 
      DEBUG("Result: " 
            << FIXED(7,3) << f1.GetParameter(0) << " +/- " 
            << FIXED(7,3) << f1.GetParError (0) << " "
            << (result->IsParameterBound(0) ? " *** bound *** " : "")  
            << (result->IsParameterFixed(0) ? " *** fixed *** " : "")  
            << endl);
 
      ASSERT(result->IsValid());
      ASSERT(!result->IsParameterBound(0));
      ASSERT(result->IsParameterFixed(0));
      ASSERT(fabs(f1.GetParameter(0) - parameter.value) < precision);
      ASSERT(fabs(f1.GetParError (0))                   < precision);
    }
    {
      NOTICE("Released parameters fit " << FIXED(7,3) << y << " +/- " << FIXED(7,3) << z << endl);
 
      TF1 f1("f1", "[0]");
 
      ASSERT(fixParameter(f1, JFitParameter_t(0, y)) == true);
      ASSERT(isParameterFixed(f1, 0) == true);
 
      ASSERT(releaseParameter(f1, 0) == true);
      ASSERT(releaseParameter(f1, 1) == false);
 
      ASSERT(!isParameterFixed(f1, 0));
 
      const TFitResultPtr result = h1.Fit(&f1, option.c_str());
 
      ASSERT(result.Get() != NULL, "TFitResultPtr");
 
      DEBUG("Result: " 
            << FIXED(7,3) << f1.GetParameter(0) << " +/- " 
            << FIXED(7,3) << f1.GetParError (0) << " "
            << (result->IsParameterBound(0) ? " *** bound *** " : "")  
            << (result->IsParameterFixed(0) ? " *** fixed *** " : "")  
            << endl);
 
      ASSERT(result->IsValid());
      ASSERT(!result->IsParameterBound(0));
      ASSERT(!result->IsParameterFixed(0));
      ASSERT(fabs(f1.GetParameter(0) - y) < f1.GetParError(0));
    }
  }
 
  return 0;
}