JRandomTimeslice.hh

Go to the documentation of this file.

#ifndef __JTIMESLICE__JRANDOMTIMESLICE__
#define __JTIMESLICE__JRANDOMTIMESLICE__
 
#include <vector>
#include <set>
 
#include <TRandom3.h>
 
#include "JDetector/JDetectorSimulator.hh"
#include "JDetector/JPMTDefaultSimulatorInterface.hh"
#include "JDetector/JCLBDefaultSimulatorInterface.hh"
#include "JDetector/JPMTAddress.hh"
#include "JDetector/JTimeRange.hh"
#include "km3net-dataformat/online/JDAQ.hh"
#include "km3net-dataformat/online/JDAQClock.hh"
#include "JTimeslice/JTimesliceL0.hh"
 
 
/**
 * \author mdejong
 */
 
namespace KM3NETDAQ {
 
  using JDETECTOR::JDetectorSimulator;
 
 
  /**
   * Timeslice with random data.
   */
  struct JRandomTimeslice :
    public JTimesliceL0
  {
    /**
     * Default constructor.
     */
    JRandomTimeslice()
    {}
 
 
    /**
     * Constructor.
     *
     * \param  chronometer        chronometer
     * \param  simbad             detector simulator
     */
    JRandomTimeslice(const JDAQChronometer&    chronometer,
                     const JDetectorSimulator& simbad)
    {
      using namespace JPP;
 
      setDAQChronometer(chronometer);
 
      if (simbad.hasK40Simulator() &&
          simbad.hasPMTSimulator() &&
          simbad.hasCLBSimulator()) {
 
        const double     Tmin = getTimeSinceRTS(chronometer.getFrameIndex());   // [ns]
 
        const JTimeRange period(Tmin, Tmin + getFrameTime());                   // [ns]
 
        // generate hits per module
 
        vector<JModuleData> buffer(simbad->size());
 
        for (size_t i = 0; i != simbad->size(); ++i) {
 
          const JModule& module = (*simbad)[i];
 
          if (!module.empty() && simbad.getCLBSimulator().hasCLB(module.getID())) {
 
            buffer[i].resize(module.size());
 
            simbad.generateHits(module, period, buffer[i]);
          }
        }
 
 
        // generate mixed L0/L1 hits per unique module pair
 
        for (size_t m1 = 0; m1 != simbad->size(); ++m1) {
 
          for (size_t m2 = 0; m2 != m1; ++m2) {
 
            const JModule& M1 = (*simbad)[m1];
            const JModule& M2 = (*simbad)[m2];
 
            if (!M1.empty() && !M2.empty()) {
 
              if (neighbours(M1, M2)) {
 
                buffer[m1].resize(M1.size());
                buffer[m2].resize(M2.size());
 
                simbad.generateHits({M1,M2}, period, {buffer[m1],buffer[m2]});
              }
            }
          }
        }
 
 
        // store data
 
        for (size_t i = 0; i != simbad->size(); ++i) {
 
          if (!(*simbad)[i].empty()) {
 
            this->push_back(JDAQSuperFrame(JDAQSuperFrameHeader(chronometer, (*simbad)[i].getID())));
 
            simbad((*simbad)[i], buffer[i], this->back());
          }
        }
      }
    }
 
 
    /**
     * Recycle time slice by randomly shuffling time intervals of data.
     *
     * Hits within one time interval are swapped with hits within another -randomly selected- time interval.\n
     * Time intervals in which a high-rate veto occurred are excluded.
     *
     * Note that the time interval should be:
     *   - (much) larger  than time differences of hits in L1 coincidence; and
     *   - (much) smaller than time covered by data (as per KM3NETDAQ::getFrameTime()).
     *
     * \param  T_ns               time interval
     */
    void recycle(const double T_ns)
    {
      using namespace std;
      using namespace JPP;
 
      typedef vector<JDAQHit>  buffer_type;
      typedef JDAQHit::JTDC_t  JTDC_t;
 
      size_t N = (size_t) (getFrameTime() / T_ns + 0.5);         // number of time intervals
 
      if (N <   100) { N =   100; }                              // allow to keep indices at time of high-rate veto
      if (N > 50000) { N = 50000; }                              // maximum expected number of hits due to high-rate veto
 
      const JTDC_t Ts = (JTDC_t) (getFrameTime() / N);           // TDC interval
 
      random_indices_t    index (N);                             // indices of time intervals for swapping
      vector<buffer_type> buffer(N);                             // data per time interval
 
      for (iterator frame = this->begin(); frame != this->end(); ++frame) {
 
        if (!frame->empty()) {
 
          for (size_t i = 0; i != N; ++i) {
            buffer[i].clear();
          }
 
          // store data per time interval
          
          vector<JTDC_t> T_max(NUMBER_OF_PMTS,
                               numeric_limits<JTDC_t>::max());   // maximal time per PMT, e.g. due to high-rate veto
 
          for (JDAQSuperFrame::const_iterator hit = frame->begin(); hit != frame->end(); ++hit) {
 
            T_max[hit->getPMT()] = hit->getT();
 
            const int i = hit->getT() / Ts;
              
            buffer[i].push_back(*hit);
          }
 
          set<size_t> keep;                                      // indices corresponding to times of high-rate veto
          
          for (int pmt = 0; pmt != NUMBER_OF_PMTS; ++pmt) {
            if (frame->testHighRateVeto(pmt)) {
              keep.insert(T_max[pmt] / Ts);
            }
          }
 
          index.random_shuffle(keep);                            // randomly shuffle values between kept indices
 
          // Wall street shuflle
 
          JDAQSuperFrame::iterator hit = frame->begin();
 
          for (size_t in = 0; in != N; ++in) {
 
            const size_t  out   = index [in];
            buffer_type&  zbuf  = buffer[out];
 
            const JTDC_t  T_in  = in  * Ts;
            const JTDC_t  T_out = out * Ts;
 
            for (buffer_type::iterator i = zbuf.begin(); i != zbuf.end(); ++i, ++hit) {
              *hit = JDAQHit(i->getPMT(), (i->getT() - T_out) + T_in, i->getToT());
            }
          }
        }
      }
    }
 
    
    /**
     * Auxiliary data structure for randomisation of indices.
     */   
    struct random_indices_t :
      public std::vector<size_t>
    {
      /**
       * Constructor.
       *
       * \param  N             number of indices
       */
      random_indices_t(const size_t N) :
        std::vector<size_t>(N)
      {}
 
      /**
       * Randomly shuffle values between fixed indices.
       *
       * \param  keep          fixed indices
       */
      inline void random_shuffle(const std::set<size_t>& keep)
      {
        for (size_t i = 0; i != this->size(); ++i) {
          (*this)[i] = i;
        }
 
        size_t i1 = 0;
 
        for (const size_t i2 : keep) {
 
          random_shuffle(i1, i2);
            
          i1 = i2 + 1;
        }
          
        random_shuffle(i1, this->size());
      }
 
    private:
      /**
       * Randomly shuffle values between given indices.
       *
       * \param  i1            first index (included)
       * \param  i2            last  index (excluded)
       */
      inline void random_shuffle(const int i1, const int i2)
      {
        for (int i = i2 - 1; i > i1; --i) {
 
          const int l = i1 + gRandom->Integer(i - i1);
 
          std::swap((*this)[i], (*this)[l]);
        }
      }
    };
  };
}
 
#endif