background_trajectory.cc

Go to the documentation of this file.

//background_trajectory
 
#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <sstream>
 
#include "TROOT.h"
#include "TFile.h"
#include "TGraphErrors.h"
 
#include "Jeep/JParser.hh"
#include "Jeep/JeepToolkit.hh"
#include "JTools/JRange.hh"
 
#include "JLang/JManip.hh"
#include "JAstronomy/JAstronomy.hh"
 
#include "slalib/slalib.h"
 
#include "flux/Flux.hh"
 
/**
 * \file
 * Auxiliary program to integrate atmospheric neutrino flux over apparent trajectory of a source.
 *
 * \author fbaas
 */
 
int main(int argc, char **argv){
  using namespace std;
  using namespace JPP;
 
  typedef JRange<double> Range;
 
  double       dec;
  double       RA;
  Range        T_start_int;
  Range        E_min_max;
  Range        detector_coord;
  double       sterrad;
  string       options;
  string       file_name;
 
  try {
 
    JParser<> zap("Auxiliary program to integrate atmospheric neutrino flux over apparent trajectory of a source.");
    
    zap['D'] = make_field(dec, "declination [deg]");
    zap['R'] = make_field(RA, "right ascension [deg]");
    zap['t'] = make_field(T_start_int, "start time + interval [MJD]");
    zap['E'] = make_field(E_min_max, "energy range [GeV]") = Range(1e3,1e7);
    zap['C'] = make_field(detector_coord, "detector coordinates [rad]") = Range(ARCA.getLatitude(),ARCA.getLongitude());
    zap['s'] = make_field(sterrad, "correction factor") = 1;
    zap['o'] = make_field(file_name) = "./formula.txt";
 
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }
 
  Flux_Atmospheric flux;
  
  double epoch_mean = 2000.0;
  
  dec *= M_PI/180.0;
  RA *= M_PI/180.0;
 
  vector<double> thetavec;
  vector<double> phivec;
 
  //setting steps for energy range and time integration
  const int NT = 10000;
  const int NE = 1000;
 
  //divide energy range in log-equidistant steps and initialize arrays
  const double xmax = log10(E_min_max.getUpperLimit());
  const double xmin = log10(E_min_max.getLowerLimit());
  double flx[NE];
  double E[NE];
  for (int i = 0; i < NE; ++i){
    flx[i] = 0;
    double x = (double(i))*(xmax-xmin)/(double(NE))+xmin;
    E[i] = pow(10.0,x);
  }
  
  //determine theta angles along path
  for (int i = 0; i < NT; ++i) {
    double RA_ap, dec_ap;
    double amprms[21];
 
    //make time range
    double time = T_start_int.getUpperLimit()*i;
    time /= (double(NT));
    time += T_start_int.getLowerLimit();
 
    //use time to determine LST, determine slalib parameters and use to transform coordinates from mean to apparent place
    double LST = slaGmst(time) + detector_coord.getUpperLimit();
    slaMappa(epoch_mean, time, amprms);
    slaMapqkz(LST-RA, dec, amprms, &RA_ap, &dec_ap);
 
    //transform to horizontal coordinate system and save values in vector
    double phi, theta;
    slaDe2h(RA_ap, dec_ap, detector_coord.getLowerLimit(), &phi, &theta);
    thetavec.push_back(theta);
    phivec.push_back(phi);
 
    //transform theta to cos zenith
    theta = cos(M_PI/2 - theta);
 
    //get fluxes for given theta and energy range
    for (int j = 0; j < NE; ++j){
      flx[j] += flux.dNdEdOmega(14, log10(E[j]), theta);
    }
  }
 
  //construct formula
  ostringstream oss_formula;
  
  //first terms separate to get bounds correct
  flx[0] *= (sterrad/NT) * pow(E[0],3.7);
  flx[1] *= (sterrad/NT) * pow(E[1],3.7);
 
  double rc = (flx[1]-flx[0])/(E[1]-E[0]);
  oss_formula << "((x>=" << E[0] << "&&x<=" << E[1] << ")*(" << rc << "*(x-" << E[0] << ")+" << flx[0] << ")+";
 
  //bulk of terms
  for (int i = 2; i < NE; ++i){
    flx[i] *= (sterrad/NT) * pow(E[i],3.7);
    rc = (flx[i]-flx[i-1])/(E[i]-E[i-1]);
    oss_formula << "(x>" << E[i-1] << "&&x<=" << E[i] << ")*(" << rc << "*(x-" << E[i-1] << ")+" << flx[i-1] << ")+";
  }
  //final term
  oss_formula << "(x>" << E[NE-1] << "&&x<=" << E_min_max.getUpperLimit() << ")*(" << rc << "*(x-" << E[NE-1] << ")+" << flx[NE-1] << "))*pow(x,-3.7)";
 
  //cout << oss_formula.str() << endl;
 
  //write formula to text file
  if (getFilenameExtension(file_name)=="txt"){
    ofstream out(file_name);
    out << oss_formula.str();
    out.close();
  }
  
  //make plots and save to rootfile
  if (getFilenameExtension(file_name)=="root"){
    TFile out(file_name.c_str(), "recreate");
 
    //formula as TF1
    auto f1 = new TF1("f1",oss_formula.str().c_str(),E_min_max.getLowerLimit(),E_min_max.getUpperLimit());
    f1->SetTitle("Linear interpolation");
    f1->Write("formula");  
 
    //plot fluxes in graph
    TGraph* graph = new TGraph();
    for (int i = 0; i < NE; ++i){
      double flx_norm = flx[i] * pow(E[i],-3.7);
      graph->SetPoint(i,E[i],flx_norm);
    }
 
    graph->GetXaxis()->SetTitle("E");
    graph->GetYaxis()->SetTitle("Flux");
    graph->Write("flux");
  
    //plot trajectory source
    TGraphErrors* graph_rec_coord = new TGraphErrors(thetavec.size(),&(phivec[0]),&(thetavec[0]));
    graph_rec_coord->SetLineWidth(0);
    graph_rec_coord->SetMarkerStyle(1);
 
    graph_rec_coord->GetXaxis()->SetTitle("Azimuth");
    graph_rec_coord->GetYaxis()->SetTitle("Altitude");
    graph_rec_coord->SetTitle("");
    graph_rec_coord->GetXaxis()->SetLimits(0,2*M_PI);
    graph_rec_coord-> Write("path");
 
    out.Write();
    out.Close();
  }
 
  return 0;
}