EgammaCoreTools/interface/PositionCalc.h

0001 #ifndef RecoEcal_EgammaCoreTools_PositionCalc_h
0002 #define RecoEcal_EgammaCoreTools_PositionCalc_h
0003
0004 /** \class PositionClac
0005  *
0006  * Finds the position and covariances for a cluster
0007  * Formerly LogPositionCalc
0008  *
0009  * \author Ted Kolberg, ND
0010  *
0011  *
0012  */
0013
0014 #include <vector>
0015 #include <map>
0016
0017 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0018 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0019 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0020 #include "DataFormats/Math/interface/Point3D.h"
0021 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0024 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0027 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0028
0029 class PositionCalc {
0030 public:
0031   typedef std::vector<std::pair<DetId, float> > HitsAndFractions;
0032   typedef std::vector<std::pair<DetId, double> > HitsAndEnergies;
0033   // You must call Initialize before you can calculate positions or
0034   // covariances.
0035
0036   PositionCalc(const edm::ParameterSet& par);
0037   PositionCalc() {}
0038
0039   const PositionCalc& operator=(const PositionCalc& rhs);
0040
0041   // Calculate_Location calculates an arithmetically or logarithmically
0042   // weighted average position of a vector of DetIds, which should be
0043   // a subset of the map used to Initialize.
0044
0045   template <typename HitType>
0046   math::XYZPoint Calculate_Location(const HitsAndFractions& iDetIds,
0047                                     const edm::SortedCollection<HitType>* iRecHits,
0048                                     const CaloSubdetectorGeometry* iSubGeom,
0049                                     const CaloSubdetectorGeometry* iESGeom = nullptr);
0050
0051 private:
0052   bool param_LogWeighted_;
0053   double param_T0_barl_;
0054   double param_T0_endc_;
0055   double param_T0_endcPresh_;
0056   double param_W0_;
0057   double param_X0_;
0058
0059   const CaloSubdetectorGeometry* m_esGeom;
0060   bool m_esPlus;
0061   bool m_esMinus;
0062 };
0063
0064 template <typename HitType>
0065 math::XYZPoint PositionCalc::Calculate_Location(const PositionCalc::HitsAndFractions& iDetIds,
0066                                                 const edm::SortedCollection<HitType>* iRecHits,
0067                                                 const CaloSubdetectorGeometry* iSubGeom,
0068                                                 const CaloSubdetectorGeometry* iESGeom) {
0069   typedef edm::SortedCollection<HitType> HitTypeCollection;
0070   math::XYZPoint returnValue(0, 0, 0);
0071
0072   // Throw an error if the cluster was not initialized properly
0073
0074   if (nullptr == iRecHits || nullptr == iSubGeom) {
0075     throw cms::Exception("PositionCalc") << "Calculate_Location() called uninitialized or wrong initialization.";
0076   }
0077
0078   if (!iDetIds.empty() && !iRecHits->empty()) {
0079     HitsAndEnergies detIds;
0080     detIds.reserve(iDetIds.size());
0081
0082     double eTot(0);
0083     double eMax(0);
0084     DetId maxId;
0085
0086     // Check that DetIds are nonzero
0087     typename HitTypeCollection::const_iterator endRecHits(iRecHits->end());
0088     HitsAndFractions::const_iterator n, endDiDs(iDetIds.end());
0089     for (n = iDetIds.begin(); n != endDiDs; ++n) {
0090       const DetId dId((*n).first);
0091       const float frac((*n).second);
0092       if (!dId.null()) {
0093         typename HitTypeCollection::const_iterator iHit(iRecHits->find(dId));
0094         if (iHit != endRecHits) {
0095           const double energy(iHit->energy() * frac);
0096           detIds.push_back(std::make_pair(dId, energy));
0097           if (0.0 < energy) {  // only save positive energies
0098             if (eMax < energy) {
0099               eMax = energy;
0100               maxId = dId;
0101             }
0102             eTot += energy;
0103           }
0104         }
0105       }
0106     }
0107
0108     if (0.0 >= eTot) {
0109       LogDebug("ZeroClusterEnergy") << "cluster with 0 energy: " << eTot << " size: " << detIds.size()
0110                                     << " , returning (0,0,0)";
0111     } else {
0112       // first time or when es geom changes set flags
0113       if (nullptr != iESGeom && m_esGeom != iESGeom) {
0114         m_esGeom = iESGeom;
0115         for (uint32_t ic(0); (ic != m_esGeom->getValidDetIds().size()) && ((!m_esPlus) || (!m_esMinus)); ++ic) {
0116           const double z(m_esGeom->getGeometry(m_esGeom->getValidDetIds()[ic])->getPosition().z());
0117           m_esPlus = m_esPlus || (0 < z);
0118           m_esMinus = m_esMinus || (0 > z);
0119         }
0120       }
0121
0122       //Select the correct value of the T0 parameter depending on subdetector
0123       auto center_cell(iSubGeom->getGeometry(maxId));
0124       const double ctreta(center_cell->getPosition().eta());
0125
0126       // for barrel, use barrel T0;
0127       // for endcap: if preshower present && in preshower fiducial,
0128       //             use preshower T0
0129       // else use endcap only T0
0130       const double preshowerStartEta = 1.653;
0131       const int subdet = maxId.subdetId();
0132       const double T0(subdet == EcalBarrel ? param_T0_barl_
0133                                            : (((preshowerStartEta < fabs(ctreta)) &&
0134                                                (((0 < ctreta) && m_esPlus) || ((0 > ctreta) && m_esMinus)))
0135                                                   ? param_T0_endcPresh_
0136                                                   : param_T0_endc_));
0137
0138       // Calculate shower depth
0139       const float maxDepth(param_X0_ * (T0 + log(eTot)));
0140       const float maxToFront(center_cell->getPosition().mag());  // to front face
0141
0142       // Loop over hits and get weights
0143       double total_weight = 0;
0144       const double eTot_inv = 1.0 / eTot;
0145       const double logETot_inv = (param_LogWeighted_ ? log(eTot_inv) : 0);
0146
0147       double xw(0);
0148       double yw(0);
0149       double zw(0);
0150
0151       HitsAndEnergies::const_iterator j, hAndE_end = detIds.end();
0152       for (j = detIds.begin(); j != hAndE_end; ++j) {
0153         const DetId dId((*j).first);
0154         const double e_j((*j).second);
0155
0156         double weight = 0;
0157         if (param_LogWeighted_) {
0158           if (e_j > 0.0) {
0159             weight = std::max(0., param_W0_ + log(e_j) + logETot_inv);
0160           } else {
0161             weight = 0;
0162           }
0163         } else {
0164           weight = e_j * eTot_inv;
0165         }
0166
0167         auto cell(iSubGeom->getGeometry(dId));
0168         const float depth(maxDepth + maxToFront - cell->getPosition().mag());
0169
0170         const GlobalPoint pos(cell->getPosition(depth));
0171
0172         xw += weight * pos.x();
0173         yw += weight * pos.y();
0174         zw += weight * pos.z();
0175
0176         total_weight += weight;
0177       }
0178       returnValue = math::XYZPoint(xw / total_weight, yw / total_weight, zw / total_weight);
0179     }
0180   }
0181   return returnValue;
0182 }
0183
0184 #endif