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 "km3net-dataformat/online/JDAQClock.hh"
 
#include "JAAnet/JAAnetToolkit.hh"
#include "JTools/JRange.hh"
#include "JMath/JMathToolkit.hh"
#include "JTrigger/JHitL0.hh"
#include "JPhysics/JConstants.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"
 
namespace {
 
  using JTRIGGER::JHitL0;
  using JDETECTOR::JModuleRouter;
 
  /**
   * Get transformation.
   *
   * \param  hit            hit
   * \param  router         module router
   * \return                hit
   */
  inline JHitL0 getHit(const Hit& hit, const JModuleRouter& router)
  {
    using namespace JPP;
 
    const JPMT& pmt = router.getModule(hit.dom_id).getPMT(hit.channel_id);
 
    return JHitL0(JDAQPMTIdentifier(hit.dom_id, hit.channel_id), pmt.getAxis(), JHit(getTime(hit.tdc, pmt.getCalibration()), getToT(hit.tot, pmt.getCalibration())));
  }
}
 
 
/**
 * \file
 *
 * Example program to analyse track fit results from Evt 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 Evt 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();
}