JAAMuonPostfit.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/offline/Trk.hh"
#include "km3net-dataformat/tools/reconstruction.hh"

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

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

#include "JAAnet/JHeadToolkit.hh"
#include "JAAnet/JAAnetToolkit.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 {

  struct JEvtWeightOption {

    static const char* const unweighted() { return "unweighted";   }
    static const char* const standard()   { return "standardFlux"; }
    static const char* const rescaled()   { return "rescaledFlux"; }
  };
}


/**
 * \file
 *
 * Example program to analyse track fit results from Evt formatted data.
 * \author bofearraigh, bjung
 */
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;
  JRange<unsigned int>   triggeredHitsRange;
  
  string                 weightOption;
  int                    debug;
  
  try {

    JProperties cuts;

    cuts.insert(gmake_property(energyRange));
    cuts.insert(gmake_property(coszenithRange));    
    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['@'] = make_field(cuts)                = JPARSER::initialised();
    zap['W'] = make_field(weightOption)        = 
      JEvtWeightOption::unweighted(),
      JEvtWeightOption::standard(),
      JEvtWeightOption::rescaled();
    zap['d'] = make_field(debug)               = 2;

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

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

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

  
  const JPMTRouter router(detector);  

  // 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(weightOption == JEvtWeightOption::unweighted() ? 
                       "Number of events" : "Rate [Hz]");

  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);
  
  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();

      if (has_reconstructed_jppmuon(*evt)) {
      
        DEBUG(get_best_reconstructed_jppmuon(*evt).comment << endl);
        
        const Trk&      best_trk = get_best_reconstructed_jppmuon(*evt);
        const JTrack3E& tb       = getTrack(best_trk);

        const double    logE           = log10(tb.getE());
        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) &&
            coszenithRange    (tb.getDZ())     &&
            energyRange       (tb.getE())) { // Apply cuts

          // Extract event-weight

          double weight = 0.0;
          
          if        (weightOption == JEvtWeightOption::unweighted()) {
            weight = 1.0;
          } else if (weightOption == JEvtWeightOption::standard()) {
            weight = scanner->getWeight(*evt);
          } else if (weightOption == JEvtWeightOption::rescaled()) {
            weight = evt->w.at(RESCALED_WEIGHT_INDEX);
          } 
      
          STATUS(RIGHT     (15)    << distance(scanners.begin(), scanner) <<
                 RIGHT     (15)    << scanner->getCounter()               <<
                 SCIENTIFIC(15, 3) << weight << '\r'); DEBUG(endl);

          // Fill histograms

          hs .Fill(NsnapshotHits,                          weight);
          hn .Fill(NtriggeredHits,                         weight);
          ho .Fill(evt->overlays,                          weight);
        
          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();
}