SdSimulation/SdAccidentalInjectorKG/Utilities.h

Go to the documentation of this file.

#ifndef _AccidentialInjectorKG_Utilities_h_
#define _AccidentialInjectorKG_Utilities_h_

#include <utl/Math.h>
#include <utl/MathConstants.h>
#include <utl/RandomEngine.h>
#include <fwk/RandomEngineRegistry.h>

#include <CLHEP/Random/Randomize.h>

#include <TH2D.h>

#include <vector>
#include <algorithm>
#include <cmath>


namespace SdAccidentalInjectorKG {

  inline
  std::pair<double, double>
  StationTopSideArea(const double radius, const double height, const double cosTheta)
  {
    const double sinTheta = std::sqrt(1 - utl::Sqr(cosTheta));
    const double areaTop  = cosTheta * utl::kPi * utl::Sqr(radius);
    const double areaSide = 2 * sinTheta * radius * height;
    return std::make_pair(areaTop, areaSide);
  }


  // station area projected into the plane perpendicular to the particle
  inline
  double
  StationArea(const double radius, const double height, const double cosTheta)
  {
    const auto ts = StationTopSideArea(radius, height, cosTheta);
    return ts.first + ts.second;
  }


  class FlatSampler {
  public:
    FlatSampler() :
      fEngine(&fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector).GetEngine())
    { }

    double operator()(const double vmin, const double vmax) const
    { return CLHEP::RandFlat::shoot(fEngine, vmin, vmax); }

    double operator()(const double vmax) const
    { return CLHEP::RandFlat::shoot(fEngine, 0, vmax); }

    double operator()() const
    { return CLHEP::RandFlat::shoot(fEngine, 0, 1); }

  protected:
    utl::RandomEngine::RandomEngineType* fEngine = nullptr;
  };


  class PoissonSampler {
  public:
    PoissonSampler() :
      fEngine(&fwk::RandomEngineRegistry::GetInstance().Get(fwk::RandomEngineRegistry::eDetector).GetEngine())
    { }

    double operator()(const double expectation) const
    { return CLHEP::RandPoisson::shoot(fEngine, expectation); }

  protected:
    utl::RandomEngine::RandomEngineType* fEngine = nullptr;
  };


  /* Strategy:
   * The PDF is available as a 2D histogram. We should use the Unuran
   * library to sample from it, but the ROOT interface TUnuran supports
   * only 1D histograms at the time of this writing.
   * We therefore apply a custom strategy.
   *
   * We integrate over the 2D histogram with a line integral. The line
   * integral is then inverted as a function of its upper bound.
   * The inverted function then allows us to convert a uniform random
   * number to the desired random 2D sample.
   *
   * Note that this approach can be generalized to higher dimensions,
   * but the numerical quality is expected to become worse as the
   * number of dimensions increases.
   */

  class HistogramSampler {
  public:
    HistogramSampler() { }

    HistogramSampler(const TH2D& h) :
      fHistogram(h)
    {
      const unsigned int nx = h.GetNbinsX();
      const unsigned int ny = h.GetNbinsY();
      fLineIntegral.reserve(nx*ny);
      double sum = 0;
      for (unsigned int ix = 1; ix <= nx; ++ix)
        for (unsigned int iy = 1; iy <= ny; ++iy) {
          sum += h.GetBinContent(ix, iy);
          fLineIntegral.push_back(sum);
        }
    }

    std::pair<double, double>
    operator()()
      const
    {
      const double z = fFlat(fLineIntegral.back());
      // this "inverts" the line integral
      const auto begin = fLineIntegral.begin();
      const int i = std::upper_bound(begin, fLineIntegral.end(), z) - begin;
      const int n = fHistogram.GetNbinsY();
      const int ix = (i / n) + 1;
      const int iy = (i % n) + 1;

      // do flat sampling inside bin
      const auto& ax = *fHistogram.GetXaxis();
      const double x = fFlat(ax.GetBinLowEdge(ix), ax.GetBinUpEdge(ix));
      const auto& ay = *fHistogram.GetYaxis();
      const double y = fFlat(ay.GetBinLowEdge(iy), ay.GetBinUpEdge(iy));
      return std::make_pair(x, y);
    }

  protected:
    TH2D fHistogram;
    std::vector<double> fLineIntegral;
    FlatSampler fFlat;
  };


  class CylinderSurfaceSampler {
  public:
    CylinderSurfaceSampler() { }

    CylinderSurfaceSampler(const double radius, const double height, const double hull) :
      // adding walls
      fRadius(radius + hull),
      fRadius2(utl::Sqr(fRadius)),
      fHeight(height + 2*hull),
      fHull(hull)
    { }

    std::tuple<double, double, double>
    operator()(const double theta, const double phi)
      const
    {
      const auto area = StationTopSideArea(fRadius, fHeight, std::cos(theta));
      const double totArea = area.first + area.second;
      // top or side?
      if (fFlat() < area.first / totArea) {
        // top
        const double r = fRadius * std::sqrt(fFlat());
        const double f = fFlat(2 * utl::kPi);
        return std::make_tuple(r * std::sin(f), r * std::cos(f), fHeight);
      }
      // side
      const double y = fFlat(-fRadius, fRadius);
      const double x = std::sqrt(fRadius2 - utl::Sqr(y));
      const double z = fFlat(fHeight);
      const double c = std::cos(phi);
      const double s = std::sin(phi);
      return std::make_tuple(x*c - y*s, x*s + y*c, z);
    }

  protected:
    double fRadius = 0;
    double fRadius2 = 0;
    double fHeight = 0;
    double fHull = 0;
    FlatSampler fFlat;
  };

}


#endif