JAnalysisSummary.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 "JGeometry3D/JCylinder3D.hh"
 
#include "JROOT/JManager.hh"
 
#include "JReconstruction/JAnalysisSummary.hh"
 
#include "JSirene/JVisibleEnergyToolkit.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 produce analysis variables in offline files. The variable definition is given in JAnalysisSummary class. 
 * \author vcarretero
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
  
  JSingleFileScanner<JTypeList<Evt, JDAQSummaryslice>>               inputFile;
  JMultipleFileScanner<JACOUSTICS::JEvt>                             calibrationFile;
  JFileRecorder<JTypeList<JMeta, JRECONSTRUCTION::JAnalysisSummary>> outputFile;
  JLimit_t&                                                          numberOfEvents = inputFile.getLimit();
  double                                                             margin;
  string                                                             detectorFile;
  double                                                             Tmax_s;
  int                                                                debug;
 
  try {
 
 
    JParser<> zap("Program to produce analysis variables in offline files. The variable definition is given in JAnalysisSummary 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['m'] = make_field(margin)              = 20.0;
    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);
  }
 
  JCylinder3D can(detector.begin(), detector.end());
  can.addMargin(margin);
 
  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::JAnalysisSummary analysis;
 
    Evt* evt = in.next();
 
        
    for (vector<Trk>::iterator track = evt->mc_trks.begin(); track != evt->mc_trks.end(); ++track) {
      track->clearusr(); // Prevent visible energy from being calculated based on `mc_usr_keys::energy_lost_in_can`       
    }
 
    if (has_primary(*evt)) {
 
      analysis.visible_energy_total  = getVisibleEnergy             (*evt, can);
      analysis.visible_energy_lepton = (has_leading_lepton(*evt) ?
                                        getVisibleEnergyLeadingLepton(*evt, can) :
                                        0.0);
      analysis.visible_energy_hadron = analysis.visible_energy_total - analysis.visible_energy_lepton;
    }
    
    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() > analysis.nMaxPMT_snapshot) {
        
        analysis.nMaxPMT_snapshot = pmt_set.size();
        analysis.nMaxPMT_snapshot_ModuleID = dom_id;
 
      }
    }
 
    for (const auto& [dom_id, pmt_set] : trigger) {
      if (pmt_set.size() > analysis.nMaxPMT_trigger) {
        
        analysis.nMaxPMT_trigger = pmt_set.size();
        analysis.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 < analysis.Dmin) {
              
              analysis.Dmin          = D;
              analysis.Dmin_ModuleID = module.getID();
              
            }
          }
        }
      }
    } 
    outputFile.put(analysis);
  } 
 
  JSingleFileScanner<JMeta> io(inputFile);
 
  io >> outputFile;
  STATUS(endl);
  outputFile.close();
 
}