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/JDAQHitSelector.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"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
 
/**
 * \author mlincett,selhedri,iagoos
 */
 
namespace JSUPERNOVA {
 
  using namespace std;
  using namespace JPP;
  using namespace KM3NETDAQ;
 
  typedef vector<double> multiplicities_t;
 
  typedef set<int> JModuleSet;
 
  /**
   * Auxiliary class to store reduced information of a coincidence on an optical module
   * This class allows storing the observables associated with searches for low-energy
   * neutrinos from Core-Collapse supernovae:
   */
  class JCoincidenceSN {
 
  public:
    double time;
    double timeslice_time;
    int    moduleID;
    int vetoIndex;
    // Observables computed using hits in 10ns window in same DOM
    int    multiplicity;
    double mean_dir_ctheta = -2; // mean PMT direction cos(theta)
    double mean_dir_norm = -1; // hit concentration |R|
    double deltaT = -1; // mean time spread
    double total_ToT = -1; // total ToT
    bool   fake      = false; // For fire drills (when simulated events are injected)
 
    void clear(){
      this->vetoIndex = 0;
      this->multiplicity = 0;
      this->mean_dir_ctheta = -2;
      this->mean_dir_norm = -1;
      this->total_ToT = -1;
      this->deltaT = -1;
      this->fake = false;
    }
 
  /**
   * Constructor:
   *
   * \param    t                 Event timing within a timeslice
   * \param    m                 Number of hits within 10ns
   * \param    dom               ID of the module
   * \param    vetoIndex         Muon veto index: 0 if pass all vetoes, 1 if pass only trigger veto, 2 if pass only correlation veto, 3 if fails all vetoes
   * \param    ctheta            Mean z position of hit cluster on DOM
   * \param    rnorm             Hit concentration |R|
   * \param    deltaT            Average time spread of hit cluster
   * \param    total_ToT         Total ToT in DOM
   * \param    timeslice_time    Absolute time of the timeslice
   * \param    fake              fake
   */
    JCoincidenceSN(double t, int m, int dom, int vetoIndex=-1, double ctheta=-2, double rnorm=-1, double deltaT=-1, double total_ToT=-1, double timeslice_time=-1, bool fake=false)
      : time(t), timeslice_time(timeslice_time), moduleID(dom), vetoIndex(vetoIndex), multiplicity(m), mean_dir_ctheta(ctheta), mean_dir_norm(rnorm), deltaT(deltaT), total_ToT(total_ToT), fake(fake)
    { }
  
  /**
   * Compute single-DOM observables |R|, cos(theta), total ToT, and Delta(t) and set the corresponding class attributes
   *
   * \param    hits              Hit vector from timeslice
   * \param    det               Single-DOM detx file used to simulate CCSN signal
   */
    bool operator()(const vector<JHitR0>& hits, const JDetector& det){
        if (hits.size() == 0) return 1;
        JVector3D mean_pmt_dir;
        double first_hit_time = 0;
        this->total_ToT = 0;
        this->multiplicity = 0;
        this->deltaT = 0;
        int hitcount=0;
        for(auto &hit: hits){
            this->multiplicity++;
            this->total_ToT += hit.getToT();
            if (hitcount > 0) {
              this->deltaT += hit.getT() - first_hit_time;
            }
            else first_hit_time = hit.getT();
            int channel_id = hit.getPMT();
            JPMTAddress pmt_address(JModuleRouter(det).getAddress(101),channel_id);
            auto pmt = det.getPMT(pmt_address).getDirection();
            mean_pmt_dir += pmt.getDirection();
            hitcount++;
        }
        if (this->multiplicity >= 2) {
            this->mean_dir_norm = mean_pmt_dir.getLength()/this->multiplicity;
            this->mean_dir_ctheta = mean_pmt_dir.getZ()/this->multiplicity;
            this->deltaT /= (this->multiplicity-1);
        }
        return 0;
    }
 
    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;
    // JDAQHitSelector& hitSelector;// = JDAQHitDefaultSelector();
  
  public:
    int frameIndex;
    JDAQUTCExtended timeUTC;
 
    /** 
     * Default constructor
     * Configures the trigger with a time window and the pretrigger threshold.
     */
    JDataSN(double window, int threshold = 4) // JDAQHitSelector& selector = JDAQHitDefaultSelector()) 
      : TMax_ns(window), min_M(threshold)     // , hitSelector(selector)
    {  };
 
    /**
     * 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
     * \param det       single-DOM detx file used to simulate CCSN neutrinos (default: empty detector)
     */
    void operator() (const vector<JHitR0>& in, const int& moduleID, const JDetector& det = JDetector()) {
    
      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) {
            JCoincidenceSN coincidence(p->getT(), M, moduleID);
            if (det.size() > 0) {
              vector<JHitR0> coincidence_hits;
              coincidence_hits.assign(p,q); 
              coincidence(coincidence_hits, det);
            }
            this->push_back(coincidence);
          }
 
          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
     * \param selector      hit selector
     * \param det           single-DOM detx file used to simulate signal events
     */
    void operator() (const JDAQTimeslice* timeslice, const JModuleRouter& moduleRouter, const JDAQHitSelector& selector = JDAQHitDefaultSelector(), const JDetector& det = JDetector()) {
 
      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),
                                                               selector); 
 
        buffer.preprocess(JPreprocessor::join_t, match);
      
        JSuperFrame1D_t& data = JSuperFrame1D_t::multiplex(buffer);
 
        if (det.size() < 0)(*this)(data, moduleID);
        else (*this)(data, moduleID, det);
    
      }
 
      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