JSirene.hh

Go to the documentation of this file.

#ifndef __JSIRENE__JSIRENE__
#define __JSIRENE__JSIRENE__

#include <limits>
#include <vector>
#include <algorithm>

#include "km3net-dataformat/offline/Hit.hh" 
#include "km3net-dataformat/offline/Trk.hh" 


/**
 * \author mdejong
 */

namespace JSIRENE {}
namespace JPP { using namespace JSIRENE; }

namespace JSIRENE {

  /**
   * Detector simulation parameters.
   */
  struct JParameters {
    /**
     * Default constructor.
     */
    JParameters()
    {
      Ecut_GeV  =   0.1;
      Emin_GeV  =   1.0;
      Dmin_m    =   0.1;
      Emax_GeV  = 250.0;
      Dmax_m    =  10.0;
      Tmax_ns   =   0.1;
      Nmax_npe  =   std::numeric_limits<int>::max();
    }

    double  Ecut_GeV;              //!< minimal energy for generation of light from shower [GeV]
    double  Emin_GeV;              //!< minimal energy of muon for shower generation [GeV]
    double  Dmin_m;                //!< minimal distance for positioning [m]
    double  Emax_GeV;              //!< maximal energy of muon below which step size is limited [GeV]
    double  Dmax_m;                //!< maximal step size when limited [m]
    double  Tmax_ns;               //!< maximal time between hits on same PMT to be merged
    int     Nmax_npe;              //!< maximal number of photo-electrons
  };


  /**
   * Auxiliary class to set-up Hit.
   *
   * This class is primarily used to limit the size of a Monte Carlo hit and 
   * thereby the memory usage of applications in case of large numbers of hits. 
   */
  struct JHit_t 
  {
    /**
     * Constructor.
     *
     * \param  id       identifier
     * \param  pmt_id   PMT identifier
     * \param  type     type
     * \param  origin   origin
     * \param  t        time [ns]
     * \param  npe      number of photo-electrons
     */
    JHit_t(const int    id,
           const int    pmt_id,
           const int    type,
           const int    origin,
           const double t,
           const int    npe)
    {
      this->id     = id;
      this->pmt_id = pmt_id;
      this->type   = type;
      this->origin = origin;
      this->t      = t;
      this->npe    = npe;
    }


    /**
     * Type conversion operator.
     *
     * \return          hit
     */
    operator const Hit& () const
    {
      static Hit hit;

      hit.id     = this->id;
      hit.pmt_id = this->pmt_id;
      hit.type   = this->type;
      hit.origin = this->origin;
      hit.t      = this->t;
      hit.a      = this->npe;

      return hit;
    }


    /**
     * Less than operator for hits.
     *
     * First hit is defined as:
     *    -# smallest PMT identifier;
     *    -# earliest time if same PMT identifier;
     *
     * \param  first            first  hit
     * \param  second           second hit
     * \return                  true if first hit earlier than second hit; else false
     */
    friend inline bool operator<(const JHit_t& first, const JHit_t& second)
    {
      if (first.pmt_id == second.pmt_id)
        return first.t      < second.t;
      else
        return first.pmt_id < second.pmt_id;
    }
    
    
    int    id;
    int    pmt_id;
    int    type;
    int    origin;
    double t;
    int    npe;
  };


  /**
   * Auxiliary data structure for list of hits with hit merging capability.
   */
  struct JHits_t :
    public std::vector<JHit_t>
  {
    /**
     * Merge hits on same PMT that are within given time window.
     *
     * The earliest hit has the sum of the number of photo-electrons of all following hits within given time window.\n
     * The hit identifiers are subsequently set in ascending order, starting at one.
     *
     * \param  Tmax_ns             maximal time difference [ns]
     */
    void merge(const double Tmax_ns)
    {
      using namespace std;

      if (!this->empty()) {

        sort(this->begin(), this->end());
      
        iterator in  = this->begin();
        iterator out = this->begin();
      
        out->id = 1;                              // set first hit identifier
      
        while (++in != this->end()) {
        
          if (out->pmt_id == in->pmt_id && in->t - out->t <= Tmax_ns) {

            out->npe += in->npe;                  // accumulate number of photo-electrons
            
          } else {
            
            int id = out->id;

            ++out;                                // increment desitination address

            *out = *in;                           // copy first new hit

            out->id = ++id;                       // set hit identifier
          }
        }

        this->erase(++out, this->end());          // remove hits
      }
    }
  };


  /**
   * Auxiliary class to set-up Trk.
   */
  struct JTrk_t : 
    public Trk
  {
    /**
     * Constructor.
     *
     * \param  id         identifier
     * \param  type       type
     * \param  mother_id  mother identifier
     * \param  pos        position
     * \param  dir        direction
     * \param  t          time [ns]
     * \param  E          energy [GeV]
     */
    JTrk_t(const int    id,
           const int    type,
           const int    mother_id,
           const Vec&   pos,
           const Vec&   dir,
           const double t,
           const double E) 
    {
      this->id        = id;
      this->type      = type;
      this->mother_id = mother_id;
      this->pos       = pos;
      this->dir       = dir;
      this->t         = t;
      this->E         = E;
    }
  };
}

#endif