EgammaClusterAlgos/src/PreshowerClusterAlgo.cc

0001 #include <vector>
0002 #include <map>
0003
0004 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0005 #include "RecoEcal/EgammaClusterAlgos/interface/PreshowerClusterAlgo.h"
0006 #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
0007 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include "TH1.h"
0012
0013 reco::PreshowerCluster PreshowerClusterAlgo::makeOneCluster(ESDetId strip,
0014                                                             HitsID* used_strips,
0015                                                             RecHitsMap* the_rechitsMap_p,
0016                                                             const CaloSubdetectorGeometry* geometry_p,
0017                                                             const CaloSubdetectorTopology* topology_p) {
0018   road_2d.clear();
0019
0020   rechits_map = the_rechitsMap_p;
0021
0022   used_s = used_strips;
0023
0024   int plane = strip.plane();
0025
0026   LogTrace("PreShowerClusterAlgo") << "Preshower Seeded Algorithm - looking for clusters";
0027   LogTrace("PreShowerClusterAlgo") << "Preshower is intersected at strip" << strip.strip() << ",at plane" << plane;
0028
0029   // create null-cluster
0030   std::vector<std::pair<DetId, float> > dummy;
0031   Point posi(0, 0, 0);
0032   LogTrace("PreShowerClusterAlgo") << " Creating a null-cluster";
0033
0034   reco::PreshowerCluster nullcluster = reco::PreshowerCluster(0., posi, dummy, plane);
0035
0036   if (strip == ESDetId(0))
0037     return nullcluster;  //works in case of no intersected strip found (e.g. in the Barrel)
0038
0039   // Collection of cluster strips
0040   EcalRecHitCollection clusterRecHits;
0041   // Map of strips for position calculation
0042   RecHitsMap recHits_pos;
0043
0044   //Make a navigator, and set it to the strip cell.
0045   EcalPreshowerNavigator navigator(strip, topology_p);
0046   navigator.setHome(strip);
0047   //search for neighbours in the central road
0048   findRoad(strip, navigator, plane);
0049   LogTrace("PreShowerClusterAlgo") << "Total number of strips in the central road:" << road_2d.size();
0050
0051   if (plane == 1) {
0052     ESDetId strip_north = navigator.north();
0053     findRoad(strip_north, navigator, plane);
0054     navigator.home();
0055     ESDetId strip_south = navigator.south();
0056     findRoad(strip_south, navigator, plane);
0057     navigator.home();
0058   }
0059   if (plane == 2) {
0060     ESDetId strip_east = navigator.east();
0061     findRoad(strip_east, navigator, plane);
0062     navigator.home();
0063     ESDetId strip_west = navigator.west();
0064     findRoad(strip_west, navigator, plane);
0065     navigator.home();
0066   }
0067   LogTrace("PreShowerClusterAlgo") << "Total number of strips in all three roads:" << road_2d.size();
0068
0069   // Start clustering from strip with max Energy in the road
0070   float E_max = 0.;
0071   bool found = false;
0072   RecHitsMap::iterator max_it;
0073   // Loop over strips:
0074   std::vector<ESDetId>::iterator itID;
0075   for (itID = road_2d.begin(); itID != road_2d.end(); itID++) {
0076     LogTrace("PreShowerClusterAlgo") << "ID =" << *itID;
0077
0078     RecHitsMap::iterator strip_it = rechits_map->find(*itID);
0079     //if ( strip_it->second.energy() < 0 ) std::cout << "           ##### E = " << strip_it->second.energy() << std::endl;
0080     if (strip_it == rechits_map->end() || !goodStrip(strip_it))
0081       continue;
0082     LogTrace("PreShowerClusterAlgo") << " strip is " << ESDetId(strip_it->first) << "E =" << strip_it->second.energy();
0083
0084     float E = strip_it->second.energy();
0085     if (E > E_max) {
0086       E_max = E;
0087       found = true;
0088       max_it = strip_it;
0089     }
0090   }
0091   // no hottest strip found ==> null cluster will be returned
0092   if (!found)
0093     return nullcluster;
0094
0095   // First, save the hottest strip
0096   clusterRecHits.push_back(max_it->second);
0097   recHits_pos.insert(std::make_pair(max_it->first, max_it->second));
0098   used_s->insert(max_it->first);
0099   LogTrace("PreShowerClusterAlgo") << "Central hottest strip" << ESDetId(max_it->first)
0100                                    << "is saved with energy E =" << E_max;
0101
0102   // Find positions of adjacent strips:
0103   ESDetId next, strip_1, strip_2;
0104   navigator.setHome(max_it->first);
0105   ESDetId startES = max_it->first;
0106
0107   if (plane == 1) {
0108     // Save two neighbouring strips to the east
0109     int nadjacents_east = 0;
0110     while ((next = navigator.east()) != ESDetId(0) && next != startES && nadjacents_east < 2) {
0111       ++nadjacents_east;
0112       LogTrace("PreShowerClusterAlgo") << "Adjacent east #" << nadjacents_east << ":" << next;
0113
0114       RecHitsMap::iterator strip_it = rechits_map->find(next);
0115
0116       if (strip_it == rechits_map->end() || !goodStrip(strip_it))
0117         continue;
0118       // Save strip for clustering if it exists, not already in use, and satisfies an energy threshold
0119       clusterRecHits.push_back(strip_it->second);
0120       // save strip for position calculation
0121       if (nadjacents_east == 1)
0122         strip_1 = next;
0123       used_s->insert(strip_it->first);
0124       LogTrace("PreShowerClusterAlgo") << "East adjacent strip #" << nadjacents_east
0125                                        << "is saved with energy E =" << strip_it->second.energy();
0126     }
0127     // Save two neighbouring strips to the west
0128     navigator.home();
0129     int nadjacents_west = 0;
0130     while ((next = navigator.west()) != ESDetId(0) && next != startES && nadjacents_west < 2) {
0131       ++nadjacents_west;
0132       LogTrace("PreShowerClusterAlgo") << "Adjacent west #" << nadjacents_west << ":" << next;
0133
0134       RecHitsMap::iterator strip_it = rechits_map->find(next);
0135       if (strip_it == rechits_map->end() || !goodStrip(strip_it))
0136         continue;
0137       clusterRecHits.push_back(strip_it->second);
0138       if (nadjacents_west == 1)
0139         strip_2 = next;
0140       used_s->insert(strip_it->first);
0141       LogTrace("PreShowerClusterAlgo") << "West adjacent strip #" << nadjacents_west
0142                                        << "is saved with energy E =" << strip_it->second.energy();
0143     }
0144   } else if (plane == 2) {
0145     // Save two neighbouring strips to the north
0146     int nadjacents_north = 0;
0147     while ((next = navigator.north()) != ESDetId(0) && next != startES && nadjacents_north < 2) {
0148       ++nadjacents_north;
0149       LogTrace("PreShowerClusterAlgo") << "Adjacent north #" << nadjacents_north << ":" << next;
0150
0151       RecHitsMap::iterator strip_it = rechits_map->find(next);
0152       if (strip_it == rechits_map->end() || !goodStrip(strip_it))
0153         continue;
0154       clusterRecHits.push_back(strip_it->second);
0155       if (nadjacents_north == 1)
0156         strip_1 = next;
0157       used_s->insert(strip_it->first);
0158       LogTrace("PreShowerClusterAlgo") << "North adjacent strip #" << nadjacents_north
0159                                        << "is saved with energy E =" << strip_it->second.energy();
0160     }
0161     // Save two neighbouring strips to the south
0162     navigator.home();
0163     int nadjacents_south = 0;
0164     while ((next = navigator.south()) != ESDetId(0) && next != startES && nadjacents_south < 2) {
0165       ++nadjacents_south;
0166       LogTrace("PreShowerClusterAlgo") << "Adjacent south #" << nadjacents_south << ":" << next;
0167
0168       RecHitsMap::iterator strip_it = rechits_map->find(next);
0169       if (strip_it == rechits_map->end() || !goodStrip(strip_it))
0170         continue;
0171       clusterRecHits.push_back(strip_it->second);
0172       if (nadjacents_south == 1)
0173         strip_2 = next;
0174       used_s->insert(strip_it->first);
0175       LogTrace("PreShowerClusterAlgo") << "South adjacent strip #" << nadjacents_south
0176                                        << "is saved with energy E =" << strip_it->second.energy();
0177     }
0178   } else {
0179     LogTrace("PreShowerClusterAlgo") << " Wrong plane number" << plane << ", null cluster will be returned! "
0180                                      << std::endl;
0181     return nullcluster;
0182   }  // end of if
0183
0184   LogTrace("PreShowerClusterAlgo") << "Total size of clusterRecHits is" << clusterRecHits.size();
0185   LogTrace("PreShowerClusterAlgo") << "Two adjacent strips for position calculation are:" << strip_1 << "and"
0186                                    << strip_2;
0187
0188   // strips for position calculation
0189   RecHitsMap::iterator strip_it1, strip_it2;
0190   if (strip_1 != ESDetId(0)) {
0191     strip_it1 = rechits_map->find(strip_1);
0192     recHits_pos.insert(std::make_pair(strip_it1->first, strip_it1->second));
0193   }
0194   if (strip_2 != ESDetId(0)) {
0195     strip_it2 = rechits_map->find(strip_2);
0196     recHits_pos.insert(std::make_pair(strip_it2->first, strip_it2->second));
0197   }
0198
0199   RecHitsMap::iterator cp;
0200   double energy_pos = 0;
0201   double x_pos = 0;
0202   double y_pos = 0;
0203   double z_pos = 0;
0204   for (cp = recHits_pos.begin(); cp != recHits_pos.end(); cp++) {
0205     double E = cp->second.energy();
0206     energy_pos += E;
0207     auto thisCell = geometry_p->getGeometry(cp->first);
0208     const GlobalPoint& position = thisCell->getPosition();
0209     x_pos += E * position.x();
0210     y_pos += E * position.y();
0211     z_pos += E * position.z();
0212   }
0213   if (energy_pos > 0.) {
0214     x_pos /= energy_pos;
0215     y_pos /= energy_pos;
0216     z_pos /= energy_pos;
0217   }
0218   Point pos(x_pos, y_pos, z_pos);
0219   LogTrace("PreShowerClusterAlgo") << "ES Cluster position ="
0220                                    << "(" << x_pos << "," << y_pos << "," << z_pos << ")";
0221
0222   EcalRecHitCollection::iterator it;
0223   double Eclust = 0;
0224
0225   std::vector<std::pair<DetId, float> > usedHits;
0226   for (it = clusterRecHits.begin(); it != clusterRecHits.end(); it++) {
0227     Eclust += it->energy();
0228     usedHits.push_back(std::pair<DetId, float>(it->id(), 1.));
0229   }
0230
0231   // ES cluster is created from vector clusterRecHits
0232   reco::PreshowerCluster cluster = reco::PreshowerCluster(Eclust, pos, usedHits, plane);
0233   LogTrace("PreShowerClusterAlgo") << " ES Cluster is created with:";
0234   LogTrace("PreShowerClusterAlgo") << " energy =" << cluster.energy();
0235   LogTrace("PreShowerClusterAlgo") << " (eta,phi) ="
0236                                    << "(" << cluster.eta() << "," << cluster.phi() << ")";
0237   LogTrace("PreShowerClusterAlgo") << " nhits =" << cluster.nhits();
0238   LogTrace("PreShowerClusterAlgo") << " radius =" << cluster.position().r();
0239   LogTrace("PreShowerClusterAlgo") << " (x,y,z) ="
0240                                    << "(" << cluster.x() << ", " << cluster.y() << "," << cluster.z() << ")";
0241
0242   // return the cluster if its energy is greater a threshold
0243   if (cluster.energy() > preshClusterEnergyCut_)
0244     return cluster;
0245   else
0246     return nullcluster;
0247 }
0248
0249 // returns true if the candidate strip fulfills the requirements to be added to the cluster:
0250 bool PreshowerClusterAlgo::goodStrip(RecHitsMap::iterator candidate_it) {
0251   if (used_s->find(candidate_it->first) != used_s->end())
0252     LogTrace("PreShowerClusterAlgo") << " This strip is in use";
0253   if (candidate_it == rechits_map->end())
0254     LogTrace("PreShowerClusterAlgo") << " No such a strip in rechits_map";
0255   if (candidate_it->second.energy() <= preshStripEnergyCut_)
0256     LogTrace("PreShowerClusterAlgo") << "Strip energy" << candidate_it->second.energy() << "is below threshold";
0257
0258   // crystal should not be included...
0259   if ((used_s->find(candidate_it->first) != used_s->end()) ||   // ...if it already belongs to a cluster
0260       (candidate_it == rechits_map->end()) ||                   // ...if it corresponds to a hit
0261       (candidate_it->second.energy() <= preshStripEnergyCut_))  // ...if it has a negative or zero energy
0262   {
0263     return false;
0264   }
0265   return true;
0266 }
0267
0268 // find strips in the road of size +/- preshSeededNstr_ from the central strip
0269 void PreshowerClusterAlgo::findRoad(ESDetId strip, EcalPreshowerNavigator theESNav, int plane) {
0270   if (strip == ESDetId(0))
0271     return;
0272
0273   ESDetId next;
0274   theESNav.setHome(strip);
0275   // First, add a central strip to the road
0276   road_2d.push_back(strip);
0277   LogTrace("PreShowerClusterAlgo") << "findRoad starts from strip" << strip;
0278
0279   if (plane == 1) {
0280     // east road
0281     int n_east = 0;
0282     LogTrace("PreShowerClusterAlgo") << " Go to the East ";
0283
0284     while (((next = theESNav.east()) != ESDetId(0) && next != strip)) {
0285       LogTrace("PreShowerClusterAlgo") << "East:" << n_east << "current strip is" << next;
0286
0287       road_2d.push_back(next);
0288       ++n_east;
0289       if (n_east == preshSeededNstr_)
0290         break;
0291     }
0292     // west road
0293     int n_west = 0;
0294     LogTrace("PreShowerClusterAlgo") << "Go to the West";
0295
0296     theESNav.home();
0297     while (((next = theESNav.west()) != ESDetId(0) && next != strip)) {
0298       LogTrace("PreShowerClusterAlgo") << "West: " << n_west << "current strip is" << next;
0299
0300       road_2d.push_back(next);
0301       ++n_west;
0302       if (n_west == preshSeededNstr_)
0303         break;
0304     }
0305     LogTrace("PreShowerClusterAlgo") << "Total number of strips found in the road at 1-st plane is" << n_east + n_west;
0306
0307   } else if (plane == 2) {
0308     // north road
0309     int n_north = 0;
0310     LogTrace("PreShowerClusterAlgo") << "Go to the North";
0311
0312     while (((next = theESNav.north()) != ESDetId(0) && next != strip)) {
0313       LogTrace("PreShowerClusterAlgo") << "North:" << n_north << "current strip is" << next;
0314
0315       road_2d.push_back(next);
0316       ++n_north;
0317       if (n_north == preshSeededNstr_)
0318         break;
0319     }
0320     // south road
0321     int n_south = 0;
0322     LogTrace("PreShowerClusterAlgo") << "Go to the South";
0323
0324     theESNav.home();
0325     while (((next = theESNav.south()) != ESDetId(0) && next != strip)) {
0326       LogTrace("PreShowerClusterAlgo") << "South:" << n_south << "current strip is" << next;
0327
0328       road_2d.push_back(next);
0329       ++n_south;
0330       if (n_south == preshSeededNstr_)
0331         break;
0332     }
0333
0334     LogTrace("PreShowerClusterAlgo") << "Total number of strips found in the road at 2-nd plane is"
0335                                      << n_south + n_north;
0336
0337   } else {
0338     LogTrace("PreShowerClusterAlgo") << " Wrong plane number, null cluster will be returned!";
0339
0340   }  // end of if
0341
0342   theESNav.home();
0343 }