JMechanics.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
 
#include "TROOT.h"
#include "TFile.h"
#include "TGraphErrors.h"
#include "TF1.h"
 
#include "JLang/JWhiteSpacesFacet.hh"
#include "JROOT/JRootToolkit.hh"
#include "JROOT/JMinimizer.hh"
 
#include "JMath/JMathlib.hh"
#include "JROOT/JRootfit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Auxiliary application to fit mechanical constants.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
 
  string    inputFile;
  string    outputFile;
  double    error;
  bool      nobuoy;
  bool      old;
  int       debug;
 
  try {
 
    JParser<> zap("Auxiliary application to fit mechanical constants.");
    
    zap['f'] = make_field(inputFile,      "input file; see https://git.km3net.de/calibration/input_tables");
    zap['o'] = make_field(outputFile)   = "mechanics.root";
    zap['e'] = make_field(error,          "uncertainty drag value")        = 0.2;
    zap['B'] = make_field(nobuoy,         "fit without buoy");
    zap['O'] = make_field(old,            "ROOT fit based on minuit");
    zap['d'] = make_field(debug)        = 1;
 
    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
 
  TGraphErrors g1;
 
  double xmax = 0.0;
  double ymax = 0.0;
 
  {
    ifstream in(inputFile.c_str());
 
    string buffer;
 
    getline(in, buffer);
  
    DEBUG("Header " << buffer << endl);
 
    const locale loc(in.getloc(), new JWhiteSpacesFacet(in.getloc(), ", \t\r\n"));
 
    in.imbue(loc);
  
    for (double x, y; in >> x >> y; ) {
 
      DEBUG(FIXED(7,3) << x << ' '<< FIXED(9,4) << y << endl);
 
      AddPoint(&g1, x, y, 0.0, error);
 
      if (x > xmax) {
        xmax = x;
        ymax = y;
      }
    }
 
    in.close();
  }
 
 
  // start -> final values
 
  double Tx  =  ymax / xmax;
  double a   = (nobuoy ? 0.0 : 0.5  / xmax);
  double b   = (nobuoy ? 0.0 : 0.25 * xmax);
 
 
  // fit
 
  if (!old) {
 
    auto f1 = P0<0>(Tx) * (X + P0<1>(b) * Log(1.0 - P0<2>(a) * X));
 
    typedef decltype(f1)  function_type;
 
    JRootfit<function_type>::debug               =  debug;
 
    const auto result = (nobuoy ? Fit(&g1, f1, {1, 2}) : Fit(&g1, f1));
 
    result.print(cout);
 
    setFormat<JMatrixND>(JFormat_t(12, 3, std::ios::scientific | std::ios::showpos));
 
    cout << "Error matrix:" << endl << result.V << endl;
  
    Tx = result.getValue(0);
    b  = result.getValue(1);
    a  = result.getValue(2);
 
  } else {
 
    TF1 f1("f1", "[0]*(x + [1]*log(1.0 - [2]*x))");
 
    f1.SetParameter(0, Tx);
    if (nobuoy) {
      f1.FixParameter(1, b);
      f1.FixParameter(2, a);
    } else {
      f1.SetParameter(1, b);
      f1.SetParameter(2, a);
    }
 
    f1.SetParError(0, 0.05);
    f1.SetParError(1, 1.0e-6);
    f1.SetParError(2, 1.0e-1);     // anomalous step size works?
    //f1.SetParError(2, 1.0e-7);
 
    g1.Fit(&f1);
 
    Tx = f1.GetParameter(0);
    b  = f1.GetParameter(1);
    a  = f1.GetParameter(2);
  }
 
 
  cout << inputFile << endl;
 
  cout << "T_x = " << FIXED(5,2) << Tx << " [s^2/m^2] * v_x^2" << endl;
 
  cout << "b [m]    " << FIXED(7,3) << b << endl;
  cout << "a [m^-1] " << FIXED(7,5) << a << endl;
 
 
  TFile out(outputFile.c_str(), "recreate");
 
  g1.SetName("g1");
 
  out << g1;
 
  out.Write();
  out.Close();
}