JVolume.hh

Go to the documentation of this file.

#ifndef __JGIZMO__JVOLUME__
#define __JGIZMO__JVOLUME__
 
#include <cmath>
 
#include "TAxis.h"
 
#include "km3net-dataformat/offline/Head.hh"
#include "JAAnet/JHead.hh"
 
 
/**
 * \author mdejong
 */
 
namespace JGIZMO {}
namespace JPP { using namespace JGIZMO; }
 
namespace JGIZMO {
 
  /**
   * Auxiliary class for histogramming of effective volume.
   * In this, it is assumed that events are generated according \f$\frac{dN}{dE} \propto E^{\alpha}\f$.
   * The weight is expressed in \f$\mathrm{km}^{3}\f$.
   */
  struct JVolume {
    /**
     * Constructor.
     *
     * \param  head          Monte Carlo run header
     * \param  Elog10        application of log10 to energy
     */
    JVolume(const Head& head, 
            const bool  Elog10 = false) :
      elog (Elog10),
      Emin (1.0),
      Emax (1.0),
      alpha(0.0),
      Wall (1.0)
    {
      using namespace JAANET;
 
      const JHead buffer(head);
 
      if        (is_valid(buffer.spectrum) && 
                 is_valid(buffer.cut_nu)   && 
                 is_valid(buffer.genvol)) {
 
        alpha = buffer.spectrum.alpha;
        Emin  = buffer.cut_nu.Emin;
        Emax  = buffer.cut_nu.Emax;
        Wall  = 1.0e-9 * (getW(Emax) - getW(Emin)) * buffer.genvol.volume / buffer.genvol.numberOfEvents;    // km^3
 
      } else if (is_valid(buffer.cut_in) &&
                 is_valid(buffer.livetime)) {
 
        Emin  = buffer.cut_in.Emin;
        Emax  = buffer.cut_in.Emax;
        Wall  = 1.0 / buffer.livetime.numberOfSeconds;    // Hz
      }
    }
 
 
    /**
     * Get spectral index of energy distribution.
     *
     * \return               spectral index
     */
    inline double getAlpha() const
    {
      return alpha;
    }
 
 
    /**
     * Get generation dependent weight.
     *
     * \return               weight
     */
    inline double getWall() const
    {
      return Wall;
    }
 
 
    /**
     * Get minimal abscissa value.
     *
     * \return               abscissa value
     */
    inline Double_t getXmin() const
    {
      return getX(Emin);
    }
 
 
    /**
     * Get maximal abscissa value.
     *
     * \return               abscissa value
     */
    inline Double_t getXmax() const
    {
      return getX(Emax);
    }
 
 
    /**
     * Get abscissa value.
     *
     * \param  E             energy
     * \param  constrain     constrain
     * \return               abscissa value
     */
    inline Double_t getX(const Double_t E, double constrain = false) const
    {
      const double x = (elog ? log10(E) : E);
 
      if (constrain) {
        if        (x < getXmin()) { 
          return getXmin();
        } else if (x > getXmax()) {
          return getXmax();
        }
      }
 
      return x;
    }
 
 
    /**
     * Get energy.
     *
     * \param  x             abscissa value
     * \param  constrain     constrain
     * \return               energy
     */
    inline Double_t getE(const Double_t x, double constrain = false) const
    {
      const double E = (elog ? pow(10.0, x) : x);
 
      if (constrain) {
        if        (E < Emin) { 
          return Emin;
        } else if (E > Emax) {
          return Emax;
        }
      }
 
      return E;
    }
 
 
    /**
     * Get bin width corrected energy spectrum dependent weight.
     *
     * \param  axis          axis
     * \param  E             energy
     * \return               weight [km^3]
     */
    inline Double_t getW(TAxis* axis, const Double_t E) const
    {
      const Int_t    index = axis->FindBin(getX(E));
 
      const Double_t xmin  = axis->GetBinLowEdge(index);
      const Double_t xmax  = axis->GetBinUpEdge (index);
      
      const Double_t Wmin  = getW(getE(xmin));
      const Double_t Wmax  = getW(getE(xmax));
 
      return Wall / (Wmax - Wmin);
    }
 
 
    /**
     * Get generation dependent integral value of given energy.
     *
     * \param  E             energy
     * \return               integral value
     */
    inline double getW(const double E) const
    {
      if (alpha != -1.0)
        return pow(E, 1.0 + alpha) / (1.0 + alpha);
      else
        return log(E);
    }
 
 
    /**
     * Multiply weight.
     *
     * \param  factor        factor
     * \return               this volume
     */
    JVolume& mul(const double factor) 
    {
      Wall *= factor;
 
      return *this;
    }
 
 
    /**
     * Divide weight.
     *
     * \param  factor        factor
     * \return               this volume
     */
    JVolume& div(const double factor) 
    {
      Wall /= factor;
 
      return *this;
    }
 
 
  protected:
    bool   elog;        //!< histogram option
    double Emin;        //!< minimal energy
    double Emax;        //!< maximal energy
    double alpha;       //!< spectral index
    double Wall;        //!< generation volume
  };
}
 
#endif