Project CMSSW displayed by LXR CMSSW_14_2_X_2024-11-08-2300/​RecoEcal/​EgammaClusterAlgos/​src/​HybridClusterAlgo.cc	Source navigation Diff markup Identifier search General search
	Version: CMSSW_14_2_X_2024-11-08-2300 CMSSW_14_2_X_2024-11-06-2300 CMSSW_14_2_X_2024-11-05-2300 CMSSW_14_2_X_2024-11-04-2300 CMSSW_14_2_X_2024-11-03-2300 Revert   Change

File indexing completed on 2024-09-07 04:37:26

 #include "RecoEcal/EgammaClusterAlgos/interface/HybridClusterAlgo.h"
 #include "RecoCaloTools/Navigation/interface/EcalBarrelNavigator.h"
 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "RecoEcal/EgammaCoreTools/interface/ClusterEtLess.h"
 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
 
 #include <iostream>
 #include <map>
 #include <vector>
 #include <set>
 
 //The real constructor
 HybridClusterAlgo::HybridClusterAlgo(double eb_str,
                                      int step,
                                      double ethres,
                                      double eseed,
                                      double xi,
                                      bool useEtForXi,
                                      double ewing,
                                      const std::vector<int>& v_chstatus,
                                      const PositionCalc& posCalculator,
                                      bool dynamicEThres,
                                      double eThresA,
                                      double eThresB,
                                      const std::vector<int>& severityToExclude,
                                      bool excludeFromCluster)
     :
 
       eb_st(eb_str),
       phiSteps_(step),
       eThres_(ethres),
       eThresA_(eThresA),
       eThresB_(eThresB),
       Eseed(eseed),
       Xi(xi),
       UseEtForXi(useEtForXi),
       Ewing(ewing),
       dynamicEThres_(dynamicEThres),
       v_chstatus_(v_chstatus),
       v_severitylevel_(severityToExclude),
       //severityRecHitThreshold_(severityRecHitThreshold),
       //severitySpikeThreshold_(severitySpikeThreshold),
       excludeFromCluster_(excludeFromCluster)
 
 {
   dynamicPhiRoad_ = false;
   LogTrace("EcalClusters") << "dynamicEThres: " << dynamicEThres_ << " : A,B " << eThresA_ << ", " << eThresB_;
   LogTrace("EcalClusters") << "Eseed: " << Eseed << " , Xi: " << Xi << ", UseEtForXi: " << UseEtForXi;
 
   //if (dynamicPhiRoad_) phiRoadAlgo_ = new BremRecoveryPhiRoadAlgo(bremRecoveryPset);
   posCalculator_ = posCalculator;
   topo_ = new EcalBarrelHardcodedTopology();
 
   std::sort(v_chstatus_.begin(), v_chstatus_.end());
 }
 
 // Return a vector of clusters from a collection of EcalRecHits:
 void HybridClusterAlgo::makeClusters(const EcalRecHitCollection* recColl,
                                      const CaloSubdetectorGeometry* geometry,
                                      reco::BasicClusterCollection& basicClusters,
                                      const EcalSeverityLevelAlgo* sevLv,
                                      bool regional,
                                      const std::vector<RectangularEtaPhiRegion>& regions) {
   //clear vector of seeds
   seeds.clear();
   //clear flagged crystals
   excludedCrys_.clear();
   //clear map of supercluster/basiccluster association
   clustered_.clear();
   //clear set of used detids
   useddetids.clear();
   //clear vector of seed clusters
   seedClus_.clear();
   //Pass in a pointer to the collection.
   recHits_ = recColl;
 
   LogTrace("EcalClusters") << "Cleared vectors, starting clusterization...";
   LogTrace("EcalClusters") << " Purple monkey aardvark.";
   //
 
   int nregions = 0;
   if (regional)
     nregions = regions.size();
 
   if (!regional || nregions) {
     EcalRecHitCollection::const_iterator it;
 
     for (it = recHits_->begin(); it != recHits_->end(); it++) {
       //Make the vector of seeds that we're going to use.
       //One of the few places position is used, needed for ET calculation.
       auto this_cell = geometry->getGeometry(it->id());
       const GlobalPoint& position = this_cell->getPosition();
 
       // Require that RecHit is within clustering region in case
       // of regional reconstruction
       bool withinRegion = false;
       if (regional) {
         std::vector<RectangularEtaPhiRegion>::const_iterator region;
         for (region = regions.begin(); region != regions.end(); region++) {
           if (region->inRegion(this_cell->etaPos(), this_cell->phiPos())) {
             withinRegion = true;
             break;
           }
         }
       }
 
       if (!regional || withinRegion) {
         //Must pass seed threshold
         float ET = it->energy() * position.basicVector().unit().perp();
 
         LogTrace("EcalClusters") << "Seed crystal: ";
 
         if (ET > eb_st) {
           // avoid seeding for anomalous channels
           if (!it->checkFlag(EcalRecHit::kGood)) {  // if rechit is good, no need for further checks
             if (it->checkFlags(v_chstatus_)) {
               if (excludeFromCluster_)
                 excludedCrys_.insert(it->id());
               continue;  // the recHit has to be excluded from seeding
             }
           }
 
           int severityFlag = sevLv->severityLevel(it->id(), *recHits_);
           std::vector<int>::const_iterator sit =
               std::find(v_severitylevel_.begin(), v_severitylevel_.end(), severityFlag);
           LogTrace("EcalClusters") << "found flag: " << severityFlag;
 
           if (sit != v_severitylevel_.end()) {
             if (excludeFromCluster_)
               excludedCrys_.insert(it->id());
             continue;
           }
           seeds.push_back(*it);
 
           LogTrace("EcalClusters") << "Seed ET: " << ET;
           LogTrace("EcalClsuters") << "Seed E: " << it->energy();
         }
       }
     }
   }
 
   //Yay sorting.
   LogTrace("EcalClusters") << "Built vector of seeds, about to sort them...";
 
   //Needs three argument sort with seed comparison operator
   sort(seeds.begin(), seeds.end(), [](auto const& x, auto const& y) { return x.energy() > y.energy(); });
 
   LogTrace("EcalClusters") << "done";
 
   //Now to do the work.
   LogTrace("EcalClusters") << "About to call mainSearch...";
   mainSearch(recColl, geometry);
   LogTrace("EcalClusters") << "done";
 
   //Hand the basicclusters back to the producer.  It has to
   //put them in the event.  Then we can make superclusters.
   std::map<int, reco::BasicClusterCollection>::iterator bic;
   for (bic = clustered_.begin(); bic != clustered_.end(); bic++) {
     reco::BasicClusterCollection bl = bic->second;
     for (int j = 0; j < int(bl.size()); ++j) {
       basicClusters.push_back(bl[j]);
     }
   }
 
   //Yay more sorting.
   sort(basicClusters.rbegin(), basicClusters.rend(), isClusterEtLess);
   //Done!
   LogTrace("EcalClusters") << "returning to producer. ";
 }
 
 void HybridClusterAlgo::mainSearch(const EcalRecHitCollection* hits, const CaloSubdetectorGeometry* geometry) {
   LogTrace("EcalClusters") << "HybridClusterAlgo Algorithm - looking for clusters";
   LogTrace("EcalClusters") << "Found the following clusters:";
 
   // Loop over seeds:
   std::vector<EcalRecHit>::iterator it;
   int clustercounter = 0;
 
   for (it = seeds.begin(); it != seeds.end(); it++) {
     std::vector<reco::BasicCluster> thisseedClusters;
     DetId itID = it->id();
 
     // make sure the current seed has not been used/will not be used in the future:
     std::set<DetId>::iterator seed_in_rechits_it = useddetids.find(itID);
 
     if (seed_in_rechits_it != useddetids.end())
       continue;
     //If this seed is already used, then don't use it again.
 
     // output some info on the hit:
     LogTrace("EcalClusters") << "*****************************************************"
                              << "Seed of energy E = " << it->energy() << " @ " << EBDetId(itID)
                              << "*****************************************************";
 
     //Make a navigator, and set it to the seed cell.
     EcalBarrelNavigatorHT navigator(itID, topo_);
 
     //Now use the navigator to start building dominoes.
 
     //Walking positive in phi:
     std::vector<double> dominoEnergyPhiPlus;                   //Here I will store the results of the domino sums
     std::vector<std::vector<EcalRecHit> > dominoCellsPhiPlus;  //These are the actual EcalRecHit for dominos.
 
     //Walking negative in phi
     std::vector<double> dominoEnergyPhiMinus;                   //Here I will store the results of the domino sums
     std::vector<std::vector<EcalRecHit> > dominoCellsPhiMinus;  //These are the actual EcalRecHit for dominos.
 
     //The two sets together.
     std::vector<double> dominoEnergy;                   //Here I will store the results of the domino sums
     std::vector<std::vector<EcalRecHit> > dominoCells;  //These are the actual EcalRecHit for dominos.
 
     //First, the domino about the seed:
     std::vector<EcalRecHit> initialdomino;
     double e_init = makeDomino(navigator, initialdomino);
     if (e_init < Eseed)
       continue;
 
     LogTrace("EcalClusters") << "Made initial domino";
 
     //
     // compute the phi road length
     double phiSteps;
     double et5x5 = et25(navigator, hits, geometry);
     if (dynamicPhiRoad_ && e_init > 0) {
       phiSteps = phiRoadAlgo_->barrelPhiRoad(et5x5);
       navigator.home();
     } else
       phiSteps = phiSteps_;
 
     //Positive phi steps.
     for (int i = 0; i < phiSteps; ++i) {
       //remember, this always increments the current position of the navigator.
       DetId centerD = navigator.north();
       if (centerD.null())
         continue;
       LogTrace("EcalClusters") << "Step ++" << i << " @ " << EBDetId(centerD);
       EcalBarrelNavigatorHT dominoNav(centerD, topo_);
 
       //Go get the new domino.
       std::vector<EcalRecHit> dcells;
       double etemp = makeDomino(dominoNav, dcells);
 
       //save this information
       dominoEnergyPhiPlus.push_back(etemp);
       dominoCellsPhiPlus.push_back(dcells);
     }
     LogTrace("EcalClusters") << "Got positive dominos";
     //return to initial position
     navigator.home();
 
     //Negative phi steps.
     for (int i = 0; i < phiSteps; ++i) {
       //remember, this always decrements the current position of the navigator.
       DetId centerD = navigator.south();
       if (centerD.null())
         continue;
 
       LogTrace("EcalClusters") << "Step --" << i << " @ " << EBDetId(centerD);
       EcalBarrelNavigatorHT dominoNav(centerD, topo_);
 
       //Go get the new domino.
       std::vector<EcalRecHit> dcells;
       double etemp = makeDomino(dominoNav, dcells);
 
       //save this information
       dominoEnergyPhiMinus.push_back(etemp);
       dominoCellsPhiMinus.push_back(dcells);
     }
 
     LogTrace("EcalClusters") << "Got negative dominos: ";
 
     //Assemble this information:
     for (int i = int(dominoEnergyPhiMinus.size()) - 1; i >= 0; --i) {
       dominoEnergy.push_back(dominoEnergyPhiMinus[i]);
       dominoCells.push_back(dominoCellsPhiMinus[i]);
     }
     dominoEnergy.push_back(e_init);
     dominoCells.push_back(initialdomino);
     for (int i = 0; i < int(dominoEnergyPhiPlus.size()); ++i) {
       dominoEnergy.push_back(dominoEnergyPhiPlus[i]);
       dominoCells.push_back(dominoCellsPhiPlus[i]);
     }
 
     //Ok, now I have all I need in order to go ahead and make clusters.
     LogTrace("EcalClusters") << "Dumping domino energies: ";
 
     //for (int i=0;i<int(dominoEnergy.size());++i){
     //       LogTrace("EcalClusters") << "Domino: " << i << " E: " << dominoEnergy[i] ;
     //      }
     //Identify the peaks in this set of dominos:
     //Peak==a domino whose energy is greater than the two adjacent dominos.
     //thus a peak in the local sense.
     double etToE(1.);
     if (!UseEtForXi) {
       etToE = 1. / e2Et(navigator, hits, geometry);
     }
     std::vector<int> PeakIndex;
     for (int i = 1; i < int(dominoEnergy.size()) - 1; ++i) {
       if (dominoEnergy[i] > dominoEnergy[i - 1] && dominoEnergy[i] >= dominoEnergy[i + 1] &&
           dominoEnergy[i] > sqrt(Eseed * Eseed + Xi * Xi * et5x5 * et5x5 * etToE *
                                                      etToE)) {  // Eseed increases ~proportiaonally to et5x5
         PeakIndex.push_back(i);
       }
     }
 
     LogTrace("EcalClusters") << "Found: " << PeakIndex.size() << " peaks.";
 
     //Order these peaks by energy:
     for (int i = 0; i < int(PeakIndex.size()); ++i) {
       for (int j = 0; j < int(PeakIndex.size()) - 1; ++j) {
         if (dominoEnergy[PeakIndex[j]] < dominoEnergy[PeakIndex[j + 1]]) {
           int ihold = PeakIndex[j + 1];
           PeakIndex[j + 1] = PeakIndex[j];
           PeakIndex[j] = ihold;
         }
       }
     }
 
     std::vector<int> OwnerShip;
     std::vector<double> LumpEnergy;
     OwnerShip.reserve(int(dominoEnergy.size()));
     for (int i = 0; i < int(dominoEnergy.size()); ++i)
       OwnerShip.push_back(-1);
 
     //Loop over peaks.
     double eThres = eThres_;
     double e5x5 = 0.0;
     for (int i = 0; i < int(PeakIndex.size()); ++i) {
       int idxPeak = PeakIndex[i];
       OwnerShip[idxPeak] = i;
       double lump = dominoEnergy[idxPeak];
 
       // compute eThres for this peak
       // if set to dynamic (otherwise uncanged from
       // fixed setting
       if (dynamicEThres_) {
         //std::cout << "i : " << i << " idxPeak " << idxPeak << std::endl;
         //std::cout << "    the dominoEnergy.size() = " << dominoEnergy.size() << std::endl;
         // compute e5x5 for this seed crystal
         //std::cout << "idxPeak, phiSteps " << idxPeak << ", " << phiSteps << std::endl;
         e5x5 = lump;
         //std::cout << "lump " << e5x5 << std::endl;
         if ((idxPeak + 1) < (int)dominoEnergy.size())
           e5x5 += dominoEnergy[idxPeak + 1];
         //std::cout << "+1 " << e5x5 << std::endl;
         if ((idxPeak + 2) < (int)dominoEnergy.size())
           e5x5 += dominoEnergy[idxPeak + 2];
         //std::cout << "+2 " << e5x5 << std::endl;
         if ((idxPeak - 1) > 0)
           e5x5 += dominoEnergy[idxPeak - 1];
         //std::cout << "-1 " << e5x5 << std::endl;
         if ((idxPeak - 2) > 0)
           e5x5 += dominoEnergy[idxPeak - 2];
         //std::cout << "-2 " << e5x5 << std::endl;
         // compute eThres
         eThres = (eThresA_ * e5x5) + eThresB_;
         //std::cout << eThres << std::endl;
         //std::cout << std::endl;
       }
 
       //Loop over adjacent dominos at higher phi
       for (int j = idxPeak + 1; j < int(dominoEnergy.size()); ++j) {
         if (OwnerShip[j] == -1 && dominoEnergy[j] > eThres && dominoEnergy[j] <= dominoEnergy[j - 1]) {
           OwnerShip[j] = i;
           lump += dominoEnergy[j];
         } else {
           break;
         }
       }
       //loop over adjacent dominos at lower phi.  Sum up energy of lumps.
       for (int j = idxPeak - 1; j >= 0; --j) {
         if (OwnerShip[j] == -1 && dominoEnergy[j] > eThres && dominoEnergy[j] <= dominoEnergy[j + 1]) {
           OwnerShip[j] = i;
           lump += dominoEnergy[j];
         } else {
           break;
         }
       }
       LumpEnergy.push_back(lump);
     }
 
     //Make the basic clusters:
     for (int i = 0; i < int(PeakIndex.size()); ++i) {
       bool HasSeedCrystal = false;
       //One cluster for each peak.
       std::vector<EcalRecHit> recHits;
       std::vector<std::pair<DetId, float> > dets;
       int nhits = 0;
       for (int j = 0; j < int(dominoEnergy.size()); ++j) {
         if (OwnerShip[j] == i) {
           std::vector<EcalRecHit> temp = dominoCells[j];
           for (int k = 0; k < int(temp.size()); ++k) {
             dets.push_back(std::pair<DetId, float>(temp[k].id(), 1.));  // by default energy fractions are 1
             if (temp[k].id() == itID)
               HasSeedCrystal = true;
             recHits.push_back(temp[k]);
             nhits++;
           }
         }
       }
       LogTrace("EcalClusters") << "Adding a cluster with: " << nhits;
       LogTrace("EcalClusters") << "total E: " << LumpEnergy[i];
       LogTrace("EcalClusters") << "total dets: " << dets.size();
 
       //Get Calorimeter position
       Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
 
       //double totChi2=0;
       //double totE=0;
       std::vector<std::pair<DetId, float> > usedHits;
       for (int blarg = 0; blarg < int(recHits.size()); ++blarg) {
         //totChi2 +=0;
         //totE+=recHits[blarg].energy();
         usedHits.push_back(std::make_pair<DetId, float>(recHits[blarg].id(), 1.0));
         useddetids.insert(recHits[blarg].id());
       }
 
       //if (totE>0)
       //totChi2/=totE;
 
       if (HasSeedCrystal) {
         // note that this "basiccluster" has the seed crystal of the hyrbid, so record it
         seedClus_.push_back(reco::BasicCluster(
             LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid, itID));
         // and also add to the vector of clusters that will be used in constructing
         // the supercluster
         thisseedClusters.push_back(reco::BasicCluster(
             LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid, itID));
       } else {
         // note that if this "basiccluster" is not the one that seeded the hybrid,
         // the seed crystal is unset in this entry in the vector of clusters that will
         // be used in constructing the super cluster
         thisseedClusters.push_back(reco::BasicCluster(
             LumpEnergy[i], pos, reco::CaloID(reco::CaloID::DET_ECAL_BARREL), usedHits, reco::CaloCluster::hybrid));
       }
     }
 
     // Make association so that superclusters can be made later.
     // but only if some BasicClusters have been found...
     if (!thisseedClusters.empty()) {
       clustered_.insert(std::make_pair(clustercounter, thisseedClusters));
       clustercounter++;
     }
   }  //Seed loop
 
 }  // end of mainSearch
 
 reco::SuperClusterCollection HybridClusterAlgo::makeSuperClusters(const reco::CaloClusterPtrVector& clustersCollection) {
   //Here's what we'll return.
   reco::SuperClusterCollection SCcoll;
 
   //Here's our map iterator that gives us the appropriate association.
   std::map<int, reco::BasicClusterCollection>::iterator mapit;
   for (mapit = clustered_.begin(); mapit != clustered_.end(); mapit++) {
     reco::CaloClusterPtrVector thissc;
     reco::CaloClusterPtr seed;  //This is not really a seed, but I need to tell SuperCluster something.
     //So I choose the highest energy basiccluster in the SuperCluster.
 
     std::vector<reco::BasicCluster> thiscoll = mapit->second;  //This is the set of BasicClusters in this
     //SuperCluster
 
     double ClusterE = 0;  //Sum of cluster energies for supercluster.
     //Holders for position of this supercluster.
     double posX = 0;
     double posY = 0;
     double posZ = 0;
 
     //Loop over this set of basic clusters, find their references, and add them to the
     //supercluster.  This could be somehow more efficient.
 
     for (int i = 0; i < int(thiscoll.size()); ++i) {
       reco::BasicCluster thisclus = thiscoll[i];  //The Cluster in question.
       for (int j = 0; j < int(clustersCollection.size()); ++j) {
         //Find the appropriate cluster from the list of references
         reco::BasicCluster cluster_p = *clustersCollection[j];
         if (thisclus == cluster_p) {  //Comparison based on energy right now.
           thissc.push_back(clustersCollection[j]);
           bool isSeed = false;
           for (int qu = 0; qu < int(seedClus_.size()); ++qu) {
             if (cluster_p == seedClus_[qu])
               isSeed = true;
           }
           if (isSeed)
             seed = clustersCollection[j];
 
           ClusterE += cluster_p.energy();
           posX += cluster_p.energy() * cluster_p.position().X();
           posY += cluster_p.energy() * cluster_p.position().Y();
           posZ += cluster_p.energy() * cluster_p.position().Z();
         }
       }  //End loop over finding references.
     }  //End loop over clusters.
 
     posX /= ClusterE;
     posY /= ClusterE;
     posZ /= ClusterE;
 
     /* //This part is moved to EgammaSCEnergyCorrectionAlgo
            double preshowerE = 0.;
            double phiWidth = 0.;
            double etaWidth = 0.;
         //Calculate phiWidth & etaWidth for SuperClusters
         reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc, preshowerE, phiWidth, etaWidth);
         SCShape_->Calculate_Covariances(suCl);
         phiWidth = SCShape_->phiWidth();
         etaWidth = SCShape_->etaWidth();
         //Assign phiWidth & etaWidth to SuperCluster as data members
         suCl.setPhiWidth(phiWidth);
         suCl.setEtaWidth(etaWidth);
          */
 
     reco::SuperCluster suCl(ClusterE, math::XYZPoint(posX, posY, posZ), seed, thissc);
     SCcoll.push_back(suCl);
 
     LogTrace("EcalClusters") << "Super cluster sum: " << ClusterE;
     LogTrace("EcalClusters") << "Made supercluster with energy E: " << suCl.energy();
 
   }  //end loop over map
   sort(SCcoll.rbegin(), SCcoll.rend(), isClusterEtLess);
   return SCcoll;
 }
 
 double HybridClusterAlgo::makeDomino(EcalBarrelNavigatorHT& navigator, std::vector<EcalRecHit>& cells) {
   //At the beginning of this function, the navigator starts at the middle of the domino,
   //and that's where EcalBarrelNavigator::home() should send it.
   //Walk one crystal in eta to either side of the initial point.  Sum the three cell energies.
   //If the resultant energy is larger than Ewing, then go an additional cell to either side.
   //Returns:  Total domino energy.  Also, stores the cells used to create domino in the vector.
   cells.clear();
   double Etot = 0;
 
   //Ready?  Get the starting cell.
   DetId center = navigator.pos();
   EcalRecHitCollection::const_iterator center_it = recHits_->find(center);
 
   if (center_it != recHits_->end()) {
     EcalRecHit SeedHit = *center_it;
     if (useddetids.find(center) == useddetids.end() && excludedCrys_.find(center) == excludedCrys_.end()) {
       Etot += SeedHit.energy();
       cells.push_back(SeedHit);
     }
   }
   //One step upwards in Ieta:
   DetId ieta1 = navigator.west();
   EcalRecHitCollection::const_iterator eta1_it = recHits_->find(ieta1);
   if (eta1_it != recHits_->end()) {
     EcalRecHit UpEta = *eta1_it;
     if (useddetids.find(ieta1) == useddetids.end() && excludedCrys_.find(ieta1) == excludedCrys_.end()) {
       Etot += UpEta.energy();
       cells.push_back(UpEta);
     }
   }
 
   //Go back to the middle.
   navigator.home();
 
   //One step downwards in Ieta:
   DetId ieta2 = navigator.east();
   EcalRecHitCollection::const_iterator eta2_it = recHits_->find(ieta2);
   if (eta2_it != recHits_->end()) {
     EcalRecHit DownEta = *eta2_it;
     if (useddetids.find(ieta2) == useddetids.end() && excludedCrys_.find(ieta2) == excludedCrys_.end()) {
       Etot += DownEta.energy();
       cells.push_back(DownEta);
     }
   }
 
   //Now check the energy.  If smaller than Ewing, then we're done.  If greater than Ewing, we have to
   //add two additional cells, the 'wings'
   if (Etot < Ewing) {
     navigator.home();  //Needed even here!!
     return Etot;       //Done!  Not adding 'wings'.
   }
 
   //Add the extra 'wing' cells.  Remember, we haven't sent the navigator home,
   //we're still on the DownEta cell.
   if (ieta2 != DetId(0)) {
     DetId ieta3 = navigator.east();  //Take another step downward.
     EcalRecHitCollection::const_iterator eta3_it = recHits_->find(ieta3);
     if (eta3_it != recHits_->end()) {
       EcalRecHit DownEta2 = *eta3_it;
       if (useddetids.find(ieta3) == useddetids.end() && excludedCrys_.find(ieta3) == excludedCrys_.end()) {
         Etot += DownEta2.energy();
         cells.push_back(DownEta2);
       }
     }
   }
   //Now send the navigator home.
   navigator.home();
   navigator.west();  //Now you're on eta1_it
   if (ieta1 != DetId(0)) {
     DetId ieta4 = navigator.west();  //Take another step upward.
     EcalRecHitCollection::const_iterator eta4_it = recHits_->find(ieta4);
     if (eta4_it != recHits_->end()) {
       EcalRecHit UpEta2 = *eta4_it;
       if (useddetids.find(ieta4) == useddetids.end() && excludedCrys_.find(ieta4) == excludedCrys_.end()) {
         Etot += UpEta2.energy();
         cells.push_back(UpEta2);
       }
     }
   }
   navigator.home();
   return Etot;
 }
 
 double HybridClusterAlgo::et25(EcalBarrelNavigatorHT& navigator,
                                const EcalRecHitCollection* hits,
                                const CaloSubdetectorGeometry* geometry) {
   DetId thisDet;
   std::vector<std::pair<DetId, float> > dets;
   dets.clear();
   EcalRecHitCollection::const_iterator hit;
   double energySum = 0.0;
 
   for (int dx = -2; dx < 3; ++dx) {
     for (int dy = -2; dy < 3; ++dy) {
       //std::cout << "dx, dy " << dx << ", " << dy << std::endl;
       thisDet = navigator.offsetBy(dx, dy);
       navigator.home();
 
       if (thisDet != DetId(0)) {
         hit = recHits_->find(thisDet);
         if (hit != recHits_->end()) {
           dets.push_back(std::pair<DetId, float>(thisDet, 1.));  // by default hit energy fraction is set to 1
           energySum += hit->energy();
         }
       }
     }
   }
 
   // convert it to ET
   //std::cout << "dets.size(), energySum: " << dets.size() << ", " << energySum << std::endl;
   Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
   double et = energySum / cosh(pos.eta());
   return et;
 }
 
 double HybridClusterAlgo::e2Et(EcalBarrelNavigatorHT& navigator,
                                const EcalRecHitCollection* hits,
                                const CaloSubdetectorGeometry* geometry) {
   DetId thisDet;
   std::vector<std::pair<DetId, float> > dets;
   dets.clear();
   EcalRecHitCollection::const_iterator hit;
 
   for (int dx = -2; dx < 3; ++dx) {
     for (int dy = -2; dy < 3; ++dy) {
       thisDet = navigator.offsetBy(dx, dy);
       navigator.home();
 
       if (thisDet != DetId(0)) {
         hit = recHits_->find(thisDet);
         if (hit != recHits_->end()) {
           dets.push_back(std::pair<DetId, float>(thisDet, 1.));  // by default hit energy fraction is set to 1
         }
       }
     }
   }
 
   // compute coefficient to turn E into Et
   Point pos = posCalculator_.Calculate_Location(dets, hits, geometry);
   return 1 / cosh(pos.eta());
 }

This page was automatically generated by the 2.2.1 LXR engine. The LXR team