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 "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;
  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['d'] = make_field(debug)        = 1;

    zap(argc, argv);
  }
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }


  const double error = 0.2;

  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; ) {

      AddPoint(&g1, x, y, 0.0, error);

      if (x > xmax) {
        xmax = x;
        ymax = y;
      }
    }

    in.close();
  }

  DEBUG("Number of points " << g1.GetN() << endl);


  TF1 f1("f1", "[0]*(x + [1]*log(1.0 - [2]*x))");

  f1.SetParameter(0, ymax / xmax);
  f1.SetParameter(1, 0.25 * xmax);
  f1.SetParameter(2, 0.50 / xmax);
  
  g1.Fit(&f1);
  g1.SetName("g1");

  cout << inputFile << endl;
  
  cout << "T_x = " << FIXED(5,2) << f1.GetParameter(0) << " [s^2/m^2] * v_x^2" << endl;

  cout << "b [m]    " << FIXED(7,3) << f1.GetParameter(1) << endl;
  cout << "a [m^-1] " << FIXED(7,5) << f1.GetParameter(2) << endl;

  TFile out(outputFile.c_str(), "recreate");

  out << g1;

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