EgammaCoreTools/src/SuperClusterShapeAlgo.cc

0001 #include <iostream>
0002
0003 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/SuperClusterShapeAlgo.h"
0005 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0006 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0007 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0008 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0009
0010 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0011 #include "DataFormats/Math/interface/Point3D.h"
0012 #include "DataFormats/Math/interface/Vector3D.h"
0013
0014 SuperClusterShapeAlgo::SuperClusterShapeAlgo(const EcalRecHitCollection* hits, const CaloSubdetectorGeometry* geo)
0015     : recHits_(hits), geometry_(geo) {}
0016
0017 void SuperClusterShapeAlgo::Calculate_Covariances(const reco::SuperCluster& passedCluster) {
0018   double numeratorEtaWidth = 0;
0019   double numeratorPhiWidth = 0;
0020
0021   double scEnergy = passedCluster.energy();
0022   double denominator = scEnergy;
0023
0024   double scEta = passedCluster.position().eta();
0025   double scPhi = passedCluster.position().phi();
0026
0027   const std::vector<std::pair<DetId, float> >& detId = passedCluster.hitsAndFractions();
0028   // Loop over recHits associated with the given SuperCluster
0029   for (std::vector<std::pair<DetId, float> >::const_iterator hit = detId.begin(); hit != detId.end(); ++hit) {
0030     EcalRecHitCollection::const_iterator rHit = recHits_->find((*hit).first);
0031     //FIXME: THIS IS JUST A WORKAROUND A FIX SHOULD BE APPLIED
0032     if (rHit == recHits_->end()) {
0033       continue;
0034     }
0035     auto this_cell = geometry_->getGeometry(rHit->id());
0036     if (this_cell == nullptr) {
0037       //edm::LogInfo("SuperClusterShapeAlgo") << "pointer to the cell in Calculate_Covariances is NULL!";
0038       continue;
0039     }
0040     GlobalPoint position = this_cell->getPosition();
0041     //take into account energy fractions
0042     double energyHit = rHit->energy() * hit->second;
0043
0044     //form differences
0045     double dPhi = position.phi() - scPhi;
0046     if (dPhi > +Geom::pi()) {
0047       dPhi = Geom::twoPi() - dPhi;
0048     }
0049     if (dPhi < -Geom::pi()) {
0050       dPhi = Geom::twoPi() + dPhi;
0051     }
0052
0053     double dEta = position.eta() - scEta;
0054
0055     if (energyHit > 0) {
0056       numeratorEtaWidth += energyHit * dEta * dEta;
0057       numeratorPhiWidth += energyHit * dPhi * dPhi;
0058     }
0059
0060     etaWidth_ = sqrt(numeratorEtaWidth / denominator);
0061     phiWidth_ = sqrt(numeratorPhiWidth / denominator);
0062   }
0063 }