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 "JGeometry3D/JAxis3D.hh"

#include "JTools/JRange.hh"

#include "JTrigger/JEventToolkit.hh"
#include "JTrigger/JHitR0.hh"
#include "JTrigger/JMatchL0.hh"
#include "JTrigger/JPreprocessor.hh"
#include "JTrigger/JSuperFrame1D.hh"
#include "JTrigger/JSuperFrame2D.hh"
#include "JTrigger/JTriggerParameters.hh"
#include "JSupport/JTriggerParametersSupportkit.hh"

#include "TH1D.h"

namespace JSUPERNOVA {

  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;

  /**
   * Get ordinal of a module according to increasing number of string and floor
   * 
   * \param location        module location
   * \param numberOfFloors  detector number of floors
   */
  inline int getModuleIndex(const JLocation& location, const int numberOfFloors) {
    return (location.getString() - 1)  * numberOfFloors + location.getFloor();
  }  

  typedef vector<double> multiplicities_t;

  typedef set<int> JModuleSet;

  /**
   * Get module indices from a set of modules
   * \param in             module set
   * \param router         module router
   * \param numberOfFloors number of floors per string
   */

  set<int> getModuleIndices(const JModuleSet& in, const JModuleRouter& router, int numberOfFloors = 18) {

    set<int> indices;

    for (JModuleSet::const_iterator p = in.begin(); p != in.end();  ++p) {
      indices.insert(getModuleIndex(router.getModule(*p).getLocation(), numberOfFloors));
    }

    return indices;
  }


  /**
   * 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 hit_type;
   
      for (JDAQEvent::const_iterator<hit_type> hit = event.begin<hit_type>();
           hit != event.end<hit_type>();
           ++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 min_M;
  
  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), min_M(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 >= min_M) {
            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 JHitR0 hit_t;

      typedef JSuperFrame2D<hit_t> JSuperFrame2D_t;

      typedef JSuperFrame1D<hit_t> JSuperFrame1D_t;

      const JMatchL0<hit_t>  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)); 

        buffer.preprocess(JPreprocessor::join_t, 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) const {
      return any_of(this->begin(), this->end(), JMatchVeto(in)); 
    }
    
  };


  /**
   * Auxiliary class to manage a cluster of coincidences
   */
  class JClusterSN : public vector<JCoincidenceSN> {

  public:

    /**
     * Finds coincidence with maximum multiplicity.
     */
    JCoincidenceSN getPeak() const {
      
      JClusterSN::const_iterator p = this->begin();
      
      for (JClusterSN::const_iterator q = p + 1; q != this->end(); q++) {
        if (q->getMultiplicity() > p->getMultiplicity()) {
          p = q;
        }
      }

      return (*p);

    }


    /* 
     * Returns 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;
    }  
    
  };


  /**
   * Interface for SN filter operator.
   * This operator is meant to determine the SN trigger count from a vector of JClusterSN
   */
  class JSNFilter {
  public:
    virtual bool operator() (const JCoincidenceSN& in) const = 0;
    virtual bool operator() (const JClusterSN&     in) const = 0;
  };


  /**
   * SN filter based on multiplicity selection
   * optional suppression of multi-module coincidences
   * WARNING: no minimum threshold for the veto
   */
  class JSNFilterM : public JSNFilter {
    
  private:
    JRange<int> A;
    bool mode;

  public:
    JSNFilterM(JRange<int> R, int m = 0)
      : A(R), mode(m)
    {}

    // select coincidence if within multiplicity acceptance
    bool operator() (const JCoincidenceSN& in) const {
      return A(in.getMultiplicity());
    } 
   
    // select cluster if any coincidence is accepted
    // optionally veto if more of one module is interested
    bool operator() (const JClusterSN& in) const {

      bool out = (*this)(in.getPeak());

      if (mode == 1) {
        out = out && (in.getModules().size() == 1);
      }

      return out;
    }  

  };

  
  /**
   * SN filter based on veto window
   */
  class JSNFilterMV : public JSNFilter {
    
  private:
    JRange<int> A;
    JVetoSet    V;

  public:
    JSNFilterMV(JRange<int>& R, JVetoSet& S)
      : A(R), V(S)
    {}

    bool operator() (const JCoincidenceSN& in) const {
      return A(in.getMultiplicity()) && !V(in);
    } 
   
    bool operator() (const JClusterSN& in) const {
      return (*this)(in.getPeak());
      // 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_ns = 1000.0; 

    JVetoSet veto;
  
  public:

    int frameIndex;
    JDAQUTCExtended timeUTC;

    /** 
     * Default constructor
     *
     * \param TMax_ns time width for coincidence clustering (track / afterpulse)
     */
    JTriggerSN(double TMax_ns) :
      TMaxCluster_ns(TMax_ns)
    {};

    void setVeto(JVetoSet &vt) {
      veto = vt;
    }


    /** 
     * Builds trigger by clustering pretrigger data
     *
     * \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_ns))) {
          cluster.push_back(*q);
          ++q; 
        } 

        p = q;

        this->push_back(cluster);      

      }
    }

    
    /**
     * 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) ) {

          // only clusters with one module pass the selection JSNFilterM(A, 1)
          triggeredModules.insert((*p)[0].getModule());

        }

      }
      
      return triggeredModules;

    }


    /**
     * Get triggered modules according to given filter
     */
    JModuleSet getModules(const JSNFilter& F) {
      
      JModuleSet out;

      for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {
        if ( F(*p) ) {
          out.insert(p->getPeak().getModule());
        }
      }
      return out;
    }

    /**
     * Fill histogram with multiplicity spectrum resulting from given filter
     */
    void fill(TH1D* out, const JSNFilter& F) {
      for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {       
        if (F(*p)) {
          out->Fill( p->getPeak().getMultiplicity() );
        }
      }
    }

    vector<double> getMultiplicities(const JSNFilter& F) {
      vector<double> out;

      for (JTriggerSN::const_iterator p = this->begin(); p != this->end(); p++) {       
        if (F(*p)) {
          out.push_back( p->getPeak().getMultiplicity() );
        }
      }
      return out;
    }
      



    /**
     * > 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 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