EgammaClusterAlgos/src/PreshowerPhiClusterAlgo.cc

0001 #include <vector>
0002 #include <map>
0003 #include <iostream>
0004
0005 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0006 #include "RecoEcal/EgammaClusterAlgos/interface/PreshowerPhiClusterAlgo.h"
0007 #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
0008 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0009 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "DataFormats/Math/interface/deltaPhi.h"
0012
0013 #include "TH1.h"
0014
0015 reco::PreshowerCluster PreshowerPhiClusterAlgo::makeOneCluster(ESDetId strip,
0016                                                                HitsID* used_strips,
0017                                                                RecHitsMap* the_rechitsMap_p,
0018                                                                const CaloSubdetectorGeometry* geometry_p,
0019                                                                double deltaEta,
0020                                                                double minDeltaPhi,
0021                                                                double maxDeltaPhi) {
0022   rechits_map = the_rechitsMap_p;
0023   used_s = used_strips;
0024
0025   int plane = strip.plane();
0026
0027   // create null-cluster
0028   std::vector<std::pair<DetId, float> > dummy;
0029   Point posi(0, 0, 0);
0030   reco::PreshowerCluster nullcluster = reco::PreshowerCluster(0., posi, dummy, plane);
0031   if (strip == ESDetId(0))
0032     return nullcluster;  // works in case of no intersected strip found (e.g. in the Barrel)
0033
0034   auto refCell = geometry_p->getGeometry(strip);
0035   const GlobalPoint& refpos = refCell->getPosition();
0036   double refEta = refpos.eta();
0037   double refPhi = refpos.phi();
0038
0039   // Collection of cluster strips
0040   EcalRecHitCollection clusterRecHits;
0041
0042   double x_pos = 0;
0043   double y_pos = 0;
0044   double z_pos = 0;
0045
0046   RecHitsMap::iterator strip_it;
0047   for (strip_it = rechits_map->begin(); strip_it != rechits_map->end(); ++strip_it) {
0048     if (!goodStrip(strip_it))
0049       continue;
0050
0051     ESDetId mystrip = (strip_it->first == DetId(0)) ? ESDetId(0) : ESDetId(strip_it->first);
0052     if (mystrip.plane() != plane)
0053       continue;
0054
0055     auto thisCell = geometry_p->getGeometry(strip_it->first);
0056     const GlobalPoint& position = thisCell->getPosition();
0057
0058     if (fabs(position.eta() - refEta) < deltaEta) {
0059       //std::cout<<"all strips : "<<mystrip.plane()<<" "<<position.phi()<<" "<<reco::deltaPhi(position.phi(), refPhi)<<std::endl;
0060
0061       if (reco::deltaPhi(position.phi(), refPhi) > 0 && reco::deltaPhi(position.phi(), refPhi) > maxDeltaPhi)
0062         continue;
0063       if (reco::deltaPhi(position.phi(), refPhi) < 0 && reco::deltaPhi(position.phi(), refPhi) < minDeltaPhi)
0064         continue;
0065
0066       //std::cout<<mystrip.zside()<<" "<<mystrip.plane()<<" "<<mystrip.six()<<" "<<mystrip.siy()<<" "<<mystrip.strip()<<" "<<position.eta()<<" "<<position.phi()<<" "<<strip_it->second.energy()<<" "<<strip_it->second.recoFlag()<<" "<<refEta<<" "<<refPhi<<" "<<reco::deltaPhi(position.phi(), refPhi)<<std::endl;
0067
0068       clusterRecHits.push_back(strip_it->second);
0069       used_s->insert(strip_it->first);
0070
0071       x_pos += strip_it->second.energy() * position.x();
0072       y_pos += strip_it->second.energy() * position.y();
0073       z_pos += strip_it->second.energy() * position.z();
0074     }
0075   }
0076
0077   EcalRecHitCollection::iterator it;
0078   double Eclust = 0;
0079
0080   std::vector<std::pair<DetId, float> > usedHits;
0081   for (it = clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
0082     Eclust += it->energy();
0083     usedHits.push_back(std::pair<DetId, float>(it->id(), 1.));
0084   }
0085
0086   if (Eclust > 0.) {
0087     x_pos /= Eclust;
0088     y_pos /= Eclust;
0089     z_pos /= Eclust;
0090   }
0091   Point pos(x_pos, y_pos, z_pos);
0092
0093   reco::PreshowerCluster cluster = reco::PreshowerCluster(Eclust, pos, usedHits, plane);
0094
0095   return cluster;
0096 }
0097
0098 // returns true if the candidate strip fulfills the requirements to be added to the cluster:
0099 bool PreshowerPhiClusterAlgo::goodStrip(RecHitsMap::iterator candidate_it) {
0100   if (used_s->find(candidate_it->first) != used_s->end())
0101     LogTrace("PreShowerPhiClusterAlgo") << " This strip is in use";
0102   if (candidate_it == rechits_map->end())
0103     LogTrace("PreShowerPhiClusterAlgo") << " No such a strip in rechits_map";
0104   if (candidate_it->second.energy() <= esStripEnergyCut_)
0105     LogTrace("PreShowerPhiClusterAlgo") << "Strip energy" << candidate_it->second.energy() << "is below threshold";
0106
0107   // crystal should not be included...
0108   if ((used_s->find(candidate_it->first) != used_s->end()) ||  // ...if it already belongs to a cluster
0109       (candidate_it == rechits_map->end()) ||                  // ...if it corresponds to a hit
0110       (candidate_it->second.energy() <= esStripEnergyCut_))    // ...if it has a negative or zero energy
0111   {
0112     return false;
0113   }
0114   return true;
0115 }