JSpark.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <set>
#include <vector>
#include <cmath>

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

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

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

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

#include "JDynamics/JDynamics.hh"

#include "JROOT/JManager.hh"

#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"

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


/**
 * \file
 *
 * Example program to check for sparks.
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTYPELIST<Evt, JMeta, TH1D, TH2D>::typelist  typelist;
  
  JMultipleFileScanner<Evt>               inputFile;
  JMultipleFileScanner<JACOUSTICS::JEvt>  calibrationFile;
  JFileRecorder<typelist>                 outputFile;
  JLimit_t&       numberOfEvents = inputFile.getLimit();
  string          detectorFile;
  double          Tmax_s;
  bool            use_weights;
  JDAQTriggerMask trigger_mask = TRIGGER_MASK_ON;
  set<double>     Dmax_m = { 0.3, 2.0 };
  int             debug;

  try {

    JProperties properties;

    properties.insert(gmake_property(trigger_mask));
    properties.insert(gmake_property(Dmax_m));

    JParser<> zap("Example program to check for sparks.");

    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['@'] = make_field(properties)          = JPARSER::initialised();
    zap['+'] = make_field(calibrationFile)     = JPARSER::initialised();
    zap['T'] = make_field(Tmax_s)              = 100.0;
    zap['o'] = make_field(outputFile);
    zap['W'] = make_field(use_weights);
    zap['d'] = make_field(debug)               = 2;

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

    
  JDetector detector;

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

  sort(detector.begin(), detector.end(), make_comparator(&JModule::getLocation));

  unique_ptr<JDynamics> dynamics;

  try {

    dynamics.reset(new JDynamics(detector, Tmax_s));

    dynamics->load(calibrationFile);
  }
  catch(const exception& error) {
    if (!calibrationFile.empty()) {
      FATAL(error.what());
    }
  }


  const floor_range   range = getRangeOfFloors(detector);
  const JStringRouter router(detector);

  TH1D h0("h0", NULL, 50, -1.0, 2.0);

  JManager<int, TH2D> H2(new TH2D("h2[%]", NULL, router.size(), -0.5, router.size() - 0.5, range.getUpperLimit(), 1 - 0.5, range.getUpperLimit() + 0.5));

  for (Int_t i = 1; i <= H2->GetXaxis()->GetNbins(); ++i) {
    H2->GetXaxis()->SetBinLabel(i, MAKE_CSTRING(router.at(i-1)));
  }
  for (Int_t i = 1; i <= H2->GetYaxis()->GetNbins(); ++i) {
    H2->GetYaxis()->SetBinLabel(i, MAKE_CSTRING(i));
  }


  outputFile.open();
  outputFile.put(JMeta(argc, argv));

  double W = 1.0;
    
  for (string file_name = ""; inputFile.hasNext(); ) {

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

    const Evt* evt = inputFile.next();

    if (JDAQTriggerMask(evt->trigger_mask).hasTriggerMask(trigger_mask)) {

      if (dynamics) {
        dynamics->update(evt->t.GetSec());
      }

      if (use_weights) {

        if (!evt->w.empty()) {
          W = evt->w[WEIGHTLIST_RUN_BY_RUN_WEIGHT];
        }
      
      } else {

        if (file_name != inputFile.getFilename()) {

          file_name = inputFile.getFilename();

          try {

            const JAANET::JHead head(getHeader(file_name));

            if (is_mupage(head))     { W = head.DAQ.livetime_s / head.livetime.numberOfSeconds; }
            if (is_pure_noise(head)) { W = head.DAQ.livetime_s / head.K40.livetime_s; }
          }
          catch(const exception& error) {}
        }
      }
    
      double    Dmin      = numeric_limits<double>::max();
      JLocation location(-1, -1);
      
      if (has_reconstructed_jppshower(*evt)) {

        const Trk      trk = get_best_reconstructed_jppshower(*evt);
        const JTrack3D ts  = getTrack(trk);
          
        if (trk.status != TRK_ST_UNDEFINED) {
      
          for (const auto& module : (dynamics ? dynamics->getDetector() : detector)) {

            if (module.getFloor() != 0) {

              const double D = getDistance(module.getPosition(), ts.getPosition());
          
              if (D < Dmin) {
                Dmin     = D;
                location = module.getLocation();
              }
            }
          }
        }

        if (Dmin != numeric_limits<double>::max()) {

          h0.Fill(log10(Dmin), W);

          const int i = distance(Dmax_m.begin(), Dmax_m.lower_bound(Dmin));
        
          H2   ->Fill((double) router.getIndex(location.getString()), (double) location.getFloor(), W);
          H2[i]->Fill((double) router.getIndex(location.getString()), (double) location.getFloor(), W);
        }
      }
    }    
    //outputFile.put(*evt);
  }
  STATUS(endl);


  outputFile.put(h0);
  outputFile.put(*H2);

  for (auto& h2 : H2) {
    outputFile.put(*h2.second);
  }
  
  outputFile.close();
}