JWifi.cc File Reference

#include <string>
#include <iostream>
#include <iomanip>
#include <chrono>
#include <map>
#include <sstream>
#include <cmath>
#include <vector>
#include <mutex>
#include <thread>
#include <condition_variable>
#include <queue>
#include "TROOT.h"
#include "TFile.h"
#include "TH2D.h"
#include "JGizmo/JGizmoToolkit.hh"
#include "JPhysics/JPDF.hh"
#include "JPhysics/JPDFSupportkit.hh"
#include "JPhysics/JPhysicsSupportkit.hh"
#include "JPhysics/KM3NeT.hh"
#include "JMath/JConstants.hh"
#include "JTools/JRange.hh"
#include "JGeometry3D/JAngle3D.hh"
#include "JGeometry3D/JOmega3D.hh"
#include "JFit/JGradient.hh"
#include "JLang/JSTDObjectReader.hh"
#include "Jeep/JProperties.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)
	Application to fit parameters to stopping muon data.

Function Documentation

main()

int main	(	int	argc,
		char **	argv )

Application to fit parameters to stopping muon data.

Author

mdejong

Definition at line 705 of file JWifi.cc.

{
  using namespace std;
  using namespace JPP;
 
  typedef std::pair<double, double>    parameter_type;
  typedef map<string, parameter_type>  option_type;
 
  experiment_type    setup;
  JAbsorptionLength  absorptionLength;
  JScatteringLength  scatteringLength;
  double             E_GeV;
  double             angle_Deg;
  size_t             threads;
  option_type        option;
  bool               normalise;
  vector<double>     xs;
  int                debug;
 
  try {
 
    JProperties properties;
 
    properties.insert(gmake_property(absorptionLength));
    properties.insert(gmake_property(scatteringLength));
 
    JParser<> zap;
 
    zap['@'] = make_field(properties, endl
                          << "possible options absorptionLength: "      << get_keys(absorptionLength)      << endl
                          << "possible options scatteringLength: "      << get_keys(scatteringLength)      << endl)  = JPARSER::initialised();
    zap['f'] = make_field(setup,
                          "douplet of histograms corresponding to scattered and direct light, "
                          << "each of which defined by <file name>:<histogram name>");
    zap['E'] = make_field(E_GeV,        "muon energy [GeV]")                             =  10.0;
    zap['G'] = make_field(angle_Deg,    "grid angle for PMT directions [deg]")           =  25.0;
    zap['M'] = make_field(JWifi::bins,  "number of points in distance and time bin")     =  JPARSER::initialised();
    zap['N'] = make_field(threads,      "number threads")                                =   0;
    zap['O'] = make_field(option,       "fit options: \"parameter <start value> <step size>\"" << endl
                          << "possible parameters: "
                          << absorption_factor_t  << ", "
                          << absorption_offset_t  << ", "
                          << scattering_factor_t  << ", "
                          << qe_factor_t          << ", "
                          << mixed_factor_a_t     << ", "
                          << mixed_factor_b_t     << ", "
                          << mixed_offset_a_t     << ", "
                          << mixed_offset_b_t);
    zap['n'] = make_field(normalise,  "normalise gradient");
    zap['x'] = make_field(xs,         "scan step sizes")                                 =  JPARSER::initialised();
    zap['d'] = make_field(debug)               = 2;
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
 
  const TH2* ha = dynamic_cast<const TH2*>(getObject(setup.Ha));
  const TH2* hb = dynamic_cast<const TH2*>(getObject(setup.Hb));
 
  if (ha == NULL) { FATAL("Missing histogram: " << setup.Ha << endl); }
  if (hb == NULL) { FATAL("Missing histogram: " << setup.Hb << endl); }
 
  const JWifi wifi(E_GeV, ha, hb, angle_Deg, threads);
 
  // fit
 
  JGradient fit(FIT::number_of_iterations, FIT::number_of_extra_steps, FIT::epsilon, 3);
 
  fit.normalise = normalise;
 
  auto fp = [&option, &fit](const string& key, double& value, const range_type& range) {
    if (option.count(key)) {
      fit.push_back(JModifier_t(key, new JEditor(value = option[key].first, range), option[key].second));
    }
  };
 
  fp(absorption_factor_t, JAbsorptionLength::get_factor(),   range_type(0.1, 1.0e1));
  fp(absorption_offset_t, JAbsorptionLength::get_offset(),   range_type(1.0, 1.0e6));
  fp(scattering_factor_t, JScatteringLength::get_factor(),   range_type(0.1, 1.0e1));
  fp(qe_factor_t,         JQE              ::get_factor(),   range_type(0.1, 1.0e1));
  fp(mixed_factor_a_t,    JMixed           ::get_factor_a(), range_type(0.1, 1.0e1));
  fp(mixed_factor_b_t,    JMixed           ::get_factor_b(), range_type(0.1, 1.0e1));
  fp(mixed_offset_a_t,    JMixed           ::get_offset_a(), range_type(1.0, 1.0e6));
  fp(mixed_offset_b_t,    JMixed           ::get_offset_b(), range_type(1.0, 1.0e6));
 
  const chrono::high_resolution_clock::time_point t0 = chrono::high_resolution_clock::now();
 
  const double chi2 = fit(wifi);
 
  const chrono::high_resolution_clock::time_point t1 = chrono::high_resolution_clock::now();
 
  cout << "chi2/NDF " << FIXED(9,5) << chi2 << "/" << setw(4) << (wifi.getN() - fit.size()) << endl;
 
  cout << "Elapsed time: " << setw(6) << (t1 - t0) / chrono::seconds(1) << " [s]" << endl;
 
  for (const auto& i : fit) {
    cout << left << setw(36) << i.name << ' ' << right << FIXED(13,3) << dynamic_cast<const JEditor*>(i.get())->getValue() << endl;
  }
 
  auto get_offset = [](const double x1, const double x2) { const double y = 1.0/x1 + 1.0/x2; return 1.0/y; };
  auto get_factor = [](const double x1, const double x2) { return x1 * x2; };
 
  cout << "Additional absorption length: " << FIXED(13,5) << get_offset(JAbsorptionLength::get_offset(), JMixed::getOffsetAbsorptionLength()) << " [m]" << endl;
  cout << "Scaling of absorption length: " << FIXED(13,5) << get_factor(JAbsorptionLength::get_factor(), JMixed::getFactorAbsorptionLength())           << endl;
  cout << "Scaling of scattering length: " << FIXED(13,5) << JScatteringLength::get_factor()                                                            << endl;
  cout << "Scaling of QE:                " << FIXED(13,5) << get_factor(JQE              ::get_factor(), JMixed::getFactorQE())                         << endl;
 
  if (!xs.empty()) {
    fit(cout, wifi, xs);
  }
}