JRandomTimeslice.hh

Go to the documentation of this file.

#ifndef __JTIMESLICE__JRANDOMTIMESLICE__
#define __JTIMESLICE__JRANDOMTIMESLICE__

#include <vector>

#include <TRandom3.h>

#include "JDetector/JDetectorSimulator.hh"
#include "JDetector/JPMTDefaultSimulatorInterface.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]
        
        JModuleData buffer;
        
        for (JDetector::const_iterator module = simbad->begin(); module != simbad->end(); ++module) {

          if (!module->empty()) {

            buffer.reset(module->size());

            simbad.generateHits(*module, period, buffer);

            this->push_back(JDAQSuperFrame(JDAQSuperFrameHeader(chronometer, module->getID())));

            simbad(*module, buffer, *(this->rbegin()));
          }
        }
      }
    }


    /**
     * 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.
     *
     * 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 JDAQHit::JTDC_t  JTDC_t;


      size_t N;                                                  // number of time intervals

      if      (T_ns < 1.0)
        N = (size_t)  getFrameTime();
      else if (T_ns < getFrameTime())
        N = (size_t) (getFrameTime() / T_ns + 0.5);
      else
        N = 1;

      const JTDC_t T_step = (JTDC_t) (getFrameTime() / N);       // TDC interval


      typedef vector<JDAQHit>  buffer_type;

      vector<size_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) {
            index [i] = i;
            buffer[i].clear();
          }


          vector<size_t> N_max(NUMBER_OF_PMTS, N);               // number of time intervals per PMT
          vector<JTDC_t> T_max(NUMBER_OF_PMTS, 0);               // 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();
          }

          for (int i = 0; i != NUMBER_OF_PMTS; ++i) {

           const size_t N_i = (T_max[i] + T_step) / T_step;

            if (N_i < N_max[i]) {
              N_max[i] = N_i;
            }
          }

          sort(N_max.begin(), N_max.end());                      // swaps by allowed time range

          for (size_t i = 0, L = 0; i != N_max.size(); ++i) {

            size_t M = N_max[i];

            if (L != M) {
              for (size_t i = M - L - 1; i != 0; --i) {
                std::swap(index[L + i], index[L + gRandom->Integer(i+1)]);
              }
            }

            L = M;
          }

          const size_t M = N_max[NUMBER_OF_PMTS - 1];


          // Wall street shuflle

          JDAQSuperFrame::const_iterator hit = frame->begin();

          for (size_t in = 0; in != M; ++in) {

            const size_t  out   = index[in];
            buffer_type&  zbuf  = buffer[out];

            const JTDC_t  T_in  = in  * T_step;
            const JTDC_t  T_out = out * T_step;

            for ( ; hit != frame->end() && hit->getT() < T_in + T_step; ++hit) {
              zbuf.push_back(JDAQHit(hit->getPMT(), (hit->getT() - T_in) + T_out, hit->getToT()));
            }
          }


          // copy back data

          for (size_t i = 0, number_of_hits = 0; i != M; ++i) {

            memcpy(frame->data() + number_of_hits, buffer[i].data(), buffer[i].size() * sizeof(JDAQHit)); 

            number_of_hits += buffer[i].size();
          }
        }
      }
    }
  };
}

#endif