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 "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;
}