Smearers.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file Smearers.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include "Acts/Utilities/BinningData.hpp"
#include "Acts/Utilities/Result.hpp"
#include "ActsExamples/Framework/RandomNumbers.hpp"
#include "ActsFatras/Digitization/DigitizationError.hpp"
#include "ActsFatras/Digitization/UncorrelatedHitSmearer.hpp"

#include <cmath>
#include <limits>
#include <random>
#include <string>
#include <tuple>
#include <utility>

namespace ActsExamples {
namespace Digitization {

struct Exact {
  Acts::Result<std::pair<double, double>> operator()(
      double value, RandomEngine& /*rnd*/) const {
    return std::pair{value, 0.0};
  }
};

struct Gauss {
  double sigma;

  Gauss(double sigma_) : sigma{sigma_} {}

  Acts::Result<std::pair<double, double>> operator()(double value,
                                                     RandomEngine& rnd) const {
    std::normal_distribution<> dist{0, sigma};
    return std::pair{value + dist(rnd), dist.stddev()};
  }
};

struct GaussTrunc {
  double sigma;
  std::pair<double, double> range = {std::numeric_limits<double>::lowest(),
                                     std::numeric_limits<double>::max()};

  GaussTrunc(double sigma_, const std::pair<double, double>& range_)
      : sigma{sigma_}, range(range_) {}

  Acts::Result<std::pair<double, double>> operator()(double value,
                                                     RandomEngine& rnd) const {
    std::normal_distribution<> dist{0., sigma};
    double svalue = value + dist(rnd);
    if (svalue >= range.first and svalue <= range.second) {
      return std::pair{svalue, dist.stddev()};
    }
    return ActsFatras::DigitizationError::SmearingOutOfRange;
  }
};

struct GaussClipped {
  double sigma;

  size_t maxAttemps = 1000;

  std::pair<double, double> range = {std::numeric_limits<double>::lowest(),
                                     std::numeric_limits<double>::max()};

  GaussClipped(double sigma_, const std::pair<double, double>& range_)
      : sigma{sigma_}, range(range_) {}

  Acts::Result<std::pair<double, double>> operator()(double value,
                                                     RandomEngine& rnd) const {
    std::normal_distribution<> dist{0., sigma};
    for (size_t attempt = 0; attempt < maxAttemps; ++attempt) {
      double svalue = value + dist(rnd);
      if (svalue >= range.first and svalue <= range.second) {
        return std::pair{svalue, dist.stddev()};
      }
    }
    return ActsFatras::DigitizationError::SmearingError;
  }
};

struct Uniform {
  Acts::BinningData binningData;

  Uniform(double pitch, const std::pair<double, double>& range_)
      : binningData(Acts::open, Acts::binX,
                    static_cast<size_t>((range_.second - range_.first) / pitch),
                    range_.first, range_.second) {}

  Uniform(Acts::BinningData&& bd) : binningData(bd) {}

  Acts::Result<std::pair<double, double>> operator()(double value,
                                                     RandomEngine& rnd) const {
    if (binningData.min < value and binningData.max > value) {
      auto bin = binningData.search(value);
      auto lower = binningData.boundaries()[bin];
      auto higher = binningData.boundaries()[bin + 1];
      std::uniform_real_distribution<> dist{0., 1.};
      double svalue = lower + (higher - lower) * dist(rnd);
      return std::pair{svalue, (higher - lower) / std::sqrt(12.)};
    }
    return ActsFatras::DigitizationError::SmearingError;
  }
};

struct Digital {
  Acts::BinningData binningData;

  Digital(double pitch, const std::pair<double, double>& range_)
      : binningData(Acts::open, Acts::binX,
                    static_cast<size_t>((range_.second - range_.first) / pitch),
                    range_.first, range_.second) {}

  Digital(Acts::BinningData&& bd) : binningData(bd) {}

  Acts::Result<std::pair<double, double>> operator()(
      double value, RandomEngine& /*rnd*/) const {
    if (binningData.min < value and binningData.max > value) {
      auto bin = binningData.search(value);
      auto lower = binningData.boundaries()[bin];
      auto higher = binningData.boundaries()[bin + 1];
      double svalue = 0.5 * (lower + higher);
      return std::pair{svalue, (higher - lower) / std::sqrt(12.)};
    }
    return ActsFatras::DigitizationError::SmearingError;
  }
};

}  // namespace Digitization
}  // namespace ActsExamples