JSupernova.hh

Go to the documentation of this file.

#ifndef __JSUPERNOVA_JSUPERNOVA__
#define __JSUPERNOVA_JSUPERNOVA__
 
#include <set>
 
#include "km3net-dataformat/online/JDAQEvent.hh"
#include "km3net-dataformat/online/JDAQTriggeredHit.hh"
#include "km3net-dataformat/online/JDAQTimeslice.hh"
#include "km3net-dataformat/online/JDAQUTCExtended.hh"
 
#include "JDetector/JDAQHitRouter.hh"
#include "JDetector/JModuleRouter.hh"
 
#include "JTools/JRange.hh"
 
#include "JTrigger/JEventToolkit.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
 
#include "TH1D.h"
 
namespace JSUPERNOVA {
 
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef set<int> JModuleSet;
 
  /**
   * Auxiliary class to store reduced information of a coincidence on an optical module
   */
  class JCoincidenceSN {
 
  private:
    double time;
    int    multiplicity;
    int    moduleID;
 
  public:
    JCoincidenceSN(double t, int m, int dom)
      : time(t), multiplicity(m), moduleID(dom)
    { }
  
    int getMultiplicity() const {
      return multiplicity;
    }
 
    int getModule() const {
      return moduleID;
    }
 
    double getTime() const {
      return time;
    }
 
    bool operator<(const JCoincidenceSN& rhs) const {
      return (time < rhs.time);
    }
 
  };
 
  /**
   * Auxiliary class to define a veto time window on a set of optical modules
   */
  class JVeto {
 
  private:
    JTimeRange timeRange;
    JModuleSet moduleSet;
 
  public:  
 
    /**
     * Default constructor
     * \param event     DAQ event
     * \param hitRouter hit router as source of hit time calibration
     * 
     */
    JVeto(const JDAQEvent& event, const JDAQHitRouter& hitRouter) {
 
      timeRange = JTimeRange::DEFAULT_RANGE;
 
      typedef JDAQTriggeredHit JHit_t;
   
      for (JDAQEvent::const_iterator<JHit_t> hit = event.begin<JHit_t>();
           hit != event.end<JHit_t>();
           ++hit) {
     
        moduleSet.insert(hit->getModuleID());
    
        timeRange.include(getTime(*hit, hitRouter.getPMT(*hit)));
   
      } 
    }
 
    
    /**
     * Get length of veto time range
     */
    double getLength() {
      return timeRange.getLength();
    }  
 
    /**
     * Determines if a coincidence is vetoed
     * \param in coincidence under test
     */
    bool operator() (const JCoincidenceSN& in) const {
      return timeRange(in.getTime()) && 
        (moduleSet.find(in.getModule()) != moduleSet.end());
    }
 
  };
 
 
  /**
   * Auxiliary class to build the supernova trigger dataset
   */
  class JDataSN
    : public vector<JCoincidenceSN> {
 
  private:
    int TMax_ns;
    int multiplicityThreshold;
  
  public:
    int frameIndex;
    JDAQUTCExtended timeUTC;
 
    /** 
     * Default constructor
     * Configures the trigger with a time window and the pretrigger threshold.
     */
    JDataSN(double window, int threshold = 4) 
      : TMax_ns(window), multiplicityThreshold(threshold)
    {  };
 
    /**
     * Builds coincidences from calibrated hit data and loads them in the internal vector.
     * \param in        hit data
     * \param moduleID  optical module ID for the hit data
     */
    void operator() (vector<JHitR0> in, int moduleID) {
    
      if (in.size() > 1) {
 
        for (vector<JHitR0>::const_iterator p = in.begin(); p != in.end(); ) {
 
          vector<JHitR0>::const_iterator q = p;
 
          while (++q != in.end() && q->getT() - p->getT() <= TMax_ns ) {}
 
          int M = distance(p, q);
 
          if (M >= multiplicityThreshold)
            { this->push_back(JCoincidenceSN(p->getT(), M, moduleID)); }
 
          p = q;
          
        } 
      }
    }
 
 
    /**
     * Builds coincidences from a timeslice and loads them into the internal vector.
     * Double pulses are merged according to the coincidence time window.
     * 
     * \param timeslice     input timeslice
     * \param moduleRouter  detector module router
     */
    void operator() (const JDAQTimeslice* timeslice, const JModuleRouter& moduleRouter) {
 
      frameIndex = timeslice->getFrameIndex();
      timeUTC = timeslice->getTimesliceStart();
 
      typedef JSuperFrame2D<JHitR0> JSuperFrame2D_t;
 
      typedef JSuperFrame1D<JHitR0> JSuperFrame1D_t;
 
      const JMatchL0<JHitR0> match(TMax_ns);
 
      for (JDAQTimeslice::const_iterator frame = timeslice->begin(); frame != timeslice->end(); ++frame) {
 
        int moduleID = frame->getModuleID();
 
        JSuperFrame2D_t& buffer = JSuperFrame2D_t::demultiplex(*frame, moduleRouter.getModule(moduleID)); 
 
        for (JSuperFrame2D_t::iterator i = buffer.begin(); i != buffer.end(); ++i) { i->join(match); }
      
        JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer); data.pop_back();
 
        (*this)(data, moduleID);
    
      }
 
      sort(this->begin(), this->end());
 
    }
    
  };
 
 
  /**
   * Auxiliary class-operator to match a JVeto with a JCoincidenceSN object
   * Provides an operator to test if a coincidence is vetoed
   */
 
  class JMatchVeto {
 
  private:
    JCoincidenceSN dut;
  
  public:
    /**
     * Default constructor
     * \param in    coincidence to be matched against veto
     */
    JMatchVeto(const JCoincidenceSN& in) : dut(in) {}
 
    /**
     * Operator
     * \param in    veto to be matched against inner coincidence
     */
    bool operator() (const JVeto& in) {
      return in(dut);
    }
 
  };
 
 
  /**
   * Auxiliary class to manage a set of vetoes
   */
  class JVetoSet : public vector<JVeto> {
    
  public:
    /** 
     * Applies the veto set to a coincidence
     * \param in    coincidence to be tested
     */ 
    bool operator() (const JCoincidenceSN& in) {
      return any_of(this->begin(), this->end(), JMatchVeto(in)); 
    } 
    
  };
 
 
  /**
   * Auxiliary class to manage a cluster of coincidences
   */
  class JClusterSN : public vector<JCoincidenceSN> {
 
  public:
    /* 
     * Return the set of modules spanned over by the cluster
     */
    JModuleSet getModules() const {
 
      JModuleSet out;
      
      for (JClusterSN::const_iterator p = this->begin(); p != this->end(); p++) {
        out.insert(p->getModule());
      }
 
      return out;
    }
 
    
    
    /* 
     * Return the set of modules spanned over by the cluster
     * \param multiplicityThreshold minimum multiplicity to count the coincidence
     */
    JModuleSet getModules(int multiplicityThreshold) const {
 
      JModuleSet out;
      
      for (JClusterSN::const_iterator p = this->begin(); p != this->end(); p++) {
        if (p->getMultiplicity() >= multiplicityThreshold) {
          out.insert(p->getModule());
        }
      }
 
      return out;
    }
   
    
  };
 
 
  /**
   * Interface for SN filter operator.
   * This operator is meant to determine the SN trigger count.
   */
  class JSNFilter {
  public:
    virtual bool operator() (const JCoincidenceSN& in) = 0;
    virtual bool operator() (const JClusterSN&     in) = 0;
  };
 
 
  /**
   * SN filter based on multiplicity selection
   * optional suppression of multi-module coincidences
   */
  class JSNFilterM : JSNFilter {
    
  private:
    JRange<int> A;
    bool mode;
 
  public:
    JSNFilterM(JRange<int> R, int m = 0) : A(R), mode(m) {}
 
    bool operator() (const JCoincidenceSN& in) {
      return A(in.getMultiplicity());
    } 
   
    bool operator() (const JClusterSN& in) {
      bool out = any_of(in.begin(), in.end(), *this);
 
      if (mode == 1) {
        out = out && (in.getModules().size() == 1);
      }
 
      return out;
    }  
 
  };
 
  
  /**
   * Select clusters without correlated coincidences
   */
  class JSNFilterNM : JSNFilter {
    int numberOfModules;
    int multiplicityThreshold;
 
  public:
    JSNFilterNM(const int n, const int m) :
      numberOfModules(n), multiplicityThreshold(m)
    {}
 
    bool operator() (const JCoincidenceSN& in) {
      return true;
    }
 
    bool operator() (const JClusterSN& in) {
      JModuleSet modules = in.getModules(multiplicityThreshold);
      return (modules.size() <= 1);
    }
 
  };
 
  
  /**
   * SN filter based on veto window
   */
  class JSNFilterMV : JSNFilter {
    
  private:
    JRange<int> A;
    JVetoSet    V;
 
  public:
    JSNFilterMV(JRange<int>& R, JVetoSet& S) : A(R), V(S) {}
 
    bool operator() (const JCoincidenceSN& in) {
      return A(in.getMultiplicity()) && !V(in);
    } 
   
    bool operator() (const JClusterSN& in) {
      return any_of(in.begin(), in.end(), *this);
    }  
 
  };
 
 
  /**
   * Auxiliary class to apply the supernova trigger to SN data
   */
  class JTriggerSN : public vector<JClusterSN> {
    
  private:
 
    double TMaxCluster_ms = 1000.0; 
 
    JVetoSet veto;
  
  public:
 
    int frameIndex;
    JDAQUTCExtended timeUTC;
 
    /** 
     * Default constructor
     *
     * \param tRange  multiplicity range for selection of supernova coincidences
     * \param TMax_ms time width for coincidence clustering (track / afterpulse)
     */
    JTriggerSN(JRange<int> tRange, double TMax_ms) :
      TMaxCluster_ms(TMax_ms)
    {};
 
    void setVeto(JVetoSet &vt) {
      veto = vt;
    }
 
 
    /** 
     * Builds trigger
     *
     * \param in pretrigger SN data
     */
    void operator() (const JDataSN& in) {
    
      frameIndex = in.frameIndex;
      timeUTC    = in.timeUTC;
 
      for (vector<JCoincidenceSN>::const_iterator p = in.begin(); p != in.end(); ) {       
      
        JClusterSN cluster;
 
        cluster.push_back(*p);
 
        vector<JCoincidenceSN>::const_iterator q = p + 1;
 
        while(q != in.end() && (q->getTime() <= (p->getTime() + TMaxCluster_ms))) {
          cluster.push_back(*q);
          ++q; 
        } 
 
        p = q;
 
        this->push_back(cluster);      
 
      }
    }
 
 
    /**
     * Benchmark different trigger steps
     *
     * 0: count every triggering cluster of coincidences
     * 1: track veto based on time-correlated coincidences
     * 2: track veto based on triggered events
     * \return vector with count as a function of trigger step
     */
    vector<int> benchmark(JRange<int> A = JRange<int>(6,10)) {
 
      const int nBenchmarks = 3;
 
      vector<int> bench(nBenchmarks, 0);
 
      bench[0] = count_if(this->begin(), this->end(), JSNFilterM(A, 0));
      bench[1] = count_if(this->begin(), this->end(), JSNFilterM(A, 1));
      bench[2] = count_if(this->begin(), this->end(), JSNFilterMV(A, veto));
 
      return bench;
 
    }
 
    
    /**
     * Get triggered modules after track veto
     *
     * \return std::set containing triggered modules IDs
     */
    JModuleSet getModules(JRange<int> A = JRange<int>(6,10)) { 
 
      JModuleSet triggeredModules;
 
      JSNFilterM F(A, 1);
      
      for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
 
        if ( F(*p) ) {
          triggeredModules.insert((*p)[0].getModule());
        }
 
      }
      
      return triggeredModules;
 
    }
 
 
    /**
     * Get modules according to a filter criterion
    JModuleSet getModules(const JSNFilter& F) {
      JModuleSet out;
      
      for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
        if (F(*p)) {
          
        }
      }
 
    }*/
 
    
    /** 
     * Fills histogram with multiplicity counts
     */
    void fill(TH1D* out, const JSNFilter& F) {
      for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
 
        JSNFilterNM F(2, 4);
        
 
        for (JClusterSN::const_iterator q = p->begin(); q != p->end(); q++) {
          
        }
      }
    }
 
    /**
     * > operator
     * used by (reverse) std:priority_queue, in absence of this operator the container will behave unpredictably
     */
    bool operator>(const JTriggerSN& rhs) const {
      return (frameIndex > rhs.frameIndex);
    }
    
    /**
     * < operator
     */
    bool operator<(const JTriggerSN& rhs) const {
      return (frameIndex < rhs.frameIndex);
    }
 
 
  };
 
 
  /**
   * SN trigger summary information
   *
   */
  class JSummarySN {
  private:
    int    DETID;
    int    activeModules;
    int    RUNNR;
    int    frame;
    int    trigger;
    JDAQUTCExtended time;
    
  public:
    JSummarySN(int d, int am, int r, int f, JDAQUTCExtended tm, int tr)
      : DETID(d), activeModules(am), RUNNR(r), frame(f), trigger(tr), time(tm)
    { }
 
    friend inline std::ostream& operator<<(std::ostream& out, const JSummarySN& summary) {
      out << summary.DETID         << " ";
      out << summary.activeModules << " ";
      out << summary.RUNNR         << " ";
      out << summary.frame         << " ";
      out << summary.time          << " ";
      out << summary.trigger;
 
      return out;
    }
 
    friend inline std::istream& operator>>(std::istream& in, JSummarySN& summary) {
      in >> summary.DETID; 
      in >> summary.activeModules;
      in >> summary.RUNNR; 
      in >> summary.frame;
      in >> summary.time ;
      in >> summary.trigger;
 
      return in;
    }
 
  };
 
  /**
   * SN trigger statistics, the information is stored in the form of a count as a function of the trigger level.
   * the livetime needs to be set manually to compute the rates for the printing.
   */
 
  class JTriggerSNStats : public vector<double> {
 
  private:
    double livetime;
 
  public:
    /**
     * default constructor
     *
     * \param nModules number of modules of the detector
     */
    JTriggerSNStats(const int nModules) : vector<double>(nModules) {}
 
    void setLiveTime(const double lt) {
      livetime = lt;
    }
 
    /**
     * put statistics into printable form
     * outputs trigger level - rate - error
     */
    string toString() {
 
      ostringstream os;
 
      os << "=== SUPERNOVA TRIGGER STATISTICS ==" << endl;
      os << "=> livetime [s] = " << livetime << endl;
      os << "Level" << '\t' << "Rate (error)" << endl;
 
      if (livetime > 0) {
        for (unsigned i = 0; i < this->size(); i++) {
 
          double count = (*this)[i];
          double rate  = count / livetime;
          double error = 100 * (sqrt(count) / count);
 
          if (rate > 0) {
            os << i << '\t';
            os << scientific << setprecision(2) << rate << "Hz "; 
            os << fixed      << setprecision(0) << "(" << setw(2) << error << "%)";
            os << endl;
          }
        }
      }
 
      return os.str();
 
    }
 
      
    /**
     * put statistics into printable form
     * outputs trigger level - rate - error
     */
    string toSummaryFile() {
 
      ostringstream os;
 
      if (livetime > 0) {
 
        os << livetime << endl;
 
        for (unsigned i = 0; i < this->size(); i++) {
 
          double count = (*this)[i];
          double rate  = count / livetime;
          double error = (sqrt(count) / count);
 
          if (rate > 0) {
            os << i << ",";
            os << scientific;
            os << setprecision(2);
            os << rate << ",";
            os << error;
            os << endl;
          }
        }
      }
      
 
      return os.str();
      
    }
    
  };
 
}
 
#endif