Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:53

0001 #include "RecoEcal/EgammaClusterAlgos/interface/HybridClusterAlgo.h"
0002 #include "RecoCaloTools/Navigation/interface/EcalBarrelNavigator.h"
0003 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0004 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "RecoEcal/EgammaCoreTools/interface/ClusterEtLess.h"
0007 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0008 
0009 #include <iostream>
0010 #include <map>
0011 #include <vector>
0012 #include <set>
0013 
0014 //The real constructor
0015 HybridClusterAlgo::HybridClusterAlgo(double eb_str,
0016                                      int step,
0017                                      double ethres,
0018                                      double eseed,
0019                                      double xi,
0020                                      bool useEtForXi,
0021                                      double ewing,
0022                                      const std::vector<int>& v_chstatus,
0023                                      const PositionCalc& posCalculator,
0024                                      bool dynamicEThres,
0025                                      double eThresA,
0026                                      double eThresB,
0027                                      const std::vector<int>& severityToExclude,
0028                                      bool excludeFromCluster)
0029     :
0030 
0031       eb_st(eb_str),
0032       phiSteps_(step),
0033       eThres_(ethres),
0034       eThresA_(eThresA),
0035       eThresB_(eThresB),
0036       Eseed(eseed),
0037       Xi(xi),
0038       UseEtForXi(useEtForXi),
0039       Ewing(ewing),
0040       dynamicEThres_(dynamicEThres),
0041       v_chstatus_(v_chstatus),
0042       v_severitylevel_(severityToExclude),
0043       //severityRecHitThreshold_(severityRecHitThreshold),
0044       //severitySpikeThreshold_(severitySpikeThreshold),
0045       excludeFromCluster_(excludeFromCluster)
0046 
0047 {
0048   dynamicPhiRoad_ = false;
0049   LogTrace("EcalClusters") << "dynamicEThres: " << dynamicEThres_ << " : A,B " << eThresA_ << ", " << eThresB_;
0050   LogTrace("EcalClusters") << "Eseed: " << Eseed << " , Xi: " << Xi << ", UseEtForXi: " << UseEtForXi;
0051 
0052   //if (dynamicPhiRoad_) phiRoadAlgo_ = new BremRecoveryPhiRoadAlgo(bremRecoveryPset);
0053   posCalculator_ = posCalculator;
0054   topo_ = new EcalBarrelHardcodedTopology();
0055 
0056   std::sort(v_chstatus_.begin(), v_chstatus_.end());
0057 }
0058 
0059 // Return a vector of clusters from a collection of EcalRecHits:
0060 void HybridClusterAlgo::makeClusters(const EcalRecHitCollection* recColl,
0061                                      const CaloSubdetectorGeometry* geometry,
0062                                      reco::BasicClusterCollection& basicClusters,
0063                                      const EcalSeverityLevelAlgo* sevLv,
0064                                      bool regional,
0065                                      const std::vector<RectangularEtaPhiRegion>& regions) {
0066   //clear vector of seeds
0067   seeds.clear();
0068   //clear flagged crystals
0069   excludedCrys_.clear();
0070   //clear map of supercluster/basiccluster association
0071   clustered_.clear();
0072   //clear set of used detids
0073   useddetids.clear();
0074   //clear vector of seed clusters
0075   seedClus_.clear();
0076   //Pass in a pointer to the collection.
0077   recHits_ = recColl;
0078 
0079   LogTrace("EcalClusters") << "Cleared vectors, starting clusterization...";
0080   LogTrace("EcalClusters") << " Purple monkey aardvark.";
0081   //
0082 
0083   int nregions = 0;
0084   if (regional)
0085     nregions = regions.size();
0086 
0087   if (!regional || nregions) {
0088     EcalRecHitCollection::const_iterator it;
0089 
0090     for (it = recHits_->begin(); it != recHits_->end(); it++) {
0091       //Make the vector of seeds that we're going to use.
0092       //One of the few places position is used, needed for ET calculation.
0093       auto this_cell = geometry->getGeometry(it->id());
0094       const GlobalPoint& position = this_cell->getPosition();
0095 
0096       // Require that RecHit is within clustering region in case
0097       // of regional reconstruction
0098       bool withinRegion = false;
0099       if (regional) {
0100         std::vector<RectangularEtaPhiRegion>::const_iterator region;
0101         for (region = regions.begin(); region != regions.end(); region++) {
0102           if (region->inRegion(this_cell->etaPos(), this_cell->phiPos())) {
0103             withinRegion = true;
0104             break;
0105           }
0106         }
0107       }
0108 
0109       if (!regional || withinRegion) {
0110         //Must pass seed threshold
0111         float ET = it->energy() * position.basicVector().unit().perp();
0112 
0113         LogTrace("EcalClusters") << "Seed crystal: ";
0114 
0115         if (ET > eb_st) {
0116           // avoid seeding for anomalous channels
0117           if (!it->checkFlag(EcalRecHit::kGood)) {  // if rechit is good, no need for further checks
0118             if (it->checkFlags(v_chstatus_)) {
0119               if (excludeFromCluster_)
0120                 excludedCrys_.insert(it->id());
0121               continue;  // the recHit has to be excluded from seeding
0122             }
0123           }
0124 
0125           int severityFlag = sevLv->severityLevel(it->id(), *recHits_);
0126           std::vector<int>::const_iterator sit =
0127               std::find(v_severitylevel_.begin(), v_severitylevel_.end(), severityFlag);
0128           LogTrace("EcalClusters") << "found flag: " << severityFlag;
0129 
0130           if (sit != v_severitylevel_.end()) {
0131             if (excludeFromCluster_)
0132               excludedCrys_.insert(it->id());
0133             continue;
0134           }
0135           seeds.push_back(*it);
0136 
0137           LogTrace("EcalClusters") << "Seed ET: " << ET;
0138           LogTrace("EcalClsuters") << "Seed E: " << it->energy();
0139         }
0140       }
0141     }
0142   }
0143 
0144   //Yay sorting.
0145   LogTrace("EcalClusters") << "Built vector of seeds, about to sort them...";
0146 
0147   //Needs three argument sort with seed comparison operator
0148   sort(seeds.begin(), seeds.end(), [](auto const& x, auto const& y) { return x.energy() > y.energy(); });
0149 
0150   LogTrace("EcalClusters") << "done";
0151 
0152   //Now to do the work.
0153   LogTrace("EcalClusters") << "About to call mainSearch...";
0154   mainSearch(recColl, geometry);
0155   LogTrace("EcalClusters") << "done";
0156 
0157   //Hand the basicclusters back to the producer.  It has to
0158   //put them in the event.  Then we can make superclusters.
0159   std::map<int, reco::BasicClusterCollection>::iterator bic;
0160   for (bic = clustered_.begin(); bic != clustered_.end(); bic++) {
0161     reco::BasicClusterCollection bl = bic->second;
0162     for (int j = 0; j < int(bl.size()); ++j) {
0163       basicClusters.push_back(bl[j]);
0164     }
0165   }
0166 
0167   //Yay more sorting.
0168   sort(basicClusters.rbegin(), basicClusters.rend(), isClusterEtLess);
0169   //Done!
0170   LogTrace("EcalClusters") << "returning to producer. ";
0171 }
0172 
0173 void HybridClusterAlgo::mainSearch(const EcalRecHitCollection* hits, const CaloSubdetectorGeometry* geometry) {
0174   LogTrace("EcalClusters") << "HybridClusterAlgo Algorithm - looking for clusters";
0175   LogTrace("EcalClusters") << "Found the following clusters:";
0176 
0177   // Loop over seeds:
0178   std::vector<EcalRecHit>::iterator it;
0179   int clustercounter = 0;
0180 
0181   for (it = seeds.begin(); it != seeds.end(); it++) {
0182     std::vector<reco::BasicCluster> thisseedClusters;
0183     DetId itID = it->id();
0184 
0185     // make sure the current seed has not been used/will not be used in the future:
0186     std::set<DetId>::iterator seed_in_rechits_it = useddetids.find(itID);
0187 
0188     if (seed_in_rechits_it != useddetids.end())
0189       continue;
0190     //If this seed is already used, then don't use it again.
0191 
0192     // output some info on the hit:
0193     LogTrace("EcalClusters") << "*****************************************************"
0194                              << "Seed of energy E = " << it->energy() << " @ " << EBDetId(itID)
0195                              << "*****************************************************";
0196 
0197     //Make a navigator, and set it to the seed cell.
0198     EcalBarrelNavigatorHT navigator(itID, topo_);
0199 
0200     //Now use the navigator to start building dominoes.
0201 
0202     //Walking positive in phi:
0203     std::vector<double> dominoEnergyPhiPlus;                   //Here I will store the results of the domino sums
0204     std::vector<std::vector<EcalRecHit> > dominoCellsPhiPlus;  //These are the actual EcalRecHit for dominos.
0205 
0206     //Walking negative in phi
0207     std::vector<double> dominoEnergyPhiMinus;                   //Here I will store the results of the domino sums
0208     std::vector<std::vector<EcalRecHit> > dominoCellsPhiMinus;  //These are the actual EcalRecHit for dominos.
0209 
0210     //The two sets together.
0211     std::vector<double> dominoEnergy;                   //Here I will store the results of the domino sums
0212     std::vector<std::vector<EcalRecHit> > dominoCells;  //These are the actual EcalRecHit for dominos.
0213 
0214     //First, the domino about the seed:
0215     std::vector<EcalRecHit> initialdomino;
0216     double e_init = makeDomino(navigator, initialdomino);
0217     if (e_init < Eseed)
0218       continue;
0219 
0220     LogTrace("EcalClusters") << "Made initial domino";
0221 
0222     //
0223     // compute the phi road length
0224     double phiSteps;
0225     double et5x5 = et25(navigator, hits, geometry);
0226     if (dynamicPhiRoad_ && e_init > 0) {
0227       phiSteps = phiRoadAlgo_->barrelPhiRoad(et5x5);
0228       navigator.home();
0229     } else
0230       phiSteps = phiSteps_;
0231 
0232     //Positive phi steps.
0233     for (int i = 0; i < phiSteps; ++i) {
0234       //remember, this always increments the current position of the navigator.
0235       DetId centerD = navigator.north();
0236       if (centerD.null())
0237         continue;
0238       LogTrace("EcalClusters") << "Step ++" << i << " @ " << EBDetId(centerD);
0239       EcalBarrelNavigatorHT dominoNav(centerD, topo_);
0240 
0241       //Go get the new domino.
0242       std::vector<EcalRecHit> dcells;
0243       double etemp = makeDomino(dominoNav, dcells);
0244 
0245       //save this information
0246       dominoEnergyPhiPlus.push_back(etemp);
0247       dominoCellsPhiPlus.push_back(dcells);
0248     }
0249     LogTrace("EcalClusters") << "Got positive dominos";
0250     //return to initial position
0251     navigator.home();
0252 
0253     //Negative phi steps.
0254     for (int i = 0; i < phiSteps; ++i) {
0255       //remember, this always decrements the current position of the navigator.
0256       DetId centerD = navigator.south();
0257       if (centerD.null())
0258         continue;
0259 
0260       LogTrace("EcalClusters") << "Step --" << i << " @ " << EBDetId(centerD);
0261       EcalBarrelNavigatorHT dominoNav(centerD, topo_);
0262 
0263       //Go get the new domino.
0264       std::vector<EcalRecHit> dcells;
0265       double etemp = makeDomino(dominoNav, dcells);
0266 
0267       //save this information
0268       dominoEnergyPhiMinus.push_back(etemp);
0269       dominoCellsPhiMinus.push_back(dcells);
0270     }
0271 
0272     LogTrace("EcalClusters") << "Got negative dominos: ";
0273 
0274     //Assemble this information:
0275     for (int i = int(dominoEnergyPhiMinus.size()) - 1; i >= 0; --i) {
0276       dominoEnergy.push_back(dominoEnergyPhiMinus[i]);
0277       dominoCells.push_back(dominoCellsPhiMinus[i]);
0278     }
0279     dominoEnergy.push_back(e_init);
0280     dominoCells.push_back(initialdomino);
0281     for (int i = 0; i < int(dominoEnergyPhiPlus.size()); ++i) {
0282       dominoEnergy.push_back(dominoEnergyPhiPlus[i]);
0283       dominoCells.push_back(dominoCellsPhiPlus[i]);
0284     }
0285 
0286     //Ok, now I have all I need in order to go ahead and make clusters.
0287     LogTrace("EcalClusters") << "Dumping domino energies: ";
0288 
0289     //for (int i=0;i<int(dominoEnergy.size());++i){
0290     //       LogTrace("EcalClusters") << "Domino: " << i << " E: " << dominoEnergy[i] ;
0291     //      }
0292     //Identify the peaks in this set of dominos:
0293     //Peak==a domino whose energy is greater than the two adjacent dominos.
0294     //thus a peak in the local sense.
0295     double etToE(1.);
0296     if (!UseEtForXi) {
0297       etToE = 1. / e2Et(navigator, hits, geometry);
0298     }
0299     std::vector<int> PeakIndex;
0300     for (int i = 1; i < int(dominoEnergy.size()) - 1; ++i) {
0301       if (dominoEnergy[i] > dominoEnergy[i - 1] && dominoEnergy[i] >= dominoEnergy[i + 1] &&
0302           dominoEnergy[i] > sqrt(Eseed * Eseed + Xi * Xi * et5x5 * et5x5 * etToE *
0303                                                      etToE)) {  // Eseed increases ~proportiaonally to et5x5
0304         PeakIndex.push_back(i);
0305       }
0306     }
0307 
0308     LogTrace("EcalClusters") << "Found: " << PeakIndex.size() << " peaks.";
0309 
0310     //Order these peaks by energy:
0311     for (int i = 0; i < int(PeakIndex.size()); ++i) {
0312       for (int j = 0; j < int(PeakIndex.size()) - 1; ++j) {
0313         if (dominoEnergy[PeakIndex[j]] < dominoEnergy[PeakIndex[j + 1]]) {
0314           int ihold = PeakIndex[j + 1];
0315           PeakIndex[j + 1] = PeakIndex[j];
0316           PeakIndex[j] = ihold;
0317         }
0318       }
0319     }
0320 
0321     std::vector<int> OwnerShip;
0322     std::vector<double> LumpEnergy;
0323     OwnerShip.reserve(int(dominoEnergy.size()));
0324     for (int i = 0; i < int(dominoEnergy.size()); ++i)
0325       OwnerShip.push_back(-1);
0326 
0327     //Loop over peaks.
0328     double eThres = eThres_;
0329     double e5x5 = 0.0;
0330     for (int i = 0; i < int(PeakIndex.size()); ++i) {
0331       int idxPeak = PeakIndex[i];
0332       OwnerShip[idxPeak] = i;
0333       double lump = dominoEnergy[idxPeak];
0334 
0335       // compute eThres for this peak
0336       // if set to dynamic (otherwise uncanged from
0337       // fixed setting
0338       if (dynamicEThres_) {
0339         //std::cout << "i : " << i << " idxPeak " << idxPeak << std::endl;
0340         //std::cout << "    the dominoEnergy.size() = " << dominoEnergy.size() << std::endl;
0341         // compute e5x5 for this seed crystal
0342         //std::cout << "idxPeak, phiSteps " << idxPeak << ", " << phiSteps << std::endl;
0343         e5x5 = lump;
0344         //std::cout << "lump " << e5x5 << std::endl;
0345         if ((idxPeak + 1) < (int)dominoEnergy.size())
0346           e5x5 += dominoEnergy[idxPeak + 1];
0347         //std::cout << "+1 " << e5x5 << std::endl;
0348         if ((idxPeak + 2) < (int)dominoEnergy.size())
0349           e5x5 += dominoEnergy[idxPeak + 2];
0350         //std::cout << "+2 " << e5x5 << std::endl;
0351         if ((idxPeak - 1) > 0)
0352           e5x5 += dominoEnergy[idxPeak - 1];
0353         //std::cout << "-1 " << e5x5 << std::endl;
0354         if ((idxPeak - 2) > 0)
0355           e5x5 += dominoEnergy[idxPeak - 2];
0356         //std::cout << "-2 " << e5x5 << std::endl;
0357         // compute eThres
0358         eThres = (eThresA_ * e5x5) + eThresB_;
0359         //std::cout << eThres << std::endl;
0360         //std::cout << std::endl;
0361       }
0362 
0363       //Loop over adjacent dominos at higher phi
0364       for (int j = idxPeak + 1; j < int(dominoEnergy.size()); ++j) {
0365         if (OwnerShip[j] == -1 && dominoEnergy[j] > eThres && dominoEnergy[j] <= dominoEnergy[j - 1]) {
0366           OwnerShip[j] = i;
0367           lump += dominoEnergy[j];
0368         } else {
0369           break;
0370         }
0371       }
0372       //loop over adjacent dominos at lower phi.  Sum up energy of lumps.
0373       for (int j = idxPeak - 1; j >= 0; --j) {
0374         if (OwnerShip[j] == -1 && dominoEnergy[j] > eThres && dominoEnergy[j] <= dominoEnergy[j + 1]) {
0375           OwnerShip[j] = i;
0376           lump += dominoEnergy[j];
0377         } else {
0378           break;
0379         }
0380       }
0381       LumpEnergy.push_back(lump);
0382     }
0383 
0384     //Make the basic clusters:
0385     for (int i = 0; i < int(PeakIndex.size()); ++i) {
0386       bool HasSeedCrystal = false;
0387       //One cluster for each peak.
0388       std::vector<EcalRecHit> recHits;
0389       std::vector<std::pair<DetId, float> > dets;
0390       int nhits = 0;
0391       for (int j = 0; j < int(dominoEnergy.size()); ++j) {
0392         if (OwnerShip[j] == i) {
0393           std::vector<EcalRecHit> temp = dominoCells[j];
0394           for (int k = 0; k < int(temp.size()); ++k) {
0395             dets.push_back(std::pair<DetId, float>(temp[k].id(), 1.));  // by default energy fractions are 1
0396             if (temp[k].id() == itID)
0397               HasSeedCrystal = true;
0398             recHits.push_back(temp[k]);
0399             nhits++;
0400           }
0401         }
0402       }
0403       LogTrace("EcalClusters") << "Adding a cluster with: " << nhits;
0404       LogTrace("EcalClusters") << "total E: " << LumpEnergy[i];
0405       LogTrace("EcalClusters") << "total dets: " << dets.size();
0406 
0407       //Get Calorimeter position
0408       Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
0409 
0410       //double totChi2=0;
0411       //double totE=0;
0412       std::vector<std::pair<DetId, float> > usedHits;
0413       for (int blarg = 0; blarg < int(recHits.size()); ++blarg) {
0414         //totChi2 +=0;
0415         //totE+=recHits[blarg].energy();
0416         usedHits.push_back(std::make_pair<DetId, float>(recHits[blarg].id(), 1.0));
0417         useddetids.insert(recHits[blarg].id());
0418       }
0419 
0420       //if (totE>0)
0421       //totChi2/=totE;
0422 
0423       if (HasSeedCrystal) {
0424         // note that this "basiccluster" has the seed crystal of the hyrbid, so record it
0425         seedClus_.push_back(reco::BasicCluster(
0426             LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid, itID));
0427         // and also add to the vector of clusters that will be used in constructing
0428         // the supercluster
0429         thisseedClusters.push_back(reco::BasicCluster(
0430             LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid, itID));
0431       } else {
0432         // note that if this "basiccluster" is not the one that seeded the hybrid,
0433         // the seed crystal is unset in this entry in the vector of clusters that will
0434         // be used in constructing the super cluster
0435         thisseedClusters.push_back(reco::BasicCluster(
0436             LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid));
0437       }
0438     }
0439 
0440     // Make association so that superclusters can be made later.
0441     // but only if some BasicClusters have been found...
0442     if (!thisseedClusters.empty()) {
0443       clustered_.insert(std::make_pair(clustercounter, thisseedClusters));
0444       clustercounter++;
0445     }
0446   }  //Seed loop
0447 
0448 }  // end of mainSearch
0449 
0450 reco::SuperClusterCollection HybridClusterAlgo::makeSuperClusters(const reco::CaloClusterPtrVector& clustersCollection) {
0451   //Here's what we'll return.
0452   reco::SuperClusterCollection SCcoll;
0453 
0454   //Here's our map iterator that gives us the appropriate association.
0455   std::map<int, reco::BasicClusterCollection>::iterator mapit;
0456   for (mapit = clustered_.begin(); mapit != clustered_.end(); mapit++) {
0457     reco::CaloClusterPtrVector thissc;
0458     reco::CaloClusterPtr seed;  //This is not really a seed, but I need to tell SuperCluster something.
0459     //So I choose the highest energy basiccluster in the SuperCluster.
0460 
0461     std::vector<reco::BasicCluster> thiscoll = mapit->second;  //This is the set of BasicClusters in this
0462     //SuperCluster
0463 
0464     double ClusterE = 0;  //Sum of cluster energies for supercluster.
0465     //Holders for position of this supercluster.
0466     double posX = 0;
0467     double posY = 0;
0468     double posZ = 0;
0469 
0470     //Loop over this set of basic clusters, find their references, and add them to the
0471     //supercluster.  This could be somehow more efficient.
0472 
0473     for (int i = 0; i < int(thiscoll.size()); ++i) {
0474       reco::BasicCluster thisclus = thiscoll[i];  //The Cluster in question.
0475       for (int j = 0; j < int(clustersCollection.size()); ++j) {
0476         //Find the appropriate cluster from the list of references
0477         reco::BasicCluster cluster_p = *clustersCollection[j];
0478         if (thisclus == cluster_p) {  //Comparison based on energy right now.
0479           thissc.push_back(clustersCollection[j]);
0480           bool isSeed = false;
0481           for (int qu = 0; qu < int(seedClus_.size()); ++qu) {
0482             if (cluster_p == seedClus_[qu])
0483               isSeed = true;
0484           }
0485           if (isSeed)
0486             seed = clustersCollection[j];
0487 
0488           ClusterE += cluster_p.energy();
0489           posX += cluster_p.energy() * cluster_p.position().X();
0490           posY += cluster_p.energy() * cluster_p.position().Y();
0491           posZ += cluster_p.energy() * cluster_p.position().Z();
0492         }
0493       }  //End loop over finding references.
0494     }    //End loop over clusters.
0495 
0496     posX /= ClusterE;
0497     posY /= ClusterE;
0498     posZ /= ClusterE;
0499 
0500     /* //This part is moved to EgammaSCEnergyCorrectionAlgo
0501            double preshowerE = 0.;
0502            double phiWidth = 0.;
0503            double etaWidth = 0.;
0504         //Calculate phiWidth & etaWidth for SuperClusters
0505         reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc, preshowerE, phiWidth, etaWidth);
0506         SCShape_->Calculate_Covariances(suCl);
0507         phiWidth = SCShape_->phiWidth();
0508         etaWidth = SCShape_->etaWidth();
0509         //Assign phiWidth & etaWidth to SuperCluster as data members
0510         suCl.setPhiWidth(phiWidth);
0511         suCl.setEtaWidth(etaWidth);
0512          */
0513 
0514     reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc);
0515     SCcoll.push_back(suCl);
0516 
0517     LogTrace("EcalClusters") << "Super cluster sum: " << ClusterE;
0518     LogTrace("EcalClusters") << "Made supercluster with energy E: " << suCl.energy();
0519 
0520   }  //end loop over map
0521   sort(SCcoll.rbegin(), SCcoll.rend(), isClusterEtLess);
0522   return SCcoll;
0523 }
0524 
0525 double HybridClusterAlgo::makeDomino(EcalBarrelNavigatorHT& navigator, std::vector<EcalRecHit>& cells) {
0526   //At the beginning of this function, the navigator starts at the middle of the domino,
0527   //and that's where EcalBarrelNavigator::home() should send it.
0528   //Walk one crystal in eta to either side of the initial point.  Sum the three cell energies.
0529   //If the resultant energy is larger than Ewing, then go an additional cell to either side.
0530   //Returns:  Total domino energy.  Also, stores the cells used to create domino in the vector.
0531   cells.clear();
0532   double Etot = 0;
0533 
0534   //Ready?  Get the starting cell.
0535   DetId center = navigator.pos();
0536   EcalRecHitCollection::const_iterator center_it = recHits_->find(center);
0537 
0538   if (center_it != recHits_->end()) {
0539     EcalRecHit SeedHit = *center_it;
0540     if (useddetids.find(center) == useddetids.end() && excludedCrys_.find(center) == excludedCrys_.end()) {
0541       Etot += SeedHit.energy();
0542       cells.push_back(SeedHit);
0543     }
0544   }
0545   //One step upwards in Ieta:
0546   DetId ieta1 = navigator.west();
0547   EcalRecHitCollection::const_iterator eta1_it = recHits_->find(ieta1);
0548   if (eta1_it != recHits_->end()) {
0549     EcalRecHit UpEta = *eta1_it;
0550     if (useddetids.find(ieta1) == useddetids.end() && excludedCrys_.find(ieta1) == excludedCrys_.end()) {
0551       Etot += UpEta.energy();
0552       cells.push_back(UpEta);
0553     }
0554   }
0555 
0556   //Go back to the middle.
0557   navigator.home();
0558 
0559   //One step downwards in Ieta:
0560   DetId ieta2 = navigator.east();
0561   EcalRecHitCollection::const_iterator eta2_it = recHits_->find(ieta2);
0562   if (eta2_it != recHits_->end()) {
0563     EcalRecHit DownEta = *eta2_it;
0564     if (useddetids.find(ieta2) == useddetids.end() && excludedCrys_.find(ieta2) == excludedCrys_.end()) {
0565       Etot += DownEta.energy();
0566       cells.push_back(DownEta);
0567     }
0568   }
0569 
0570   //Now check the energy.  If smaller than Ewing, then we're done.  If greater than Ewing, we have to
0571   //add two additional cells, the 'wings'
0572   if (Etot < Ewing) {
0573     navigator.home();  //Needed even here!!
0574     return Etot;       //Done!  Not adding 'wings'.
0575   }
0576 
0577   //Add the extra 'wing' cells.  Remember, we haven't sent the navigator home,
0578   //we're still on the DownEta cell.
0579   if (ieta2 != DetId(0)) {
0580     DetId ieta3 = navigator.east();  //Take another step downward.
0581     EcalRecHitCollection::const_iterator eta3_it = recHits_->find(ieta3);
0582     if (eta3_it != recHits_->end()) {
0583       EcalRecHit DownEta2 = *eta3_it;
0584       if (useddetids.find(ieta3) == useddetids.end() && excludedCrys_.find(ieta3) == excludedCrys_.end()) {
0585         Etot += DownEta2.energy();
0586         cells.push_back(DownEta2);
0587       }
0588     }
0589   }
0590   //Now send the navigator home.
0591   navigator.home();
0592   navigator.west();  //Now you're on eta1_it
0593   if (ieta1 != DetId(0)) {
0594     DetId ieta4 = navigator.west();  //Take another step upward.
0595     EcalRecHitCollection::const_iterator eta4_it = recHits_->find(ieta4);
0596     if (eta4_it != recHits_->end()) {
0597       EcalRecHit UpEta2 = *eta4_it;
0598       if (useddetids.find(ieta4) == useddetids.end() && excludedCrys_.find(ieta4) == excludedCrys_.end()) {
0599         Etot += UpEta2.energy();
0600         cells.push_back(UpEta2);
0601       }
0602     }
0603   }
0604   navigator.home();
0605   return Etot;
0606 }
0607 
0608 double HybridClusterAlgo::et25(EcalBarrelNavigatorHT& navigator,
0609                                const EcalRecHitCollection* hits,
0610                                const CaloSubdetectorGeometry* geometry) {
0611   DetId thisDet;
0612   std::vector<std::pair<DetId, float> > dets;
0613   dets.clear();
0614   EcalRecHitCollection::const_iterator hit;
0615   double energySum = 0.0;
0616 
0617   for (int dx = -2; dx < 3; ++dx) {
0618     for (int dy = -2; dy < 3; ++dy) {
0619       //std::cout << "dx, dy " << dx << ", " << dy << std::endl;
0620       thisDet = navigator.offsetBy(dx, dy);
0621       navigator.home();
0622 
0623       if (thisDet != DetId(0)) {
0624         hit = recHits_->find(thisDet);
0625         if (hit != recHits_->end()) {
0626           dets.push_back(std::pair<DetId, float>(thisDet, 1.));  // by default hit energy fraction is set to 1
0627           energySum += hit->energy();
0628         }
0629       }
0630     }
0631   }
0632 
0633   // convert it to ET
0634   //std::cout << "dets.size(), energySum: " << dets.size() << ", " << energySum << std::endl;
0635   Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
0636   double et = energySum / cosh(pos.eta());
0637   return et;
0638 }
0639 
0640 double HybridClusterAlgo::e2Et(EcalBarrelNavigatorHT& navigator,
0641                                const EcalRecHitCollection* hits,
0642                                const CaloSubdetectorGeometry* geometry) {
0643   DetId thisDet;
0644   std::vector<std::pair<DetId, float> > dets;
0645   dets.clear();
0646   EcalRecHitCollection::const_iterator hit;
0647 
0648   for (int dx = -2; dx < 3; ++dx) {
0649     for (int dy = -2; dy < 3; ++dy) {
0650       thisDet = navigator.offsetBy(dx, dy);
0651       navigator.home();
0652 
0653       if (thisDet != DetId(0)) {
0654         hit = recHits_->find(thisDet);
0655         if (hit != recHits_->end()) {
0656           dets.push_back(std::pair<DetId, float>(thisDet, 1.));  // by default hit energy fraction is set to 1
0657         }
0658       }
0659     }
0660   }
0661 
0662   // compute coefficient to turn E into Et
0663   Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
0664   return 1 / cosh(pos.eta());
0665 }