JAAPostfit.cc File Reference

Example program to analyse track fit results from Evt formatted data. More...

#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/offline/Trk.hh"
#include "km3net-dataformat/tools/reconstruction.hh"
#include "JTools/JRange.hh"
#include "JPhysics/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/JTriggeredFileScanner.hh"
#include "JSupport/JSupport.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Example program to analyse track fit results from Evt formatted data.

Author

bofearraigh

Definition in file JAAPostfit.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 37 of file JAAPostfit.cc.

{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  JParallelFileScanner<Evt> inputFile;
  JLimit_t&  numberOfEvents = inputFile.getLimit();
  string     detectorFile;
  string     outputFile;
  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['d'] = make_field(debug)               = 2;

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

  const double rad_to_deg   =  180/M_PI;
  const double R_m          =  700.0;
  const int    MAX_OVERLAYS = 200;
  const double zenith_min   = -1.0;
  const double zenith_max   = 1.0;
    
  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 job("job", NULL, 100,   0.5, 100.5);
  TH1D hz  ("hz",   "; cos(#theta)"                  , 21 , zenith_min, zenith_max);                                                  
  TH1D ho  ("ho",   "; #overlays"                    , MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);                             
  TH2D hzo ("hzo",  "; cos(#theta) ; #overlays"      , 21, zenith_min, zenith_max, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);
  TH2D hxy ("hxy",  "; x position ; y position"      , 201, -R_m, R_m, 201, -R_m, R_m);               
  TH1D hq  ("hq",   "; quality parameter"            , 100, -200.0, 800);
  TH1D hb0 ("hb0",  "; log(beta0)"                   , 60, -2, 3.5);                                               
  TH1D he  ("he",   "; log(Ereco [GeV])"             , 75, -2, 4);                                                 
  TH2D heo ("heo",  "; log(Ereco [GeV]) ; #overlays" , 75, -2, 4, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5); 
  TH2D hzq ("hzq",  "; cos(#theta); quality"         , 21, zenith_min, zenith_max, 50, -10.0, 500.0);                    
  TH2D hze ("hze",  "; cos(#theta); log(Ereco [GeV])", 21, zenith_min, zenith_max, 75, -2, 4);                    
  TH2D hzb0("hzb0", "; cos(#theta); log(beta0)"      , 21, zenith_min, zenith_max, 60, -2, 3.5);     

  while (inputFile.hasNext()) {

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

    const Evt * evt = inputFile.next();

    if ( has_reconstructed_jppmuon(*evt)) {
      
      DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl);
      
      Trk best_trk = get_best_reconstructed_track<JPP_RECONSTRUCTION_TYPE>(*evt );

      JTrack3E tb = getTrack(get_best_reconstructed_jppmuon(*evt));

      const double Efit = tb.getE(); 

        ho .Fill(evt->overlays);
        hz .Fill(tb.getDZ());
        hzo.Fill(tb.getDZ(), evt->overlays) ;
        hxy.Fill(tb.getX(), tb.getY());
        hq .Fill(best_trk.lik);
        hzq.Fill(tb.getDZ(), best_trk.lik );
        
        if (best_trk.fitinf[JGANDALF_BETA0_RAD]){
        hb0.Fill(log10(rad_to_deg* best_trk.fitinf[JGANDALF_BETA0_RAD]) );
        hzb0.Fill(tb.getDZ(),  log10(rad_to_deg* best_trk.fitinf[JGANDALF_BETA0_RAD]) );
        }
        he .Fill(log10(Efit));
        hze.Fill(tb.getDZ(), log10(Efit) );
        heo.Fill(log10(Efit), evt->overlays );
      }
    }   
  STATUS(endl);

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