JSireneToolkit.hh

Go to the documentation of this file.

#ifndef __JSIRENE__JSIRENETOOLKIT__
#define __JSIRENE__JSIRENETOOLKIT__

#include <vector>
#include <cmath>

#include "TRandom3.h"

#include "JPhysics/JConstants.hh"
#include "JPhysics/JPDFTypes.hh"
#include "JPhysics/JGeane.hh"
#include "JAAnet/JAAnetToolkit.hh"
#include "JLang/JComparator.hh"


/**
 * \author mdejong
 */


/**
 * Detector simulations.
 */
namespace JSIRENE {}
namespace JPP { using namespace JSIRENE; }

namespace JSIRENE {

  using JPHYSICS::getSpeedOfLight;
  using JPHYSICS::MASS_MUON;
  using JPHYSICS::getIndexOfRefraction;
  using JPHYSICS::JGeane;
  using JPHYSICS::gWater;
  using JPHYSICS::JPDFType_t;
  using JAANET::JHitType_t;


  /**
   * Get number of photo-electrons of a hit given the expectation 
   * values of the number of photo-electrons on a module and PMT.
   *
   * The return value is evaluated by pick-and-drop statistics from
   * the generated number of photo-electrons when the expectation 
   * value of the number of photo-electrons on a module deviates 
   * less than 5 sigmas from 0 (i.e.\ when it is less than 25).
   * Otherwise, the return value is evaluated by Poisson statistics
   * from the expectation value of the number of photo-electrons on PMT.
   *
   * \param  NPE       expectation value of npe on module
   * \param  N         generated   value of npe on module
   * \param  npe       expectation value of npe on PMT
   * \return           number of photo-electrons on PMT
   */
  inline int getNumberOfPhotoElectrons(const double NPE,
                                       const int    N,
                                       const double npe)
  {
    if (NPE < 25.0) {
      
      int n = 0;
                            
      for (int i = N; i != 0; --i) {
        if (gRandom->Rndm() * NPE < npe) {
          ++n;
        }
      }
      
      return n;

    } else {
      
      return gRandom->Poisson(npe);
    }
  }


  /**
   * Get number of photo-electrons of a hit given number of photo-electrons on PMT.
   *
   * The number of photo-electrons of a hit may be larger than unity
   * to limit the overall number of hits and consequently the memory coonsumption as well as 
   * the number of times the arrival time needs to be evaluated which is CPU intensive.
   *
   * \param  npe       number of photo-electrons on PMT
   * \return           number of photo-electrons of hit
   */
  inline int getNumberOfPhotoElectrons(const int npe)
  {
    static const int NPE = 20;
    
    if (npe > NPE) {
      
      const int n = (int) (NPE * log10((double) npe / (double) NPE));
    
      if      (n == 0)
        return 1;
      else
        return n;
    }

    return 1;
  }


  /**
   * Get hit type corresponding to given PDF type.
   *
   * \param  pdf       PDF type
   * \param  shower    force origin from neutrino interaction shower
   * \return           hit type
   */
  inline JHitType_t getHitType(const JPDFType_t pdf, const bool shower = false)
  {
    using namespace JAANET;
    using namespace JPHYSICS;
    
    switch (pdf) {

    case DIRECT_LIGHT_FROM_MUON:
      return HIT_TYPE_MUON_DIRECT;

    case SCATTERED_LIGHT_FROM_MUON:
      return HIT_TYPE_MUON_SCATTERED;

    case DIRECT_LIGHT_FROM_DELTARAYS:
      return HIT_TYPE_DELTARAYS_DIRECT;

    case SCATTERED_LIGHT_FROM_DELTARAYS:
      return HIT_TYPE_DELTARAYS_SCATTERED;

    case DIRECT_LIGHT_FROM_EMSHOWER:
      return (shower ? HIT_TYPE_SHOWER_DIRECT : HIT_TYPE_BREMSSTRAHLUNG_DIRECT );

    case SCATTERED_LIGHT_FROM_EMSHOWER:
      return (shower ? HIT_TYPE_SHOWER_SCATTERED : HIT_TYPE_BREMSSTRAHLUNG_SCATTERED);

    default:
      return HIT_TYPE_UNKNOWN;      
    }    
  }
  

  /**
   * Vertex of energy loss of muon.
   */
  struct JVertex {
    /**
     * Default constructor.
     */
    JVertex() :
      z (0.0),
      t (0.0),
      E (0.0),
      Es(0.0)
    {}


    /**
     * Constructor.
     *
     * \param  z         position [m]
     * \param  t         time     [ns]
     * \param  E         energy   [GeV]
     */
    JVertex(const double z,
            const double t,
            const double E)
    {
      this->z  = z;
      this->t  = t;
      this->E  = E;
      this->Es = 0.0;
    }


    /**
     * Get position.
     *
     * \return           position [m]
     */
    double getZ() const
    {
      return z;
    }


    /**
     * Get time.
     *
     * \return           time [ns]
     */
    double getT() const
    {
      return t;
    }


    /**
     * Get muon energy.
     *
     * \return           energy [GeV]
     */
    double getE() const
    {
      return E;
    }


    /**
     * Get shower energy.
     *
     * \return           energy [GeV]
     */
    double getEs() const
    {
      return Es;
    }


    /**
     * Get range of muon.
     * This method applies only ionisation energy loss.
     *
     * \return           range [m]
     */
    double getRange() const
    {
      if (E > MASS_MUON * getIndexOfRefraction())
        return (E - MASS_MUON * getIndexOfRefraction()) / gWater.getA();
      else
        return 0.0;
    }
    /**
     * Get range of muon using precalculated a parameter
     * This method applies only ionisation energy loss.
     * \param a      ionization a parameter[GeV/m]
     * \return           range [m]
     */
    double getRange(double a) const                                                                                                 {
      if (E > MASS_MUON * getIndexOfRefraction())
        return (E - MASS_MUON * getIndexOfRefraction()) / a;
      else
        return 0.0;                                                                                                         }


    /**
     * Get range of muon.
     *
     * \param  geane     energy loss
     * \return           range [m]
     */
    double getRange(const JGeane& geane) const
    {
      return geane(E);
    }


    /**
     * Apply energy loss energy.
     *
     * \param  Es        energy [GeV]
     */
    void applyEloss(const double Es) 
    {
      this->E  -= Es;
      this->Es  = Es;
    }


    /**
     * Step.
     * This method applies only ionisation energy loss.
     *
     * \param  ds        step [m]
     * \return           this vertex
     */
    JVertex& step(const double ds)
    {
      z += ds;
      t += ds / getSpeedOfLight();
      E -= ds * gWater.getA();

      return *this;
    }
    /* This method applies only ionisation energy loss given a value of a.
     *
     * \param  ds        step [m]
     * \param a        a  ionization parameter [GeV/m]
     * \return           this vertex
     */
    JVertex& step( const double ds, const double a)
    {
      z += ds;
      t += ds / getSpeedOfLight();
      E -= ds * a;

      return *this;
    }
    /**
     * Step.
     *
     * \param  geane     energy loss
     * \param  ds        step [m]
     * \return           this vertex
     */
    JVertex& step(const JGeane& geane, const double ds)
    {
      z += ds;
      t += ds / getSpeedOfLight();
      E  = geane(E, ds);

      return *this;
    }

  protected:
    double z;
    double t;
    double E;    
    double Es;
  };


  /**
   * Muon trajectory.
   */
  struct JTrack :
    public std::vector<JVertex>
  {
    /**
     * Constructor.
     *
     * \param  vertex    muon starting point
     */
    JTrack(const JVertex& vertex) :
      std::vector<JVertex>(1, vertex)
    {}


    /**
     * Get muon energy at given position.
     *
     * \param  z         position [m]
     * \return           energy [GeV]
     */
    double getE(const double z) const
    {
      if (!this->empty()) {
        
        const_iterator pos = lower_bound(this->begin(), this->end(), z, JLANG::make_comparator(&JVertex::getZ));

        if (pos != this->end() && pos != this->begin()) {
          
          --pos;
        
          return pos->getE() -  (z - pos->getZ()) * gWater.getA();
        }
      }

      return MASS_MUON;
    }
  };


  /**
   * Get kinetic energy of particle with given mass.
   *
   * \param  E         energy [GeV]
   * \param  m         mass   [GeV]
   * \return           energy [GeV]
   */
  inline double getKineticEnergy(const double E, const double m)
  {
    if (E > m)
      return sqrt((E - m) * (E + m));
    else
      return 0.0;
  }
}

#endif