software/JReconstruction/JConvertEvt.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <algorithm>
#include <memory>
#include <set>
#include <map>

#include "TROOT.h"
#include "TFile.h"

#include "km3net-dataformat/offline/Head.hh"
#include "km3net-dataformat/offline/MultiHead.hh"
#include "km3net-dataformat/offline/Evt.hh"
#include "km3net-dataformat/offline/io_online.hh"
#include "km3net-dataformat/definitions/trkmembers.hh"

#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JPMTRouter.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JPMTParametersMap.hh"
#include "JDetector/JPMTAnalogueSignalProcessor.hh"

#include "JDynamics/JDynamics.hh"

#include "JCompass/JEvt.hh"
#include "JCompass/JSupport.hh"

#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JSupport.hh"

#include "JTrigger/JHit.hh"
#include "JTrigger/JHitToolkit.hh"

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

#include "JDAQ/JDAQEventIO.hh"
#include "JDAQ/JDAQSummarysliceIO.hh"

#include "JSupport/JParallelFileScanner.hh"
#include "JSupport/JTreeScanner.hh"
#include "JSupport/JFileRecorder.hh"
#include "JSupport/JMonteCarloFileSupportkit.hh"
#include "JSupport/JSupportToolkit.hh"
#include "JSupport/JSupport.hh"
#include "JSupport/JMeta.hh"

#include "JReconstruction/JEvt.hh"
#include "JReconstruction/JEvtToolkit.hh"

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


/**
 * \file
 * Auxiliary program to convert fit results to Evt format.
 *
 * \author mdejong
 */
int main(int argc, char **argv)
{
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  typedef JTYPELIST<JCOMPASS::JOrientation,
                    JACOUSTICS::JEvt>::typelist                    calibration_types;
  typedef JMultipleFileScanner<calibration_types>                  JCalibration_t;

  JMultipleFileScanner<> inputFile;
  JFileRecorder<JTYPELIST<JAAnetTypes_t, JDAQSummaryslice, JMeta>::typelist> outputFile;
  JLimit_t&              numberOfEvents = inputFile.getLimit();
  string                 detectorFile;
  JCalibration_t         calibrationFile;
  double                 Tmax_s;
  JPMTParametersMap      pmtParameters;
  JQualitySorter         qualitySorter;
  size_t                 numberOfFits;
  int                    debug;

  try {

    JParser<> zap("Auxiliary program to convert fit results to Evt format.\
                  \nThe option -L corresponds to the name of a shared library \
                  \nand function so to rearrange the order of fit results.");

    zap['f'] = make_field(inputFile);
    zap['o'] = make_field(outputFile)          = "aanet.root";
    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['P'] = make_field(pmtParameters,       "PMT simulation data (or corresponding file name)")                 = JPARSER::initialised();
    zap['L'] = make_field(qualitySorter)       = JPARSER::initialised();
    zap['N'] = make_field(numberOfFits)        = 0;
    zap['d'] = make_field(debug)               = 2;

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

  if (detectorFile == "" && !calibrationFile.empty()) {
    FATAL("Missing detector file." << endl);
  }

  typedef JParallelFileScanner< JTypeList<JDAQEvent, JFIT::JEvt> >  JParallelFileScanner_t;
  typedef JParallelFileScanner_t::multi_pointer_type                multi_pointer_type;

  JDetector detector;

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

  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 JPMTRouter    pmt_router(dynamics ? dynamics->getDetector() : detector);
  const JModuleRouter mod_router(dynamics ? dynamics->getDetector() : detector);


  outputFile.open();

  Head header;

  try {
    header = getHeader(inputFile);
  } catch(const exception& error) {}

  JAANET::JHead buffer(header);

  // DAQ livetime is not yet set for real data, for Monte Carlo, it is set in JTriggerEfficiency

  if (buffer.pull(&JAANET::JHead::DAQ) == buffer.end()) {

    buffer.DAQ.livetime_s = getLivetime(inputFile.begin(), inputFile.end());
    buffer.push(&JAANET::JHead::DAQ);

    copy(buffer, header);
  }

  if (detectorFile != "") {

    buffer.calibration.buffer = (dynamics ? calibration::dynamical() : calibration::statical());
    buffer.push(&JAANET::JHead::calibration);  

    copy(buffer, header);
  }

  
  outputFile.put(header);


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

  map<int, size_t> miss_mod;
  map<int, size_t> miss_pmt;

  counter_type number_of_events = 0;

  for (JMultipleFileScanner<>::const_iterator i = inputFile.begin(); i != inputFile.end(); ++i) {
    
    STATUS("Processing: " << *i << endl);
    
    JParallelFileScanner_t in(*i);
    JTreeScanner<Evt>      mc(*i);

    in.setLimit(inputFile.getLimit());

    Vec center(0,0,0);

    int mc_run_id = 0;
    
    try {

      const JAANET::JHead head(getHeader(*i));
    
      center    = get<Vec>(head);
      mc_run_id = head.start_run.run_id;
      
    } catch(const exception& error) {}

    
    while (in.hasNext()) {

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

      multi_pointer_type ps = in.next();

      JDAQEvent*  tev   = ps;
      JFIT::JEvt* evt   = ps;

      if (dynamics) {
        dynamics->update(*tev);                                                      // update detector
      }

      JFIT::JEvt::iterator __end = evt->end();
      
      if (numberOfFits > 0) {
        advance(__end = evt->begin(), min(numberOfFits, evt->size()));
      }

      if (qualitySorter.is_valid()) {
        partial_sort(evt->begin(), __end, evt->end(), qualitySorter);
      }
      
      Evt out;                                                                       // output event

      if (mc.getEntries() != 0) {

        Evt* event = mc.getEntry(tev->getCounter());

        if (event != NULL) {

          out = *event;                                                              // import Monte Carlo true data

          if (out.mc_run_id == 0) {
            out.mc_run_id = mc_run_id;
          }

          for (vector<Trk>::iterator i = out.mc_trks.begin(); i != out.mc_trks.end(); ++i) {
            i->pos += center;                                                        // offset true positions
          }
        }
      }
      
      read(out, *tev);                                                               // import DAQ data

      if (!pmt_router->empty()) {                                                    // import calibration data

        for (vector<Hit>::iterator i = out.mc_hits.begin(); i != out.mc_hits.end(); ++i) {

          if (pmt_router.hasPMT(i->pmt_id)) {

            const JPMTAddress    address = pmt_router.getAddress(i->pmt_id);
            const JPMTIdentifier id      = pmt_router.getIdentifier(address);
            const JPMT&          pmt     = pmt_router.getPMT(address);

            i->dom_id     = id.getID();
            i->channel_id = id.getPMTAddress();
            i->pos        = getPosition (pmt);
            i->dir        = getDirection(pmt);

          } else {

            miss_pmt[i->pmt_id] += 1;
          }
        }
      }

      if (!mod_router->empty()) {                                                    // import calibration data

        for (vector<Hit>::iterator i = out.hits.begin(); i != out.hits.end(); ++i) {

          if (mod_router.hasModule(i->dom_id)) {

            const JPMTIdentifier id(i->dom_id, i->channel_id);

            const JPMTParameters&             parameters = pmtParameters.getPMTParameters(id);
            const JPMTAnalogueSignalProcessor cpu(parameters);

            const JPMT& pmt = mod_router.getPMT(id);

            JTRIGGER::JHit::setSlewing(parameters.slewing);
            
            const JTRIGGER::JHit hit(getTime(i->tdc, pmt.getCalibration()), i->tot);
            
            i->pmt_id = pmt.getID();
            i->pos    = getPosition (pmt);
            i->dir    = getDirection(pmt);
            i->t      = hit.getT();
            i->tot    = hit.getToT();
            i->a      = cpu.getNPE(i->tot);

          } else {

            miss_mod[i->dom_id] += 1;
          }
        }
      }

      struct : public set<JUUID> {

        inline int get_index(const JUUID& element) const
        {
          const_iterator i = this->find(element);
          
          if (i != this->end())
            return std::distance(this->begin(), i);
          else
            return -1;
        }
      } uuid;
      
      for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {
        uuid.insert(fit->getUUID());
      }

      for (JFIT::JEvt::const_iterator fit = evt->begin(); fit != __end; ++fit) {     // import reconstruction data
        
        Trk trk;
        
        trk.id        = uuid.get_index(fit->getUUID());
        trk.pos       = Vec(fit->getX(),  fit->getY(),  fit->getZ());
        trk.dir       = Vec(fit->getDX(), fit->getDY(), fit->getDZ());
        trk.t         = fit->getT();
        trk.E         = fit->getE();
        trk.lik       = fit->getQ();
        trk.rec_type  = JPP_RECONSTRUCTION_TYPE;
        trk.status    = (fit->getStatus() == COMPLETE_CHAIN ? TRK_ST_FINALSTATE : TRK_ST_UNDEFINED);
        
        if (fit->hasParentUUID()) {
          trk.mother_id = uuid.get_index(fit->getParentUUID());
        }

        for (JHistory::const_iterator i = fit->getHistory().begin(); i != fit->getHistory().end(); ++i) {
          trk.rec_stages.push_back(i->type);
        }
        
        for (int i = 0; i != fit->getN(); ++i) {
          trk.fitinf.push_back(fit->getW(i));
        }

        trk.error_matrix = fit->getV();
        
        out.trks.push_back(trk);
      }

      out.id = ++number_of_events;

      outputFile.put(out);
    }
  }

  STATUS(endl);

  for (const auto& i : miss_pmt) { ERROR("Misses PMT    "  << setw(8) << i.first << ' ' << setw(8) << i.second << endl); }
  for (const auto& i : miss_mod) { ERROR("Misses module "  << setw(8) << i.first << ' ' << setw(8) << i.second << endl); }

  JMultipleFileScanner<JTYPELIST<MultiHead, JDAQSummaryslice, JMeta>::typelist> io(inputFile);

  io >> outputFile;

  outputFile.close();
}