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        (buffer.is_valid(&JHead::spectrum) && 
                 buffer.is_valid(&JHead::cut_nu)   && 
                 buffer.is_valid(&JHead::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 (buffer.is_valid(&JHead::cut_in) &&
                 buffer.is_valid(&JHead::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