PositionCalc.h

Go to the documentation of this file.

#ifndef RecoEcal_EgammaCoreTools_PositionCalc_h
#define RecoEcal_EgammaCoreTools_PositionCalc_h

#include <vector>
#include <map>

#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/Math/interface/Point3D.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

class PositionCalc {
public:
  typedef std::vector<std::pair<DetId, float> > HitsAndFractions;
  typedef std::vector<std::pair<DetId, double> > HitsAndEnergies;
  // You must call Initialize before you can calculate positions or
  // covariances.

  PositionCalc(const edm::ParameterSet& par);
  PositionCalc() {}

  const PositionCalc& operator=(const PositionCalc& rhs);

  // Calculate_Location calculates an arithmetically or logarithmically
  // weighted average position of a vector of DetIds, which should be
  // a subset of the map used to Initialize.

  template <typename HitType>
  math::XYZPoint Calculate_Location(const HitsAndFractions& iDetIds,
                                    const edm::SortedCollection<HitType>* iRecHits,
                                    const CaloSubdetectorGeometry* iSubGeom,
                                    const CaloSubdetectorGeometry* iESGeom = nullptr);

private:
  bool param_LogWeighted_;
  double param_T0_barl_;
  double param_T0_endc_;
  double param_T0_endcPresh_;
  double param_W0_;
  double param_X0_;

  const CaloSubdetectorGeometry* m_esGeom;
  bool m_esPlus;
  bool m_esMinus;
};

template <typename HitType>
math::XYZPoint PositionCalc::Calculate_Location(const PositionCalc::HitsAndFractions& iDetIds,
                                                const edm::SortedCollection<HitType>* iRecHits,
                                                const CaloSubdetectorGeometry* iSubGeom,
                                                const CaloSubdetectorGeometry* iESGeom) {
  typedef edm::SortedCollection<HitType> HitTypeCollection;
  math::XYZPoint returnValue(0, 0, 0);

  // Throw an error if the cluster was not initialized properly

  if (nullptr == iRecHits || nullptr == iSubGeom) {
    throw cms::Exception("PositionCalc") << "Calculate_Location() called uninitialized or wrong initialization.";
  }

  if (!iDetIds.empty() && !iRecHits->empty()) {
    HitsAndEnergies detIds;
    detIds.reserve(iDetIds.size());

    double eTot(0);
    double eMax(0);
    DetId maxId;

    // Check that DetIds are nonzero
    typename HitTypeCollection::const_iterator endRecHits(iRecHits->end());
    HitsAndFractions::const_iterator n, endDiDs(iDetIds.end());
    for (n = iDetIds.begin(); n != endDiDs; ++n) {
      const DetId dId((*n).first);
      const float frac((*n).second);
      if (!dId.null()) {
        typename HitTypeCollection::const_iterator iHit(iRecHits->find(dId));
        if (iHit != endRecHits) {
          const double energy(iHit->energy() * frac);
          detIds.push_back(std::make_pair(dId, energy));
          if (0.0 < energy) {  // only save positive energies
            if (eMax < energy) {
              eMax = energy;
              maxId = dId;
            }
            eTot += energy;
          }
        }
      }
    }

    if (0.0 >= eTot) {
      LogDebug("ZeroClusterEnergy") << "cluster with 0 energy: " << eTot << " size: " << detIds.size()
                                    << " , returning (0,0,0)";
    } else {
      // first time or when es geom changes set flags
      if (nullptr != iESGeom && m_esGeom != iESGeom) {
        m_esGeom = iESGeom;
        for (uint32_t ic(0); (ic != m_esGeom->getValidDetIds().size()) && ((!m_esPlus) || (!m_esMinus)); ++ic) {
          const double z(m_esGeom->getGeometry(m_esGeom->getValidDetIds()[ic])->getPosition().z());
          m_esPlus = m_esPlus || (0 < z);
          m_esMinus = m_esMinus || (0 > z);
        }
      }

      //Select the correct value of the T0 parameter depending on subdetector
      auto center_cell(iSubGeom->getGeometry(maxId));
      const double ctreta(center_cell->getPosition().eta());

      // for barrel, use barrel T0;
      // for endcap: if preshower present && in preshower fiducial,
      //             use preshower T0
      // else use endcap only T0
      const double preshowerStartEta = 1.653;
      const int subdet = maxId.subdetId();
      const double T0(subdet == EcalBarrel ? param_T0_barl_
                                           : (((preshowerStartEta < fabs(ctreta)) &&
                                               (((0 < ctreta) && m_esPlus) || ((0 > ctreta) && m_esMinus)))
                                                  ? param_T0_endcPresh_
                                                  : param_T0_endc_));

      // Calculate shower depth
      const float maxDepth(param_X0_ * (T0 + log(eTot)));
      const float maxToFront(center_cell->getPosition().mag());  // to front face

      // Loop over hits and get weights
      double total_weight = 0;
      const double eTot_inv = 1.0 / eTot;
      const double logETot_inv = (param_LogWeighted_ ? log(eTot_inv) : 0);

      double xw(0);
      double yw(0);
      double zw(0);

      HitsAndEnergies::const_iterator j, hAndE_end = detIds.end();
      for (j = detIds.begin(); j != hAndE_end; ++j) {
        const DetId dId((*j).first);
        const double e_j((*j).second);

        double weight = 0;
        if (param_LogWeighted_) {
          if (e_j > 0.0) {
            weight = std::max(0., param_W0_ + log(e_j) + logETot_inv);
          } else {
            weight = 0;
          }
        } else {
          weight = e_j * eTot_inv;
        }

        auto cell(iSubGeom->getGeometry(dId));
        const float depth(maxDepth + maxToFront - cell->getPosition().mag());

        const GlobalPoint pos(cell->getPosition(depth));

        xw += weight * pos.x();
        yw += weight * pos.y();
        zw += weight * pos.z();

        total_weight += weight;
      }
      returnValue = math::XYZPoint(xw / total_weight, yw / total_weight, zw / total_weight);
    }
  }
  return returnValue;
}

#endif