JAAPostfit.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/definitions/weightlist.hh"

#include "km3net-dataformat/tools/reconstruction.hh"

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

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

#include "JDetector/JDetector.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JDetectorToolkit.hh"

#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JEvtWeightToolkit.hh"

#include "JPhysics/JConstants.hh"

#include "JGeometry3D/JGeometry3DToolkit.hh"

#include "JSupport/JLimit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JEvtWeightFileScannerSet.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"

#include "JTools/JRange.hh"


namespace {

  inline const char* const track_t()  { return "track"; }  
  inline const char* const shower_t() { return "shower"; }
}


/**
 * \file
 *
 * Example program to analyse fit results from Evt formatted data.
 *
 * \author bofearraigh, bjjung
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  JMultipleFileScanner_t     inputFiles;
  
  JLimit_t                   numberOfEvents;
  
  string                     detectorFile;
  string                     outputFile;

  JRange<double>             energyRange;
  JRange<double>             coszenithRange;
  JCylinder3D                containmentVolume = getMaximumContainmentVolume();
  JRange<unsigned int>       triggeredHitsRange;

  bool                       useWeights;
  JFluxMap                   fluxMaps;
  
  string                     oscProbTable;
  JOscParameters<>           oscParameters;

  string                     recotype;
    
  int                        debug;
  
  try {

    JProperties cuts;

    cuts.insert(gmake_property(energyRange));
    cuts.insert(gmake_property(coszenithRange));    
    cuts.insert(gmake_property(containmentVolume));
    cuts.insert(gmake_property(triggeredHitsRange));

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

    zap['f'] = make_field(inputFiles);
    zap['n'] = make_field(numberOfEvents)
      = JLimit::max();
    zap['a'] = make_field(detectorFile)
      = "";
    zap['o'] = make_field(outputFile)
      = "histogram.root";
    zap['R'] = make_field(recotype) =
      track_t(),    
      shower_t();
    zap['W'] = make_field(useWeights);
    zap['@'] = make_field(fluxMaps)
      = JPARSER::initialised();
    zap['P'] = make_field(oscProbTable)
      = JPARSER::initialised();
    zap['#'] = make_field(oscParameters)
      = JPARSER::initialised();
    zap['%'] = make_field(cuts)
      = JPARSER::initialised();
    zap['d'] = make_field(debug)
      = 2;

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

  if (useWeights && numberOfEvents != JLimit::max()) {
    FATAL("Cannot apply weighting to limited number of events.");
  }


  // Load detector file
  
  JDetector detector;

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

  
  // Set reco getter function pointer

  bool (*has_reconstruction)(const Evt&);
  const Trk& (*get_best_fit)(const Evt&);

  if        (recotype == track_t()) {
    has_reconstruction = &has_reconstructed_jppmuon;
    get_best_fit = &get_best_reconstructed_jppmuon;
  } else if (recotype == shower_t()) {
    has_reconstruction = &has_reconstructed_jppshower;
    get_best_fit = &get_best_reconstructed_jppshower;
  } else {
    FATAL("Invalid recotype: " << recotype);
  }
  

  // Histogram bins

  const Int_t    N_Nhits       =               400;
  const Double_t Nhits_min     =         0.0 - 0.5;
  const Double_t Nhits_max     =      1000.0 - 0.5;
  
  const Int_t    N_Noverlays   =                40;
  const Double_t Noverlays_min =         0.0 - 0.5;
  const Double_t Noverlays_max =      1000.0 - 0.5;

  const Int_t    N_radius      =               201;
  const Double_t radius_min    =      -700.0 - 0.5;
  const Double_t radius_max    =      +700.0 + 0.5;

  const Int_t    N_height      =                18;
  const Double_t height_min    =              20.0;
  const Double_t height_max    =             200.0;

  const Int_t    N_ToTfrac     =                20;
  const Double_t ToTfrac_min   =               0.0;
  const Double_t ToTfrac_max   =               1.0;

  const Int_t    N_dt          =              1500;
  const Double_t dt_min        =      -500.0 - 0.5;
  const Double_t dt_max        =      1000.0 - 0.5;

  const Int_t    N_zenith      =                21;
  const Double_t zenith_min    =              -1.0;
  const Double_t zenith_max    =              +1.0;
  
  const Int_t    N_energy      =                60;
  const Double_t energy_min    =              -2.0;
  const Double_t energy_max    =               4.0;

  const Int_t    N_quality     =               100;
  const Double_t quality_min   =            -200.0;
  const Double_t quality_max   =            +800.0;

  const Int_t    N_beta0       =                60;
  const Double_t beta0_min     =              -2.0;
  const Double_t beta0_max     =              +3.5;
  

  // Set of histograms for data/MC-comparisons

  string ordinateLabel(useWeights ? "Rate [Hz]" : "Number of events");

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

  TH1D hs  ("hs",   MAKE_CSTRING("; #snapshot hits; "                                  << ordinateLabel),
            N_Nhits,     Nhits_min,     Nhits_max);
  TH1D hn  ("hn",   MAKE_CSTRING("; #triggered hits; "                                 << ordinateLabel),
            N_Nhits,     Nhits_min,     Nhits_max);
  TH1D ho  ("ho",   MAKE_CSTRING("; #overlays; "                                       << ordinateLabel),
            N_Noverlays, Noverlays_min, Noverlays_max);
  
  TH1D hz  ("hz",   MAKE_CSTRING("; cos(#theta); "                                     << ordinateLabel),
            N_zenith,    zenith_min,    zenith_max);
  TH1D he  ("he",   MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); "                      << ordinateLabel),
            N_energy,    energy_min,    energy_max);
  TH1D ht  ("ht",   MAKE_CSTRING("; #Delta t [ns]; "                                   << ordinateLabel),
            N_dt,        dt_min,        dt_max);

  TH1D hZ  ("hZ",   MAKE_CSTRING("; longitudinal ToT-barycenter [m]; "                 << ordinateLabel),
            N_height,    height_min,    height_max);
  TH1D hT0 ("hT0",  MAKE_CSTRING("; ToT-fraction below earliest triggered hit [ns]; "  << ordinateLabel),
            N_ToTfrac,   ToTfrac_min,   ToTfrac_max);
  TH1D hT1 ("hT1",  MAKE_CSTRING("; ToT-fraction above earliest triggered hit [ns]; "  << ordinateLabel),
            N_ToTfrac,   ToTfrac_min,   ToTfrac_max);  

  TH1D hq  ("hq",   MAKE_CSTRING("; quality; "                                         << ordinateLabel),
            N_quality,   quality_min,   quality_max);
  TH1D hb0 ("hb0",  MAKE_CSTRING("; #beta_{0}; "                                       << ordinateLabel),
            N_beta0,     beta0_min,     beta0_max);
  
  TH2D hzo ("hzo",  MAKE_CSTRING("; cos(#theta); #overlays; "                          << ordinateLabel),
            N_zenith,    zenith_min,    zenith_max,
            N_Noverlays, Noverlays_min, Noverlays_max);
  TH2D hze ("hze",  MAKE_CSTRING("; cos(#theta); log_{10}(E_{reco} / 1 GeV); "         << ordinateLabel),
            N_zenith,    zenith_min,    zenith_max,
            N_energy,    energy_min,    energy_max);
  TH2D heo ("heo",  MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #overlays; "           << ordinateLabel),
            N_energy,    energy_min,    energy_max,
            N_Noverlays, Noverlays_min, Noverlays_max);
  TH2D hzt ("hzt",  MAKE_CSTRING("; cos(#theta); #Delta t [ns]; "                      << ordinateLabel),
            N_zenith,    zenith_min,    zenith_max,
            N_dt,        dt_min,        dt_max);
  TH2D het ("het",  MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); #Delta t [ns]; "       << ordinateLabel),
            N_energy,    energy_min,    energy_max,
            N_dt,        dt_min,        dt_max);
  TH2D hot ("hot",  MAKE_CSTRING("; #overlays; #Delta t [ns]; "                        << ordinateLabel),
            N_Noverlays, Noverlays_min, Noverlays_max,
            N_dt,        dt_min,        dt_max);
  
  TH2D hzq ("hzq",  MAKE_CSTRING("; cos(#theta); quality; "                            << ordinateLabel),
            N_zenith,    zenith_min,    zenith_max,
            N_quality,   quality_min,   quality_max);
  TH2D heq ("heq",  MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); quality; "             << ordinateLabel),
            N_energy,    energy_min,    energy_max,
            N_quality,   quality_min,   quality_max);
  TH2D hoq ("hoq",  MAKE_CSTRING("; #overlays; quality; "                              << ordinateLabel),
            N_Noverlays, Noverlays_min, Noverlays_max,
            N_quality,   quality_min,   quality_max);
  
  TH2D hob0("hob0", MAKE_CSTRING("; #overlays; log_{10}(#beta_{0}); "                  << ordinateLabel),
            N_Noverlays, Noverlays_min, Noverlays_max,
            N_beta0,     beta0_min,     beta0_max);
  TH2D heb0("heb0", MAKE_CSTRING("; log_{10}(E_{reco} / 1 GeV); log_{10}(#beta_{0}); " << ordinateLabel),
            N_energy,    energy_min,    energy_max,
            N_beta0,     beta0_min,     beta0_max);     
  TH2D hzb0("hzb0", MAKE_CSTRING("; cos(#theta); log_{10}(#beta_{0}); "                << ordinateLabel),
            N_zenith,    zenith_min,    zenith_max,
            N_beta0,     beta0_min,     beta0_max);
  
  TH2D hxy ("hxy",  MAKE_CSTRING("; x [m]; y [m]; "                                    << ordinateLabel),
            N_radius,    radius_min,    radius_max,
            N_radius,    radius_min,    radius_max);
  

  JEvtWeightFileScannerSet<> scanners(inputFiles);

  if (!oscProbTable.empty()) {

    JOscProbInterpolator<> interpolator(oscProbTable.c_str(), oscParameters);

    fluxMaps.configure(make_oscProbFunction(interpolator));
  }

  if (scanners.setFlux(fluxMaps) == 0) {
    WARNING("No flux function set." << endl);
  }
  
  STATUS(RIGHT(15) << "header ID" << RIGHT(15) << "event" << RIGHT(15) << "weight" << endl);
  
  for (JEvtWeightFileScannerSet<>::iterator scanner = scanners.begin(); scanner != scanners.end(); ++scanner) {
    
    scanner->setLimit(numberOfEvents);
    
    while (scanner->hasNext()) {
      
      const Evt* evt = scanner->next();
      
      const double weight = (useWeights ? scanner->getWeight(*evt) : 1.0);

      const size_t NsnapshotHits  = evt->hits.size();
      const size_t NtriggeredHits = count_if(evt->hits.begin(), evt->hits.end(),
                                             make_predicate(&Hit::trig, (long long unsigned int) 0, JComparison::gt()));
      
      if (!triggeredHitsRange(NtriggeredHits)) { continue; }
      
      hs .Fill(NsnapshotHits,  weight);
      hn .Fill(NtriggeredHits, weight);
      ho .Fill(evt->overlays,  weight);
      
      if (!has_reconstruction(*evt)) { continue; }
      
      const Trk&      best_trk = (*get_best_fit)(*evt);
      const JTrack3E& tb       = getTrack(best_trk);

      const double    logE     = log10(tb.getE());
        
      if (!energyRange   (tb.getE())                    ||
          !coszenithRange(tb.getDZ())                   ||
          !containmentVolume.is_inside(tb.getPosition())) { // Apply cuts
        continue;
      }

      // Extract event-weight
          
      STATUS(RIGHT     (15)    << distance(scanners.begin(), scanner) <<
             RIGHT     (15)    << scanner->getCounter()               <<
             SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);

      // Fill histograms
        
      hz .Fill(tb.getDZ(),                             weight);
      he .Fill(logE,                                   weight);
      hq .Fill(best_trk.lik,                           weight);

      hzo.Fill(tb.getDZ(),            evt->overlays,   weight);
      hze.Fill(tb.getDZ(),            logE,            weight);
      heo.Fill(logE,                  evt->overlays,   weight);

      hoq.Fill(evt->overlays,         best_trk.lik,    weight);
      hzq.Fill(tb.getDZ(),            best_trk.lik,    weight);
      heq.Fill(logE,                  best_trk.lik,    weight);

      hxy.Fill(tb.getX(),             tb.getY(),       weight);

      if (best_trk.fitinf[JGANDALF_BETA0_RAD]) {

        const double logb0 = log10(180.0 / M_PI * best_trk.fitinf[JGANDALF_BETA0_RAD]);
          
        hb0 .Fill(logb0,                               weight);
        hob0.Fill(evt->overlays,      logb0,           weight);   
        heb0.Fill(logE,               logb0,           weight);
        hzb0.Fill(tb.getDZ(),         logb0,           weight);
      }

      double ToTbelow     = 0.0;
      double ToTabove     = 0.0;
      double ToTtotal     = 0.0;
      double ToTweightedZ = 0.0;
        
      vector<Hit>::const_iterator firstHit = min_element(evt->hits.cbegin(), evt->hits.cend(),
                                                         make_comparator(&Hit::t));
        
      for (vector<Hit>::const_iterator hit = evt->hits.cbegin(); hit != evt->hits.cend(); ++hit) {
          
        if (router.hasPMT(hit->pmt_id)) {
            
          const JPMT&     pmt  = router.getPMT(hit->pmt_id);
            
          const double    t0   = tb.getT(pmt);
          const double    t1   = getTime(*hit);
            
          ht .Fill(t1 - t0,                            weight);
          hot.Fill(evt->overlays,     t1 - t0,         weight);
          hzt.Fill(tb.getDZ(),        t1 - t0,         weight);
          het.Fill(log10(best_trk.E), t1 - t0,         weight);
              
          if (hit->trig > 0) {

            ToTtotal     += hit->tot;
            ToTweightedZ += hit->tot * hit->pos.z;
                
            if        (hit->pos.z < firstHit->pos.z)  {
              ToTbelow   += hit->tot;
            } else if (hit->pos.z >= firstHit->pos.z) {
              ToTabove   += hit->tot;
            }
          }
        }
      }
          
      hT0.Fill(ToTbelow / ToTtotal,                     weight);
      hT1.Fill(ToTabove / ToTtotal,                     weight);
      hZ .Fill(ToTweightedZ / ToTtotal,                 weight);
    }
  }
  
  STATUS(endl);

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

  return 0;
}