EgammaClusterAlgos/src/CosmicClusterAlgo.cc

0001
0002 #include "RecoEcal/EgammaClusterAlgos/interface/CosmicClusterAlgo.h"
0003
0004 #include <vector>  //TEMP JHAUPT 4-27
0005
0006 // Geometry
0007 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0008 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0009 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0010 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0011 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0012 #include "RecoEcal/EgammaCoreTools/interface/ClusterEtLess.h"
0013
0014 // Return a vector of clusters from a collection of EcalRecHits:
0015 //
0016 std::vector<reco::BasicCluster> CosmicClusterAlgo::makeClusters(const EcalRecHitCollection *hits,
0017                                                                 const EcalUncalibratedRecHitCollection *uncalibhits,
0018                                                                 const CaloSubdetectorGeometry *geometry_p,
0019                                                                 const CaloSubdetectorTopology *topology_p,
0020                                                                 const CaloSubdetectorGeometry *geometryES_p,
0021                                                                 EcalPart ecalPart,
0022                                                                 bool regional,
0023                                                                 const std::vector<RectangularEtaPhiRegion> &regions) {
0024   seeds.clear();
0025   used_s.clear();
0026   canSeed_s.clear();
0027   clusters_v.clear();
0028
0029   recHits_ = hits;
0030   uncalibRecHits_ = uncalibhits;
0031
0032   inEB = true;
0033   double threshold = 0;
0034   std::string ecalPart_string;
0035   if (ecalPart == endcap) {
0036     threshold = ecalEndcapSeedThreshold;
0037     ecalPart_string = "EndCap";
0038     inEB = false;
0039   }
0040   if (ecalPart == barrel) {
0041     threshold = ecalBarrelSeedThreshold;
0042     ecalPart_string = "Barrel";
0043     inEB = true;
0044   }
0045
0046   if (verbosity < pINFO) {
0047     std::cout << "-------------------------------------------------------------" << std::endl;
0048     std::cout << "Island algorithm invoked for ECAL" << ecalPart_string << std::endl;
0049     std::cout << "Looking for seeds, threshold used = " << threshold << " ADC" << std::endl;
0050   }
0051
0052   int nregions = 0;
0053   if (regional)
0054     nregions = regions.size();
0055
0056   if (!regional || nregions) {
0057     EcalRecHitCollection::const_iterator it;
0058     for (it = recHits_->begin(); it != recHits_->end(); it++) {
0059       if (!uncalibRecHits_) {
0060         if (verbosity < pINFO) {
0061           std::cout << "-------------------------------------------------------------" << std::endl;
0062           std::cout << "No Uncalibrated RecHits no Uncalibrated rec hit collection available" << std::endl;
0063         }
0064         continue;
0065       }
0066
0067       // if rechit affected by features other than these, do not allow if seeding
0068       // features are hard coded, for now; add severity?
0069       uint32_t rhFlag = (*it).recoFlag();
0070       if (!(rhFlag == EcalRecHit::kGood || rhFlag == EcalRecHit::kOutOfTime || rhFlag == EcalRecHit::kPoorCalib))
0071         continue;
0072
0073       EcalUncalibratedRecHitCollection::const_iterator itt = uncalibRecHits_->find(it->id());
0074
0075       if (itt == uncalibRecHits_->end()) {
0076         if (verbosity < pINFO) {
0077           std::cout << "-------------------------------------------------------------" << std::endl;
0078           std::cout << "No Uncalibrated RecHit associated with the RecHit Probably no Uncalibrated rec hit collection "
0079                        "available"
0080                     << std::endl;
0081         }
0082         continue;
0083       }
0084
0085       EcalUncalibratedRecHit uhit_p = *itt;
0086
0087       // looking for cluster seeds
0088       if (uhit_p.amplitude() < (inEB ? ecalBarrelSeedThreshold : ecalEndcapSeedThreshold))
0089         continue;  //
0090
0091       auto thisCell = geometry_p->getGeometry(it->id());
0092
0093       // Require that RecHit is within clustering region in case
0094       // of regional reconstruction
0095       bool withinRegion = false;
0096       if (regional) {
0097         std::vector<RectangularEtaPhiRegion>::const_iterator region;
0098         for (region = regions.begin(); region != regions.end(); region++) {
0099           if (region->inRegion(thisCell->etaPos(), thisCell->phiPos())) {
0100             withinRegion = true;
0101             break;
0102           }
0103         }
0104       }
0105
0106       if (!regional || withinRegion) {
0107         seeds.push_back(
0108             *it);  // JHaupt 4-27-2008 Et -> energy, most likely not needed as there is already a threshold requirement.
0109       }
0110     }
0111   }
0112
0113   sort(seeds.begin(), seeds.end(), [](auto const &x, auto const &y) { return x.energy() > y.energy(); });
0114
0115   if (verbosity < pINFO) {
0116     std::cout << "JH Total number of seeds found in event = " << seeds.size() << std::endl;
0117     for (EcalRecHitCollection::const_iterator ji = seeds.begin(); ji != seeds.end(); ++ji) {
0118       //std::cout << "JH Seed Energy " << ji->energy() << " hashed " << ((EBDetId)ji->id()).hashedIndex()  << std::endl;
0119     }
0120   }
0121
0122   mainSearch(geometry_p, topology_p, geometryES_p, ecalPart);
0123   std::sort(clusters_v.rbegin(), clusters_v.rend(), isClusterEtLess);
0124
0125   if (verbosity < pINFO) {
0126     std::cout << "---------- end of main search. clusters have been sorted ----" << std::endl;
0127   }
0128
0129   return clusters_v;
0130 }
0131
0132 // Search for clusters
0133 //
0134 void CosmicClusterAlgo::mainSearch(const CaloSubdetectorGeometry *geometry_p,
0135                                    const CaloSubdetectorTopology *topology_p,
0136                                    const CaloSubdetectorGeometry *geometryES_p,
0137                                    EcalPart ecalPart) {
0138   if (verbosity < pINFO) {
0139     std::cout << "Building clusters............" << std::endl;
0140   }
0141
0142   // Loop over seeds:
0143   std::vector<EcalRecHit>::iterator it;
0144   for (it = seeds.begin(); it != seeds.end(); it++) {
0145     // check if this crystal is able to seed
0146     // (event though it is already used)
0147     bool usedButCanSeed = false;
0148     if (canSeed_s.find(it->id()) != canSeed_s.end())
0149       usedButCanSeed = true;
0150
0151     // make sure the current seed does not belong to a cluster already.
0152     if ((used_s.find(it->id()) != used_s.end()) && (usedButCanSeed == false)) {
0153       if (it == seeds.begin()) {
0154         if (verbosity < pINFO) {
0155           std::cout << "##############################################################" << std::endl;
0156           std::cout << "DEBUG ALERT: Highest energy seed already belongs to a cluster!" << std::endl;
0157           std::cout << "##############################################################" << std::endl;
0158         }
0159       }
0160
0161       // seed crystal is used or is used and cannot seed a cluster
0162       // so continue to the next seed crystal...
0163       continue;
0164     }
0165
0166     // clear the vector of hits in current clusterconst EcalRecHitCollection* hits,
0167     current_v9.clear();
0168     current_v25.clear();
0169     current_v25Sup.clear();
0170
0171     // Create a navigator at the seed
0172     CaloNavigator<DetId> navigator(it->id(), topology_p);
0173     DetId seedId = navigator.pos();
0174     navigator.setHome(seedId);
0175
0176     // Is the seed a local maximum?
0177     bool localMaxima = checkMaxima(navigator);
0178
0179     if (localMaxima) {
0180       // build the 5x5 taking care over which crystals //JHaupt 4-27-08 3x3 is a good option...
0181       // can seed new clusters and which can't
0182       prepareCluster(navigator, geometry_p);
0183     }
0184
0185     // If some crystals in the current vector then
0186     // make them into a cluster
0187     if (!current_v25.empty()) {
0188       makeCluster(geometry_p, geometryES_p, seedId);
0189     }
0190
0191   }  // End loop on seed crystals
0192 }
0193
0194 void CosmicClusterAlgo::makeCluster(const CaloSubdetectorGeometry *geometry,
0195                                     const CaloSubdetectorGeometry *geometryES,
0196                                     DetId seedId) {
0197   double energy = 0;
0198   double energySecond = 0.;  //JHaupt 4-27-08 Added for the second crystal stream
0199   double energyMax = 0.;     //JHaupt 4-27-08 Added for the max crystal stream
0200   DetId detFir;
0201   DetId detSec;
0202   //bool goodCluster = false; //JHaupt 4-27-08 Added so that some can be earased.. used another day Might not be needed as seeds are energy ordered...
0203
0204   std::vector<DetId>::iterator it;
0205   for (it = current_v9.begin(); it != current_v9.end(); it++) {
0206     // Martina - check recoFlag for crystals sorrounding the good one
0207     EcalRecHitCollection::const_iterator itt = recHits_->find(*it);
0208     // double-check that iterator can be dereferenced
0209     if (itt == recHits_->end())
0210       continue;
0211     // if rechit affected by features other than these, do not allow 2 crystal seeding either
0212     // features are hard coded, for now; add severity?
0213     uint32_t rhFlag = (*itt).recoFlag();
0214     if (!(rhFlag == EcalRecHit::kGood || rhFlag == EcalRecHit::kOutOfTime || rhFlag == EcalRecHit::kPoorCalib))
0215       continue;
0216
0217     ///
0218
0219     EcalUncalibratedRecHitCollection::const_iterator ittu = uncalibRecHits_->find(*it);
0220     EcalUncalibratedRecHit uhit_p = *ittu;
0221
0222     if (uhit_p.amplitude() > energySecond) {
0223       energySecond = uhit_p.amplitude();
0224       detSec = uhit_p.id();
0225     }
0226     if (energySecond > energyMax) {
0227       std::swap(energySecond, energyMax);
0228       std::swap(detFir, detSec);
0229     }
0230   }
0231
0232   //if ((detFir.det() == DetId::Ecal) && (detFir.subdetId() == EcalEndcap)) inEB = false;
0233   //We ignore the possiblity that the cluster may span into the EE
0234
0235   if ((energyMax < (inEB ? ecalBarrelSingleThreshold : ecalEndcapSingleThreshold)) &&
0236       (energySecond < (inEB ? ecalBarrelSecondThreshold : ecalEndcapSecondThreshold)))
0237     return;
0238
0239   for (it = current_v25.begin(); it != current_v25.end(); it++) {
0240     EcalRecHitCollection::const_iterator itt = recHits_->find(*it);
0241     EcalRecHit hit_p = *itt;
0242
0243     EcalUncalibratedRecHitCollection::const_iterator ittu = uncalibRecHits_->find(*it);
0244     EcalUncalibratedRecHit uhit_p = *ittu;
0245
0246     if (uhit_p.amplitude() > (inEB ? ecalBarrelSupThreshold : ecalEndcapSupThreshold)) {
0247       // energy fraction = 1
0248       current_v25Sup.push_back(std::pair<DetId, float>(hit_p.id(), 1.));
0249       energy += hit_p.energy();  //Keep the fully corrected energy
0250     }
0251   }
0252
0253   Point position;
0254   position = posCalculator_.Calculate_Location(current_v25Sup, recHits_, geometry, geometryES);
0255
0256   //don't write empty clusters
0257   if (energy == 0 && position == Point(0, 0, 0))
0258     return;
0259
0260   if (verbosity < pINFO) {
0261     std::cout << "JH******** NEW CLUSTER ********" << std::endl;
0262     std::cout << "JHNo. of crystals = " << current_v25Sup.size() << std::endl;
0263     std::cout << "JH     Energy     = " << energy << std::endl;
0264     std::cout << "JH     Phi        = " << position.phi() << std::endl;
0265     std::cout << "JH     Eta        = " << position.eta() << std::endl;
0266     std::cout << "JH*****************************" << std::endl;
0267     // specialize this
0268     //       std::cout << "JH****Emax****  "<<energyMax << " ieta " <<detFir.ieta() <<" iphi "<<detFir.ieta()  << std::endl;
0269     //       std::cout << "JH****Esec****  "<<energySecond << " ieta " <<detSec.ieta() <<" iphi "<<detSec.ieta() << std::endl;
0270   }
0271   clusters_v.push_back(
0272       reco::BasicCluster(energy, position, reco::CaloID(), current_v25Sup, reco::CaloCluster::multi5x5, seedId));
0273 }
0274
0275 bool CosmicClusterAlgo::checkMaxima(CaloNavigator<DetId> &navigator) {
0276   bool maxima = true;
0277   EcalRecHitCollection::const_iterator thisHit;
0278   EcalRecHitCollection::const_iterator seedHit = recHits_->find(navigator.pos());
0279   double thisEnergy = 0.;
0280   double seedEnergy = seedHit->energy();
0281
0282   std::vector<DetId> swissCrossVec;
0283   swissCrossVec.clear();
0284
0285   swissCrossVec.push_back(navigator.west());
0286   navigator.home();
0287   swissCrossVec.push_back(navigator.east());
0288   navigator.home();
0289   swissCrossVec.push_back(navigator.north());
0290   navigator.home();
0291   swissCrossVec.push_back(navigator.south());
0292   navigator.home();
0293
0294   for (unsigned int i = 0; i < swissCrossVec.size(); ++i) {
0295     thisHit = recHits_->find(swissCrossVec[i]);
0296     if ((swissCrossVec[i] == DetId(0)) || thisHit == recHits_->end())
0297       thisEnergy = 0.0;
0298     else
0299       thisEnergy = thisHit->energy();
0300     if (thisEnergy > seedEnergy) {
0301       maxima = false;
0302       break;
0303     }
0304   }
0305
0306   return maxima;
0307 }
0308
0309 void CosmicClusterAlgo::prepareCluster(CaloNavigator<DetId> &navigator, const CaloSubdetectorGeometry *geometry) {
0310   DetId thisDet;
0311   std::set<DetId>::iterator setItr;
0312
0313   // now add the 5x5 taking care to mark the edges
0314   // as able to seed and where overlapping in the central
0315   // region with crystals that were previously able to seed
0316   // change their status so they are not able to seed
0317   //std::cout << std::endl;
0318   for (int dx = -2; dx < 3; ++dx)  //for (int dx = -2; dx < 3; ++dx)
0319   {
0320     for (int dy = -2; dy < 3; ++dy)  //for (int dy = -2; dy < 3; ++ dy)
0321     {
0322       // navigate in free steps forming
0323       // a full 5x5
0324       thisDet = navigator.offsetBy(dx, dy);
0325       navigator.home();
0326
0327       // add the current crystal
0328       //std::cout << "adding " << dx << ", " << dy << std::endl;
0329
0330       // now consider if we are in an edge (outer 16)
0331       // or central (inner 9) region
0332       if ((abs(dx) > 1) || (abs(dy) > 1)) {
0333         // this is an "edge" so should be allowed to seed
0334         // provided it is not already used
0335         //std::cout << "   setting can seed" << std::endl;
0336         addCrystal(thisDet, false);  //These are in the V25
0337         canSeed_s.insert(thisDet);
0338       }  // end if "edge"
0339       else {
0340         // or else we are in the central 3x3
0341         // and must remove any of these crystals from the canSeed set
0342         setItr = canSeed_s.find(thisDet);
0343         addCrystal(thisDet, true);  //These are in the V9
0344         if (setItr != canSeed_s.end()) {
0345           //std::cout << "   unsetting can seed" << std::endl;
0346           canSeed_s.erase(setItr);
0347         }
0348       }  // end if "centre"
0349
0350     }  // end loop on dy
0351
0352   }  // end loop on dx
0353
0354   //std::cout << "*** " << std::endl;
0355   //std::cout << " current_v contains " << current_v.size() << std::endl;
0356   //std::cout << "*** " << std::endl;
0357 }
0358
0359 void CosmicClusterAlgo::addCrystal(const DetId &det, const bool in9) {
0360   EcalRecHitCollection::const_iterator thisIt = recHits_->find(det);
0361   if ((thisIt != recHits_->end()) && (thisIt->id() != DetId(0))) {
0362     if ((used_s.find(thisIt->id()) == used_s.end())) {
0363       EcalUncalibratedRecHitCollection::const_iterator thisItu = uncalibRecHits_->find(det);
0364       used_s.insert(det);
0365       if ((thisIt->energy() >= -1.) && !(thisItu->chi2() < -1.)) {
0366         if (in9)
0367           current_v9.push_back(det);
0368         current_v25.push_back(det);
0369       }
0370     }
0371   }
0372 }