TOPPmtTTSPar.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <TObject.h>
#include <framework/logging/Logger.h>
#include <vector>
 
#include <TRandom.h>
 
namespace Belle2 {
  class TOPPmtTTSPar : public TObject {
 
  public:
    struct Gaussian {
      float fraction = 0; 
      float mean = 0;     
      float sigma = 0;    
    };
 
    enum {c_NumPmtPixels = 16};
 
    TOPPmtTTSPar()
    {}
 
    explicit TOPPmtTTSPar(const std::string& serialNumber):
      m_serialNumber(serialNumber)
    {}
 
 
    void appendGaussian(unsigned pmtPixel, const Gaussian& gaus)
    {
      pmtPixel--;
      if (pmtPixel >= c_NumPmtPixels) {
        B2ERROR("TOPPmtTTSPar::appendGaussian: invalid PMT pixel "
                << LogVar("PMT pixel", pmtPixel + 1));
        return;
      }
      m_gaussians[pmtPixel].push_back(gaus);
    }
 
 
    void appendGaussian(unsigned pmtPixel, double fraction, double mean, double sigma)
    {
      Gaussian gaus;
      gaus.fraction = fraction;
      gaus.mean = mean;
      gaus.sigma = sigma;
      appendGaussian(pmtPixel, gaus);
    }
 
    int getNumOfPixels() const {return c_NumPmtPixels;}
 
    const std::string& getSerialNumber() const {return m_serialNumber;}
 
 
    const std::vector<Gaussian>& getGaussians(unsigned pmtPixel) const
    {
      pmtPixel--;
      if (pmtPixel >= c_NumPmtPixels) {
        B2ERROR("TOPPmtTTSPar::getGaussians: invalid PMT pixel. "
                "Returning data for pixel 1."
                << LogVar("PMT pixel", pmtPixel + 1));
        pmtPixel = 0;
      }
      return m_gaussians[pmtPixel];
    }
 
 
    double getRandomTime(unsigned pmtPixel) const
    {
      pmtPixel--;
      if (pmtPixel >= c_NumPmtPixels) return 0;
 
      double prob = gRandom->Rndm();
      double s = 0;
      for (const auto& gaus : m_gaussians[pmtPixel]) {
        s = s + gaus.fraction;
        if (prob < s) {
          return gRandom->Gaus(gaus.mean, gaus.sigma);
        }
      }
      return 0;  // this should not happen, if fractions are properly normalized
    }
 
 
    void normalizeFractions()
    {
      double sum = 0;
      for (int ich = 0 ; ich < c_NumPmtPixels ; ich++) {
        for (const auto& gaus : m_gaussians[ich]) {
          sum = sum + gaus.fraction;
        }
        if (sum == 0) return;
        for (auto& gaus : m_gaussians[ich]) {
          gaus.fraction = gaus.fraction / sum;
        }
      }
      return;
    }
 
 
 
  private:
 
    std::string m_serialNumber;                          
    std::vector<Gaussian> m_gaussians[c_NumPmtPixels];    
    ClassDef(TOPPmtTTSPar, 2); 
  };
 
} // end namespace Belle2