JAcoustics/JEvent.hh

Go to the documentation of this file.

#ifndef __JACOUSTICS__JEVENT__
#define __JACOUSTICS__JEVENT__

#include <string>
#include <istream>
#include <ostream>
#include <iomanip>
#include <vector>
#include <algorithm>

#include <TROOT.h>
#include <TObject.h>

#include "JAcoustics/JCounter.hh"
#include "JAcoustics/JTransmission.hh"

#include "JTools/JQuantile.hh"


/**
 * \file
 *
 * Acoustic event.
 * \author mdejong
 */
namespace JACOUSTICS {}
namespace JPP { using namespace JACOUSTICS; }

namespace JACOUSTICS {

  /**
   * Acoustic event.
   */
  struct JEvent :
    public JCounter,
    public std::vector<JTransmission>,
    public TObject
  {
    /**
     * Auxiliary class to determine value of acoustic events.\n
     * This class can be used with JSUPPORT::JTreeScanner so to read acoustics events in order of time-of-emission.
     */
    struct JEvaluator {
      /**
       * Default constructor.
       */
      JEvaluator()
      {}


      /**
       * Get value of object.
       *
       * \param  event           event
       * \return                 value
       */
      inline double operator()(const JEvent& event) const
      {
        return event.begin()->getToE();
      }
    };


    /**
     * Default constructor.
     */
    JEvent() :
      JCounter(),
      overlays(0),
      id      (-1)
    {}


    /**
     * Constructor.
     *
     * The transmissions will be sorted to ensure proper functioning of method JEvent::merge and class JEventOverlap.
     * 
     * \param  oid            detector identifier
     * \param  counter        counter
     * \param  id             identifier
     * \param  __begin        begin of data
     * \param  __end          end   of data
     */
    template<class T>
    JEvent(const std::string& oid,
           const int counter,
           const int id,
           T         __begin,
           T         __end) :   
      JCounter(counter),
      oid     (oid),
      overlays(0),
      id      (id)
    {
      using namespace std;

      for (T i = __begin; i != __end; ++i) {
        push_back(*i);
      }

      sort(this->begin(), this->end());
    }


    /**
     * Virtual destructor.
     */
    virtual ~JEvent()
    {}


    /**
     * Get detector identifier.
     *
     * \return                detector identifier.
     */
    const std::string& getOID() const
    {
      return oid;
    }


    /**
     * Get counter.
     *
     * \return                counter
     */
    int getCounter() const
    {
      return counter;
    }


    /**
     * Get overlays.
     *
     * \return                overlays
     */
    int getOverlays() const
    {
      return overlays;
    }


    /**
     * Get identifier.
     *
     * \return                identifier
     */
    int getID() const
    {
      return id;
    }


    /**
     * Merge event.
     *
     * It is assumed that the transmissions in both events are ordered
     * according the default less-than operator.
     *
     * \param  event          event
     */
    void merge(const JEvent& event)
    {
      using namespace std;

      vector<JTransmission> buffer;

      const_iterator __hit1 = this ->begin();
      const_iterator __end1 = this ->end();

      const_iterator __hit2 = event.begin();
      const_iterator __end2 = event.end();

      buffer.resize(this->size() + event.size());

      iterator out = buffer.begin();

      while (__hit1 != __end1 && __hit2 != __end2) {

        if        (*__hit1 < *__hit2) {

          *out = *__hit1;
          ++__hit1;

        } else if (*__hit2 < *__hit1) {

          *out = *__hit2;
          ++__hit2;

        } else {

          *out  = *__hit1;

          ++__hit1;
          ++__hit2;
        }

        ++out;
      }

      // append remaining hits from either set

      out = copy(__hit1, __end1, out);
      out = copy(__hit2, __end2, out);

      buffer.resize(distance(buffer.begin(), out));

      this->swap(buffer);

      ++overlays;
    }


    /**
     * Remove outliers.
     *
     * First, the transmissions are sorted according to their time-of-emission.
     * The mean and standard deviation of the times of emission is then determined,
     * excluding the first and last possible number of outliers from either side.
     * From the excluded transmissions, the ones with a time-of-emission less than 
     * the given number of standard deviations away from the mean are retained.
     *
     * \param  ns             number of outliers
     * \param  stdev          number of standard deviations
     */
    void remove(const size_t ns, const double stdev)
    {
      using namespace std;
      using namespace JPP;

      if (ns > 0 && 2*ns < this->size()) {

        sort(this->begin(), this->end());

        iterator         p = this->begin();
        reverse_iterator q = this->rbegin();

        advance(p, ns);
        advance(q, ns);

        JQuantile Q;

        for (const_iterator i = p; i != q.base(); ++i) {
          Q.put(i->getToE());
        }

        const double t1 = Q.getMean()  -  stdev * Q.getSTDev();
        const double t2 = Q.getMean()  +  stdev * Q.getSTDev();

        for (size_t i = ns; i != 0 && p->getToE() >= t1; --i, --p) {}
        for (size_t i = ns; i != 0 && q->getToE() <= t2; --i, --q) {}

        vector<JTransmission> buffer(p, q.base());

        this->swap(buffer);
      }
    }


    /**
     * Write event to output stream.
     *
     * \param  out            output stream
     * \param  event          event
     * \return                output stream
     */
    friend inline std::ostream& operator<<(std::ostream& out, const JEvent& event)
    {
      using namespace std;

      out << event.getOID()                 << endl;
      out << setw(8) << event.getCounter()  << endl;
      out << setw(2) << event.getOverlays() << endl;
      out << setw(3) << event.getID()       << endl;

      for (const_iterator i = event.begin(); i != event.end(); ++i) {

        out << setw(10) << i->getID()        << ' ' 
            << setw(10) << i->getRunNumber() << ' ' 
            << fixed << setw(12) << setprecision(6) << i->getToA()  << ' '
            << fixed << setw(12) << setprecision(6) << i->getToE()  << ' '
            << fixed << setw(8)  << setprecision(0) << i->getQ()    << endl;
      }

      return out;
    }

    ClassDef(JEvent, 2);
    
  protected:
    std::string oid;
    int         overlays;
    int         id;
  };


  /**
   * Less than operator for acoustics events.
   *
   * The less than operator is applied to the first hit in the events.\n
   * If there are no hits in either event, the counter of the events is used.
   *
   * \param  first          first  event
   * \param  second         second event
   * \return                true if first event earliear than second; else false
   */
  inline bool operator<(const JEvent& first, const JEvent& second)
  {
    if (!first.empty() && !second.empty())
      return first.begin()->getToE() < second.begin()->getToE();
    else
      return first.getCounter() < second.getCounter();
  }
}

#endif