JPrintOscProb.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <string>
#include <map>

#include "JLang/JVectorize.hh"

#include "Jeep/JPrint.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

#include "JOscProb/JOscChannel.hh"
#include "JOscProb/JOscParameters.hh"
#include "JOscProb/JOscProbInterpolator2D.hh"


/**
 * \file
 * Example program to print oscillation probabilities interpolated in energy and zenith-angle.
 *
 * \author bjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  string                  inputFile;
  JOscParameters          parameters;
  map<double, double>     interpolationValues;
  JOscChannel             channel;
  int                     debug;

  try {
    
    JParser<> zap;

    zap['f'] = make_field(inputFile,           "oscillation probability table file");
    zap['x'] = make_field(interpolationValues, "Map of interpolation values in (energy [GeV], cos(theta)") = JPARSER::initialised();
    zap['@'] = make_field(parameters,          "oscillation parameters")                                   = JPARSER::initialised();
    zap['#'] = make_field(channel,             "Oscillation channel")                                      = JOscChannel(JOscChannel::MUON, JOscChannel::MUON, JOscChannel::PARTICLE);
    zap['d'] = make_field(debug)                                                                           = 3;

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

  const JOscProbTable6D<> g6(inputFile.c_str());
  const JOscProbTable2D<> g2(g6, parameters);
  
  const JOscProbInterpolator2D<> interpolator(&g2);

  NOTICE(RIGHT(15) << "Energy [GeV]" <<
         RIGHT(15) << "cos(theta)"   <<
         RIGHT(15) << "P"            << endl);
    
  for (map<double, double>::const_iterator i = interpolationValues.cbegin(); i != interpolationValues.cend(); ++i) {

    const double E       = i->first;
    const double costh   = i->second;
    const double oscProb = interpolator(channel, log10(E), costh);
      
    NOTICE(FIXED(15,3) << E       <<
           FIXED(15,3) << costh   <<
           FIXED(15,3) << oscProb << '\r'); DEBUG(endl);
  }
}