JAAnet.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <cmath>

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

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"

#include "JTools/JRange.hh"
#include "JTools/JConstants.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JMath/JMathToolkit.hh"
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSupport.hh"

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


/**
 * \file
 *
 * Example program to analyse track fit results from AAnet formatted data.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTOOLS::JRange<double> JRange_t;

  JMultipleFileScanner<Evt> inputFile;
  JLimit_t&  numberOfEvents = inputFile.getLimit();
  string     detectorFile;
  string     outputFile;
  JRange_t   E;
  int        debug;

  try {

    JParser<> zap("Example program to analyse track fit results from AAnet formatted data.");

    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile)        = "";
    zap['o'] = make_field(outputFile)          = "histogram.root";
    zap['E'] = make_field(E)                   = JRange_t();
    zap['d'] = make_field(debug)               = 2;

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

    
  JDetector detector;

  if (detectorFile != "") {
    try {
      load(detectorFile, detector);
    }
    catch(const JException& error) {
      FATAL(error);
    }
  }

  const JModuleRouter router(detector);

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

  TH1D hx("hx",  NULL,  100,   -3.0,  +2.3);    // [log(deg)]
  TH1D hd("hd",  NULL,  100,    0.0,  10.0);    // [m]
  TH1D ht("ht",  NULL,  100, -100.0, 100.0);    // [ns]
  TH1D he("he",  NULL,  100,   -5.0,  +5.0);    // [log10(E)]
  TH1D h1("h1",  NULL,  100,  -50.0, +50.0);    // [ns]

  
  while (inputFile.hasNext()) {

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

    const Evt* evt = inputFile.next();

    //if (has_muon(*evt) && has_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt, rec_stages_range(JMUONBEGIN, JMUONEND))) {
    if (has_muon(*evt) && has_reconstructed_jppmuon(*evt)) {

      DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl);

      JTrack3E ta = getTrack(get_muon(*evt));
      //JTrack3E tb = getTrack(get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt, rec_stages_range(JMUONBEGIN, JMUONEND)));
      JTrack3E tb = getTrack(get_best_reconstructed_jppmuon(*evt));

      ta.add(getTimeSinceRTS(evt->mc_t));

      ta.move(ta.getIntersection(tb), getSpeedOfLight(), gWater);
      tb.move(tb.getIntersection(ta), getSpeedOfLight(), gWater);

      if (E(ta.getE())) {

        hx.Fill(log10(getAngle(ta.getDirection(), tb.getDirection())));
        hd.Fill((ta.getPosition() - tb.getPosition()).getLength());
        ht.Fill( ta.getT() - tb.getT());
        he.Fill(log10(tb.getE()/ta.getE()));
      }

      for (vector<Hit>::const_iterator i = evt->hits.begin(); i != evt->hits.end(); ++i) {

        if (router.hasModule(i->dom_id)) {

          const JHitL0 hit = getHit(*i, router);

          h1.Fill(hit.getT() - tb.getT(hit.getPosition()));
        }
      }
    }
  }
  STATUS(endl);

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