CDCDedx2DCell.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 <framework/logging/Logger.h>
 
#include <TObject.h>
#include <TH2F.h>
 
namespace Belle2 {
  class CDCDedx2DCell: public TObject {
 
  public:
 
    CDCDedx2DCell(): m_version(), m_twodgains() {};
 
    CDCDedx2DCell(short version, const std::vector<TH2F>& twodgains): m_version(version), m_twodgains(twodgains) {};
 
    ~CDCDedx2DCell() {};
 
    CDCDedx2DCell& operator*=(CDCDedx2DCell const& rhs)
    {
      if (m_version != rhs.getVersion()) {
        B2WARNING("2D cell gain parameters do not match, cannot merge!");
        return *this;
      }
      for (unsigned int layer = 0; layer < getSize(); ++layer) {
        const TH2F* newhist = rhs.getHist(layer);
        if (newhist->GetEntries() > 0)m_twodgains[layer].Multiply(newhist);
        else B2ERROR("ERROR! constant histograms is empty");
      }
      return *this;
    }
 
    short getVersion() const {return m_version; };
 
    unsigned int getSize() const { return m_twodgains.size(); };
 
    const TH2F* getHist(unsigned int layer) const
    {
      if (m_twodgains.size() == 0) {
        B2ERROR("CDCDedx2DCell: vector too short");
        return nullptr;
      }
 
      if (m_twodgains.size() == 2) {
        if (layer == 0) return &m_twodgains[0];
        else if (layer == 1) return &m_twodgains[1];
        else {
          B2ERROR("CDCDedx2DCell: const histograms NOT found");
          return nullptr;
        }
      } else if (m_twodgains.size() == 56) {
        if (layer < m_twodgains.size()) return &m_twodgains[layer];
      } else {
        B2ERROR("CDCDedx2DCell: Something is wrong # of const histograms");
        return nullptr;
      }
      B2ERROR("CDCDedx2DCell: wrong layer number");
      return nullptr;
    }
 
    double getMean(unsigned int layer, int dbin, int ebin) const
    {
      unsigned int mylayer = 0;
      if (layer >= 8 && m_twodgains.size() == 2) mylayer = 1;
      else if (m_twodgains.size() == 56) mylayer = layer;
 
      if (mylayer < m_twodgains.size()) return m_twodgains[mylayer].GetBinContent(dbin, ebin);
 
      B2ERROR("CDCDedx2DCell: wrong layer number, returning 1");
      return 1.0;
    }
 
    double getMean(unsigned int layer, double doca, double enta) const
    {
      if (layer > 56) {
        B2ERROR("Asking for a CDC layer beyond 56!");
        return 0;
      }
 
      unsigned int mylayer = 0;
      if (layer >= 8 && m_twodgains.size() == 2) mylayer = 1;
      else if (m_twodgains.size() == 56) mylayer = layer;
 
      if (mylayer >= m_twodgains.size()) {
        B2ERROR("CDCDedx2DCell: wrong layer number, returning 1");
        return 1.0;
      }
 
      // assume rotational symmetry
      if (enta < -3.1416 / 2.0) enta += 3.1416 / 2.0;
      if (enta > 3.1416 / 2.0) enta -= 3.1416 / 2.0;
 
      // iterate bin number by one since TH2F bins start at 1 not 0
      double dbinsize = m_twodgains[mylayer].GetXaxis()->GetBinCenter(2) - m_twodgains[mylayer].GetXaxis()->GetBinCenter(1);
      int dbin = std::floor((doca - m_twodgains[mylayer].GetXaxis()->GetBinLowEdge(1)) / dbinsize) + 1;
 
      double ebinsize = m_twodgains[mylayer].GetYaxis()->GetBinCenter(2) - m_twodgains[mylayer].GetYaxis()->GetBinCenter(1);
      // int ebin = std::floor((std::sin(enta) - m_twodgains[mylayer].GetYaxis()->GetBinLowEdge(1)) / ebinsize) + 1;
      int ebin = std::floor((enta - m_twodgains[mylayer].GetYaxis()->GetBinLowEdge(1)) / ebinsize) + 1;
 
      double mean = 1.0;
      if (dbin > 0 && dbin <= m_twodgains[mylayer].GetNbinsX() && ebin > 0 && ebin <= m_twodgains[mylayer].GetNbinsY())
        mean = m_twodgains[mylayer].GetBinContent(dbin, ebin);
      else if (dbin > m_twodgains[mylayer].GetNbinsX() && ebin > 0 && ebin <= m_twodgains[mylayer].GetNbinsY())
        mean = m_twodgains[mylayer].GetBinContent(m_twodgains[mylayer].GetNbinsX(), ebin);
      else
        // B2WARNING("Problem with 2D CDC dE/dx calibration! " << doca << "\t" << dbin << "\t" << std::sin(enta) << "\t" << ebin);
        B2WARNING("Problem with 2D CDC dE/dx calibration! " << doca << "\t" << dbin << "\t" << enta << "\t" << ebin);
 
      return mean;
    };
 
  private:
    short m_version; 
    std::vector<TH2F> m_twodgains; 
    ClassDef(CDCDedx2DCell, 5); 
  };
} // end namespace Belle2