JRate.cc

Go to the documentation of this file.

#include <string>
#include <iostream>

#include "TROOT.h"
#include "TFile.h"
#include "TH1D.h"
#include "TProfile.h"

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JDAQ/JDAQEventIO.hh"

#include "JAstronomy/JAstronomy.hh"
#include "JAstronomy/JAstronomyToolkit.hh"

#include "JSupport/JTriggeredFileScanner.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupport.hh"

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


namespace {

  /**
   * Neutrino flux for RXJ1713.
   *
   * \param  E     neutrino energy [GeV]
   * \return       flux            [GeV^-1 s^-1 m^-2]
   */
  inline double rxj1713(const double E)
  {
    static const double k = 16.80e-15;    // [GeV^-1 s^-1 cm^-2]
    static const double a = 1.72;
    static const double e = 2.10;         // [TeV]
    
    return 1e4 * k * pow(E*1e-3, -a) * exp(-sqrt(E*1e-3/e));
  }
}  

  
/**
 * \file
 *
 * Example program to plot event rates using JASTRONOMY::JStarTrek.
 * \author mdejong
 */
int main(int argc, char **argv)
{
 using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JTriggeredFileScanner<> inputFile;
  JLimit_t&      numberOfEvents = inputFile.getLimit();
  string         outputFile;
  double         ct;
  double         numberOfBlocks;
  bool           join;
  int            debug;

  try {

    JParser<> zap("Example program to plot event rates using JASTRONOMY::JStarTrek.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "volume.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['Z'] = make_field(ct)                  = -0.25;
    zap['N'] = make_field(numberOfBlocks)      =  1.0;
    zap['J'] = make_field(join);                          // join neutrino and anti-neutrino
    zap['d'] = make_field(debug)               = 1;

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

  const double radius = 1.0;                                   // size of source [deg]
  const double omega  = 2*PI * (1.0 - cos(radius*PI/180));     // solid angle of source

  const JStarTrek startrek(Sicily, RXJ1713);


  Head head;

  try {
    head = getHeader(inputFile);
  }
  catch(const JException& error) {
    FATAL(error);
  }


  double Wall = numberOfBlocks;

  {
    JHead buffer(head);

    if (!buffer.is_valid(&JHead::genvol) || buffer.genvol.numberOfEvents == 0) {
      FATAL("No generated events." << endl);
    }

    Wall /= buffer.genvol.numberOfEvents;

    STATUS("Generation: " << buffer.genvol.numberOfEvents << endl);

    if (buffer.is_valid(&JHead::cut_nu)) {
      STATUS("Solid angle: " << 2 * (buffer.cut_nu.cosTmax - buffer.cut_nu.cosTmin) << "pi" << endl);
    }
  }

  if (join) {
    Wall *= 0.5;   // join neutrino and anti-neutrino
  }


  // output

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

  TH1D hc("ct [live time]", NULL, 1000, -1.0, +1.0);

  for (int i = 1; i != hc.GetNbinsX(); ++i) {

    const double x = hc.GetBinCenter(i);
    const double y = startrek(x);

    hc.SetBinContent(i, y);
  }

  TH1D h0("h0 [background]", NULL, 20, +1.0, +6.0);
  TH1D h1("h1 [signal]",     NULL, 20, +1.0, +6.0);

  TH1D*  H[] = { &h0, &h1 };
  double W[] = { 0.0, 0.0 };

  for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
    H[i]->Sumw2();
  }


  while (inputFile.hasNext()) {

    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl);

    JTriggeredFileScanner<>::multi_pointer_type ps = inputFile.next();

    //const JDAQEvent* tev   = ps;
    const Evt*       event = ps;

    if (has_neutrino(*event)) {
      
      const Trk& neutrino = get_neutrino(*event);
      
      const double E  = neutrino.E;         // GeV
      const double dz = neutrino.dir.z;     // cos(zenith angle)
      const double P  = 1.0;                // P(Earth)?
      const double W2 = event->w[1];        // generation weight GeV x m^2 x sr x (sec/year)
      const double W3 = event->w[2];        // number of atmospheric neutrino events / year
      
      // number of events / year
      
      const double x   = log10(E);
      const double y[] = { P * W3 * omega      * startrek(dz),
                           P * W2 * rxj1713(E) * startrek(dz) };
      
      if (dz >= ct) {   
        for (int i = 0; i != sizeof(H)/sizeof(H[0]); ++i) {
          H[i]->Fill(x, y[i] * Wall);   
          W[i] +=       y[i] * Wall;
        }
      }
    }
  }
  STATUS(endl);

  for (int i = 0; i != sizeof(W)/sizeof(W[0]); ++i) {
    NOTICE("W[" << i << "] = " << FIXED(9,2) << W[i] << endl);
  }

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