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;
  JRange<double>                 quality_range;
  JRange<double>                 energy_range;
  JRange<double>                 coszenith_range;
  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['Q'] = make_field(quality_range)       =  JRange<double>();
    zap['E'] = make_field(energy_range)        =  JRange<double>();
    zap['c'] = make_field(coszenith_range)     =  JRange<double>(-1.0,1.0);
    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 =  300;
  const double MAX_TRIGGERED_HITS   = 2500;
  const double E_RANGE_MIN   = -1.0;
  const double E_RANGE_MAX   = 7;
  const double Z_RANGE_MIN   = coszenith_range.getLowerLimit();
  const double Z_RANGE_MAX   = coszenith_range.getUpperLimit();
 
  TFile out(outputFile.c_str(), "recreate");
  TH1D job("job", NULL, 100,   0.5, 100.5);
  TH1D hz  ("hz",   "; cos(#theta_{zenith})"                     , 21 , Z_RANGE_MIN, Z_RANGE_MAX);                                                  
  TH1D ho  ("ho",   "; # of overlays"                            , MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5);                             
  TH2D hzo ("hzo",  "; cos(#theta_{zenith}) ; # of overlays"     , 21, Z_RANGE_MIN, Z_RANGE_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"                      , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
  TH1D hN  ("hN",  ";  JEnergy # of hits"                        , 600, -0.5, MAX_TRIGGERED_HITS-0.5);
  TH1D hq  ("hq",   "; quality parameter"                        , 200, -200.0, 800);
  TH1D hb0 ("hb0",  "; log(beta0)"                               , 60, -2.0, 3.5);                                               
  TH1D he  ("he",   "; log(Ereco [GeV])"                         , 75, E_RANGE_MIN, E_RANGE_MAX);                                                 
  TH2D heo ("heo",  "; log(Ereco [GeV]) ; # of overlays"         , 75, E_RANGE_MIN, E_RANGE_MAX, MAX_OVERLAYS , -0.5, MAX_OVERLAYS-0.5); 
  TH2D hen ("hen",  "; log(Ereco [GeV]) ; # of triggered hits"   , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5); 
  TH2D heN ("heN",  "; log(Ereco [GeV]) ; JEnergy # of hits"     , 75, E_RANGE_MIN, E_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5); 
  TH2D hzq ("hzq",  "; cos(#theta_{zenith}); quality"            , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -200.0, 800.0);                    
  TH2D hzn ("hzn",  "; cos(#theta_{zenith}); # of triggered hits", 21, Z_RANGE_MIN, Z_RANGE_MAX, 600 , -0.5, MAX_TRIGGERED_HITS-0.5);
  TH2D hzN ("hzN",  "; cos(#theta_{zenith}); JEnergy # of hits"  , 21, Z_RANGE_MIN, Z_RANGE_MAX, 600, -0.5, MAX_TRIGGERED_HITS-0.5);
  TH2D hze ("hze",  "; cos(#theta_{zenith}); log(Ereco [GeV])"   , 21, Z_RANGE_MIN, Z_RANGE_MAX, 75, E_RANGE_MIN, E_RANGE_MAX);                    
  TH2D hzb0("hzb0", "; cos(#theta_{zenith}); log(beta0)"         , 21, Z_RANGE_MIN, Z_RANGE_MAX, 60, E_RANGE_MIN, 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_application(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);
        
        if (coszenith_range(track.getDZ())  &&
            quality_range  (fit->getQ())    &&
            energy_range   (fit->getE())) {
 
          const double overlays = tev->getOverlays();
          const double Efit  = fit->getE();
 
          ho .Fill(overlays);
          hn .Fill((double) tev->size<JDAQTriggeredHit>());
          hz .Fill(track.getDZ());
          hzn.Fill(track.getDZ(), (double) tev->size<JDAQTriggeredHit>());
          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) );
            heN .Fill(log10(Efit), fit->getW(JENERGY_NUMBER_OF_HITS) );
            hzN .Fill(track.getDZ(), 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)) );
          }
          he  .Fill(log10(Efit));
          hen .Fill(log10(Efit), (double) tev->size<JDAQTriggeredHit>());
          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();
}