JSupernovaProcessor.cc

Go to the documentation of this file.

#include <iostream>
#include <iomanip>
#include <queue>
#include <string>
 
#include "nlohmann/json.hpp"
 
#include "JDAQ/JDAQTimesliceIO.hh"
#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"
 
#include "JDAQ/JDAQTags.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
#include "JNet/JControlHostObjectIterator.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
#include "JROOT/JROOTClassSelector.hh"
 
#include "JTrigger/JDAQHitToTSelector.hh"
 
#include "JSupport/JSupport.hh"
 
#include "JSupernova.hh"
 
#include "JSupernovaProcessor.hh"
 
using json = nlohmann::json;
 
/**
 * \file
 *
 * Online supernova monitor
 * \author mlincett,selhedri,iagoos
 */
int main(int argc, char* argv[]) {
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  string    controlhost;
  string    ligier;
  string    outputFileName;
  bool      testMode = false;
  int       nTestEntries;
  JTag      outputTag;
  JTimeval  timeout_us;
  int       numberOfTimeouts;
  int       debug;
  int       queueLength;
  int       windowLength;
  int       statPrintInterval_s;
 
  string detectorFile;
  string dummyDetectorFile;
  JROOTClassSelector selector;
  string summaryFile;
 
  int       TMax_ns;
  int       preTriggerThreshold;
  JRange<int> M, M_BDT;
  double TVeto_ns;
  JDAQHitToTSelector     totSelector_ns;
 
  try {
 
    JParser<> zap("Supernova realtime monitor");
 
    zap['H'] = make_field(controlhost, "CH server (input)")         = "localhost";
    zap['L'] = make_field(ligier, "Ligier server (output)")         = "";
    zap['T'] = make_field(outputTag, "Output tag for the json messages") = JTag("SNT_part");
    zap['O'] = make_field(testMode, "Test mode if true (no output to Ligier)");
    zap['N'] = make_field(nTestEntries, "Number of test entries to write") = 1;
    zap['o'] = make_field(outputFileName, "Output json file name if run in test mode") = "";
    zap['t'] = make_field(timeout_us)       = 1e6;
    zap['n'] = make_field(numberOfTimeouts) = 1e3;
    zap['a'] = make_field(detectorFile);
    zap['C'] = make_field(selector)         = getROOTClassSelection<JDAQTimesliceTypes_t>();
    zap['Q'] = make_field(queueLength,  "number of timeslices of trigger queue")          = 100;
    zap['W'] = make_field(windowLength, "number of timeslices of trigger sliding window") = 5;
    zap['T'] = make_field(TMax_ns, "coincidence time window [ns]") = 10;
    zap['M'] = make_field(M, "multiplicity range for SN coincidences") = JRange<int>(6,10);
    zap['S'] = make_field(preTriggerThreshold, "muon veto multiplicity threshold") = 4;
    zap['V'] = make_field(TVeto_ns, "muon veto time interval") = 1000;
    zap['s'] = make_field(summaryFile, "summary output file") = "";
    zap['P'] = make_field(statPrintInterval_s, "statistics & file print interval [s]") = 30;
    zap['d'] = make_field(debug)            = 1;
    zap['b'] = make_field(dummyDetectorFile)   = "";
    zap['B'] = make_field(M_BDT, "multiplicity range for SN coincidences (with BDT)") = JRange<int>(5,10);
    zap['m'] = make_field(totSelector_ns) = JDAQHitToTSelector(4, 255);
 
    zap(argc, argv);
 
  } 
 
  catch(const exception &error) {
    FATAL(error.what() << endl);
  }
 
  if (queueLength < windowLength) {
    FATAL("Length of the trigger window must be smaller than the queue.");
  }
 
 
  setDAQLongprint(debug >= JEEP::debug_t);
 
  using namespace JSUPERNOVA;
 
  // -------------------------------
  // load detector and build routers
  // -------------------------------
 
  JDetector detector;
  JDetector dummyDetector;
 
  try {
    load(detectorFile, detector);
  }
  catch(const JException& error) {
    FATAL(error);
  }
 
  try {
    load(dummyDetectorFile, dummyDetector);
  }
  catch(const JException& error) {
    dummyDetector = JDetector();
  }
 
  const JDAQHitRouter hitRouter(detector);
 
  const JModuleRouter& moduleRouter = hitRouter; 
 
  const int DETID = detector.getID();
 
  const int detectorSize = detector.size();
 
  JSNFilterM F_M1(M, 1);
 
  // -------------------------------
  // initialize processing queue
  // -------------------------------
 
  typedef JTriggerSN trigger_type; 
 
  typedef priority_queue<trigger_type, vector<trigger_type>, greater<trigger_type> > queue_type;
 
  typedef deque<trigger_type> window_type; 
 
 
  
  typedef map<int, map<int, double> > rates_type; // (frameindex, DOMID) -> DOM rate [kHz]
  typedef map<int, int>               npmt_type ; // (frameindex)        -> # active PMTs
  typedef map<int, JVetoSet>          veto_type ; // (frameindex)        -> vetoes
 
  queue_type  trgQueue;
  window_type trgWindow;
  rates_type  rates;
  npmt_type   pmts;
  veto_type   veto;
 
  JTriggerSNStats stats(detectorSize);
 
  long int counter_live_ts = 0;
  long int counter_lost_ts = 0;
 
  double frameTime_s = getFrameTime() / 1.0e9;
 
  // Initialize output file name and number of entries for testing
  ofstream outputFile;
  int nWrittenEntries = 0;
  if (testMode){
    outputFile.open(outputFileName);
  }
 
  // -------------------------------
  // main processing
  // -------------------------------
  
  try {
    
    // setup input
 
    typedef JDAQTimesliceSN  data_type;
    typedef JDAQSummaryslice summary_type;
    typedef JDAQEvent        event_type;
 
    JControlHostObjectIterator<data_type>    in(controlhost, timeout_us, true);
 
    // timeout for the asynchronous reading of summary and event data
    // needs to be smaller than the timeslice duration
    const double asyncTimeout_us = 1000.0;
 
    JControlHostObjectIterator<summary_type> sm(controlhost, asyncTimeout_us, true);
    JControlHostObjectIterator<event_type>   ev(controlhost, asyncTimeout_us, true);
 
    JSinglePointer<JControlHost> out;
 
    if (!testMode && ligier != "") {
      out.reset(new JControlHost(ligier));
      out->MyId(argv[0]); // pid
    }
 
    // setup state
 
    int RUN = 0;
 
    for (int i = 0; i != numberOfTimeouts; ) {
 
      if (in.hasNext()) {
 
        data_type* timeslice = in.next();
        
        DEBUG(timeslice->getDAQHeader() << endl);
 
        int timesliceSize = timeslice->size();
        
        // --------------------------------
        // run number initialise and update
        // --------------------------------
 
        const int r = timeslice->getRunNumber();
 
        if (r != RUN) {
 
          if (RUN != 0) {
 
            NOTICE("RUN CHANGE" << endl);
 
            while (trgQueue.size() > 0) { trgQueue.pop(); }
 
            trgWindow.clear();  
 
            rates.clear();
 
            pmts.clear();
 
            veto.clear();
 
          }
 
          RUN = r;
          
        }
        
        // ------------------------------
        // process pending summary slices
        // ------------------------------       
 
        while ( sm.hasNext() ) {
 
          JDAQSummaryslice* summary = sm.next();
 
          DEBUG("SUM " << summary->getDAQHeader() << endl);
 
          int frame_index = summary->getFrameIndex();
 
          for (JDAQSummaryslice::const_iterator summary_frame = summary->begin();
               summary_frame != summary->end();
               ++summary_frame) {
 
            int DOMID = summary_frame->getModuleID();
 
            for (int ipmt = 0 ; ipmt < NUMBER_OF_PMTS ; ipmt++) {
              rates[frame_index][DOMID] += summary_frame->getRate(ipmt, 1.0/1000);
            }
 
            pmts[frame_index] += summary_frame->countActiveChannels();
          }
        
        }
 
        DEBUG("LOADING EVENTS" << endl);
 
        while ( ev.hasNext() ) {
 
          JDAQEvent* event = ev.next();
 
          DEBUG("EVT " << event->getDAQHeader() << endl);
          
          int frame_index = event->getFrameIndex();
 
          veto[frame_index].push_back(JVeto(*event, hitRouter));
 
        }
 
        // -----------------
        // process timeslice
        // -----------------
 
        JDataSN preTrigger(TMax_ns, preTriggerThreshold);
 
        preTrigger(timeslice, moduleRouter, totSelector_ns, dummyDetector);
 
        JTriggerSN trigger(TVeto_ns);
 
        trigger(preTrigger);
 
        trgQueue.push(trigger);
 
        //----------------
        // compute trigger
        //----------------
 
        vector<double> observables;
        vector<int> multiplicities;
 
        if (trgQueue.size() >= (unsigned) queueLength) {
          
          while (trgWindow.size() <= (unsigned) windowLength) {
 
            trigger_type pending = trgQueue.top();
              
            if ( trgWindow.size() == 0 || pending > trgWindow.back() ) {
              
              trgWindow.push_back( pending );
 
              counter_live_ts++;
 
            } else {
              // latecoming (out of order) timeslice 
              counter_lost_ts++;
            }
 
            trgQueue.pop();
          
          }
 
          // build triggered modules
          
          int trg_cc_counts   = 0;
          int trg_cc_modules  = 0;
          set<int> cc_modules;
 
          int trg_ev_counts   = 0;
          int trg_ev_modules  = 0;
          set<int> ev_modules;
 
          // loop over the trigger window and count the triggers
          for (int its = 0; its < windowLength; its++) {
            
            const int frame_index = trgWindow[its].frameIndex;
 
            JVetoSet vetoSet;
            if (veto.count(frame_index)) { 
              vetoSet = veto.at(frame_index);
            }
 
            JSNFilterMV F_MV(M, vetoSet);
            JSNFilterMV F_MV_BDT(M_BDT, vetoSet);
 
            set<int> cc_vec = trgWindow[its].getModules(F_M1);
            set<int> ev_vec = trgWindow[its].getModules(F_MV);
 
            cc_modules.insert(cc_vec.begin(), cc_vec.end());
            ev_modules.insert(ev_vec.begin(), ev_vec.end());
 
            trg_cc_counts += count_if(trgWindow[its].begin(), trgWindow[its].end(), F_M1);
            trg_ev_counts += count_if(trgWindow[its].begin(), trgWindow[its].end(), F_MV);
 
            // Write table of observables
            for (auto &trg: trgWindow[its]){
              auto sn_candidate = trg.getPeak();
              if (F_MV_BDT(sn_candidate)){
                multiplicities.push_back(sn_candidate.multiplicity);
                observables.push_back(sn_candidate.total_ToT);
                observables.push_back(sn_candidate.deltaT);
                observables.push_back(sn_candidate.mean_dir_norm);
                observables.push_back(sn_candidate.mean_dir_ctheta);
              }
            }
          }
 
          trg_cc_modules = cc_modules.size();
          trg_ev_modules = ev_modules.size();
          
          // trigger window slide of one element
 
          int currentFrame             = trgWindow[0].frameIndex;
          JDAQUTCExtended currentTime  = trgWindow[0].timeUTC;
 
          trgWindow.pop_front();
 
          // calculate trigger
          
          ++stats[trg_cc_counts];
 
          // calculate active modules
 
          int    activeModules =  -1;
          double detectorRate  = 0.0;
 
          if (!rates.empty()          &&
                rates.count(currentFrame)) {
              
            activeModules = 0;
 
            for (map<int, double>::const_iterator p = rates.at(currentFrame).begin(); 
                 p != rates.at(currentFrame).end(); p++ ) {
 
              detectorRate  += p->second;
 
              activeModules += (p->second > 0);
 
            }
          } else {
           
            activeModules = timesliceSize;
 
          }
 
          // build summary message
                    
          json jd;
 
          jd[detid_name]                     = DETID;
          jd[active_doms_name]               = activeModules;
          jd[detector_rate_name]             = int(detectorRate / 1000.0);
          jd[run_number_name]                = RUN;
          jd[frame_index_name]               = currentFrame;
          jd[daq_time_name]                  = to_string(currentTime);
          jd[trigger_name][cc_name][c_name]  = trg_cc_counts;
          jd[trigger_name][cc_name][m_name]  = trg_cc_modules;
          jd[trigger_name][ev_name][c_name]  = trg_ev_counts;
          jd[trigger_name][ev_name][m_name]  = trg_ev_modules;
          jd[active_pmts_name]               = pmts[currentFrame];
          jd[observables_name]               = observables;
          jd[multiplicities_name]            = multiplicities;
 
          string msg = jd.dump();
 
          DEBUG(msg << endl);
 
          if (testMode){
            if (nWrittenEntries < nTestEntries){
              // write summary information to output file (first nTestEntries)
              outputFile << msg << endl;
              nWrittenEntries++;
            } else {
              // When test entries written, close output file and exit loop
              outputFile.close();
              break;
            }
          } else {
            // send summary information to output ligier
 
            if (out != NULL && trg_ev_counts > 0) {
              out->PutFullString(outputTag, msg);
            }
          }
          
          // print stats
 
          if ( (counter_live_ts % ((int)(statPrintInterval_s / frameTime_s)) == 0 ) ) {
            
            double livetime = counter_live_ts * frameTime_s;
            
            stats.setLiveTime(livetime);
 
            NOTICE(endl);
            NOTICE(stats.toString());
            NOTICE("=> discarded out-of-order timeslices = " << counter_lost_ts << endl);          
 
            if (summaryFile != "") {
              ofstream of(summaryFile.c_str());
              of << stats.toSummaryFile();
              of.close();
            }
 
          }
 
        } else { 
          NOTICE("Filling trigger queue: " << trgQueue.size() << "/" << queueLength << '\r'); 
        }
 
      }  else {
 
        NOTICE("timeout " << setw(3) << i << endl);
 
        ++i;
      }
    }
  }
  
  catch(const JSocketException& error) {
    ERROR(error.what() << endl);
  }
 
}