software/JReconstruction/JSpark.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <set>
#include <map>
 
#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/Evt.hh"
 
 
#include "JAAnet/JAAnetToolkit.hh"
#include "JAAnet/JHead.hh"
#include "JAAnet/JHeadToolkit.hh"
 
#include "JDAQ/JDAQSummarysliceIO.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JDynamics/JDynamics.hh"
 
#include "JROOT/JManager.hh"
 
#include "JReconstruction/JSparkSummary.hh"
 
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JSummaryFileRouter.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"
 
#include "JTrigger/JTriggerToolkit.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
 
/**
 * \file
 *
 * Program to monitor variables for spark removal in offline files. The variable definition is given in JSparkSummary class.
 * \author vcarretero
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  //typedef JTYPELIST<JMeta, TTree>::typelist  typelist;
  
  JSingleFileScanner<JTypeList<Evt, JDAQSummaryslice>>            inputFile;
  JMultipleFileScanner<JACOUSTICS::JEvt>                          calibrationFile;
  JFileRecorder<JTypeList<JMeta, JRECONSTRUCTION::JSparkSummary>> outputFile;
  JLimit_t&                                                       numberOfEvents = inputFile.getLimit();
  string                                                          detectorFile;
  double                                                          Tmax_s;
  int                                                             debug;
 
  try {
 
 
    JParser<> zap("Program to monitor variables for spark removal in offline files. The variable definition is given in JSparkSummary class.");
 
    zap['f'] = make_field(inputFile);
    zap['n'] = make_field(numberOfEvents)      = JLimit::max();
    zap['a'] = make_field(detectorFile);
    zap['+'] = make_field(calibrationFile)     = JPARSER::initialised();
    zap['T'] = make_field(Tmax_s)              = 100.0;
    zap['o'] = make_field(outputFile);
    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());
    }
  }
 
  JSummaryFileRouter  summary(inputFile);
 
  const JModuleRouter router(detector);
 
  outputFile.open();
  outputFile.put(JMeta(argc, argv));
 
  for (JTreeScanner<Evt> in(inputFile, numberOfEvents); in.hasNext(); ) {
 
    STATUS("event: " << setw(10) << in.getCounter() << '\r'); DEBUG(endl);
    
    JRECONSTRUCTION::JSparkSummary spark;
 
    const Evt* evt = in.next();
    
    summary.update(JDAQChronometer(evt->det_id,
                                   evt->run_id,
                                   evt->frame_index,
                                   JDAQUTCExtended(evt->t.GetSec(), evt->t.GetNanoSec() / 16))); 
 
    // Computation of max pmt getting a hit in a DOM for a given event, snapshot and triggered                                
    std::map<int,std::set<int>> snapshot, trigger;
 
    for (const Hit& hit : evt->hits) {
      bool valid = true;
      if(router.hasModule(hit.dom_id)){
        
        const JDAQSummaryFrame& frame = summary.getSummaryFrame(hit.dom_id);
        const JModule& module = router.getModule(hit.dom_id);
 
        valid = (getDAQStatus(frame, module, hit.channel_id)                &&         // check UDP
                 getPMTStatus(frame, module, hit.channel_id)                &&         // check if disabled, HRV or FIFO (almost) full
                 frame[hit.channel_id].is_valid()                           &&         // check if valid pmt
                 !module.getPMT(hit.channel_id).has(PMT_DISABLE));                     // check if disabled pmt
        
      }
      if(valid){
        snapshot[hit.dom_id].insert(hit.channel_id);
 
        if(hit.trig) trigger[hit.dom_id].insert(hit.channel_id);
      }
     
    }
    
    for (const auto& [dom_id, pmt_set] : snapshot) {
      if (pmt_set.size() > spark.nMaxPMT_snapshot) {
        
        spark.nMaxPMT_snapshot = pmt_set.size();
        spark.nMaxPMT_snapshot_ModuleID = dom_id;
 
      }
    }
 
    for (const auto& [dom_id, pmt_set] : trigger) {
      if (pmt_set.size() > spark.nMaxPMT_trigger) {
        
        spark.nMaxPMT_trigger = pmt_set.size();
        spark.nMaxPMT_trigger_ModuleID = dom_id;
      
      }
    }
 
    //Computation of the minimal distance from jshower vertex to a module
 
    if (dynamics) {
      dynamics->update(evt->t.GetSec());
    }
      
    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 < spark.Dmin) {
              
              spark.Dmin          = D;
              spark.Dmin_ModuleID = module.getID();
              
            }
          }
        }
      }
    } 
    outputFile.put(spark);
  } 
 
  JSingleFileScanner<JMeta> io(inputFile);
 
  io >> outputFile;
  STATUS(endl);
  outputFile.close();
 
}