JDataPostfit.cc

Go to the documentation of this file.

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

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

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/online/JDAQModuleIdentifier.hh"
#include "km3net-dataformat/definitions/reconstruction.hh"
#include "km3net-dataformat/definitions/fitparameters.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"

#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JVolume.hh"

#include "JGeometry3D/JTrack3D.hh"
#include "JGeometry3D/JGeometry3DToolkit.hh"
#include "JGeometry3D/JCylinder3D.hh"

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

#include "JDAQ/JDAQEventIO.hh"

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

/**
 * \file
 *
 * Program to histogram fit results from reconstructed data
 * \author bofearraigh, rgracia
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;

  using namespace KM3NETDAQ;

  typedef JTriggeredFileScanner<JEvt>                       JTriggeredFileScanner_t;
  typedef JTriggeredFileScanner_t::multi_pointer_type       multi_pointer_type;
  
  JTriggeredFileScanner<JEvt>    inputFile;                               
  JLimit_t&                      numberOfEvents = inputFile.getLimit();
  string                         outputFile;
  string                         option;
  size_t                         numberOfPrefits;                         // The number of reconstructed tracks to plot per event (1 by default)
  JQualitySorter                 quality_sorter;                          // The criterion by which to pick the "best" reconstructed track for each event.
  int                            application;
  int                            debug;

  try {

    JParser<> zap("Program to histogram fit results from reconstructed data.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "postfit.root";
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['N'] = make_field(numberOfPrefits)     = 1;
    zap['A'] = make_field(application)         = JMUONENERGY, JMUONPREFIT, JMUONSIMPLEX, JMUONGANDALF, JMUONSTART; //makes histograms after this stage of the reconstruction 
    zap['L'] = make_field(quality_sorter)      = JPARSER::initialised();
    zap['d'] = make_field(debug)               = 2;

    zap(argc, argv);
  }

  catch(const JException& error) {
    FATAL(error);
  }

  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;
  const double MAX_TRIGGERED_HITS   = 400;

  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 hn  ("hn",  ";  # of triggered hits"          , MAX_TRIGGERED_HITS, -0.5, MAX_TRIGGERED_HITS-0.5);                                               
  TH1D hN  ("hN",  ";  JEnergy # of hits"            , MAX_TRIGGERED_HITS, -0.5, MAX_TRIGGERED_HITS-0.5);
  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()) {               // Loop over each reconstructed event
    
    STATUS("event: " << setw(10) << inputFile.getCounter() << '\r'); DEBUG(endl); // if debug > set value, print out this line

    multi_pointer_type ps = inputFile.next();

    JDAQEvent* tev   = ps;
    JEvt*      evt   = ps;

    job.Fill(1.0);

    if (!evt->empty()) {

      JEvt::iterator __end = partition(evt->begin(), evt->end(), JHistory::is_event(application));

      if (evt->begin() == __end) {
        continue;
      }
      
      if (numberOfPrefits > 0) {        

        JEvt::iterator __q = __end;

        advance(__end = evt->begin(), min(numberOfPrefits, (size_t) distance(evt->begin(), __q)));
        
        partial_sort(evt->begin(), __end, __q, quality_sorter);

      } else {

        sort(evt->begin(), __end, quality_sorter);
      }
    
      if (numberOfPrefits > 0) {
        advance(__end = evt->begin(), min(numberOfPrefits, evt->size()));
      }

      for (JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
        
        JTrack3D track = getTrack(*fit);
        
        const double overlays = tev->getOverlays();

        ho .Fill(overlays);
        hn.Fill((double) tev->size<JDAQTriggeredHit>());
        hz .Fill(track.getDZ());
        hzo.Fill(track.getDZ(), overlays) ;
        hxy.Fill(track.getX(), track.getY());
        hq .Fill(fit->getQ());
        hzq.Fill(track.getDZ(), fit->getQ() );
        
        if (fit->hasW(JENERGY_NUMBER_OF_HITS)) { 
          hN .Fill( fit->getW(JENERGY_NUMBER_OF_HITS) );
        }

        if (fit->hasW(JGANDALF_BETA0_RAD)) { 
          hb0 .Fill(log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD) ));
          hzb0.Fill(track.getDZ(), log10( rad_to_deg* fit->getW(JGANDALF_BETA0_RAD)) );
        }
        const double Efit  = fit->getE();
         he .Fill(log10(Efit));
         hze.Fill(track.getDZ(), log10(Efit) );
         heo.Fill(log10(Efit), overlays );
      }
    }
  }
  STATUS("Number of events input             " << setw(8) << right << job.GetBinContent(1) << endl);
  out.Write();
  out.Close();
}