HcalHits/src/SimG4HcalHitJetFinder.cc

0001 ///////////////////////////////////////////////////////////////////////////////
0002 // File: SimG4HcalHitJetFinder.cc
0003 // Description: Jet finder class for SimG4HcalValidation
0004 ///////////////////////////////////////////////////////////////////////////////
0005 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Validation/HcalHits/interface/SimG4HcalHitJetFinder.h"
0008
0009 #include <algorithm>
0010 #include <cmath>
0011 #include <iostream>
0012
0013 SimG4HcalHitJetFinder::SimG4HcalHitJetFinder(double cone) : jetcone(cone) {}
0014
0015 SimG4HcalHitJetFinder::~SimG4HcalHitJetFinder() { edm::LogVerbatim("ValidHcal") << "SimG4HcalHitJetFinder:: Deleting"; }
0016
0017 void SimG4HcalHitJetFinder::setCone(double cone) { jetcone = cone; }
0018
0019 void SimG4HcalHitJetFinder::setInput(std::vector<CaloHit> *hhit) { input = *hhit; }
0020
0021 std::vector<SimG4HcalHitCluster> *SimG4HcalHitJetFinder::getClusters(bool hcal_only) {
0022   clusvector.erase(clusvector.begin(), clusvector.end());
0023   if (input.empty()) {
0024     return &clusvector;
0025   }
0026
0027   std::vector<CaloHit>::iterator itr;
0028   for (itr = input.begin(); itr != input.end(); itr++) {
0029     edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder::getClusters_1 - input :  e " << itr->e() << "  eta "
0030                                   << itr->eta() << "  phi " << itr->phi() << "  subdet " << itr->det();
0031   }
0032
0033   sort(input.begin(), input.end());  // sort input in descending order
0034
0035   for (itr = input.begin(); itr != input.end(); itr++) {
0036     edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder::getClusters_2 - input :  e " << itr->e() << "  eta "
0037                                   << itr->eta() << "  phi " << itr->phi() << "  subdet " << itr->det();
0038   }
0039
0040   std::vector<SimG4HcalHitCluster> temp;  // dummy container for clusters
0041
0042   //  first input hit -> first cluster
0043
0044   CaloHit hit;
0045   SimG4HcalHitCluster cluster;
0046
0047   std::vector<CaloHit>::iterator itr_hits;
0048
0049   int j, first_seed = 0;
0050   for (j = 0, itr_hits = input.begin(); itr_hits != input.end(); j++, itr_hits++) {
0051     int h_type = itr_hits->det();  // if desired HCAL hits (only) clusterfinding
0052     if (((h_type == static_cast<int>(HcalBarrel) || h_type == static_cast<int>(HcalEndcap) ||
0053           h_type == static_cast<int>(HcalForward)) &&
0054          hcal_only) ||
0055         (!hcal_only)) {
0056       cluster += input[j];
0057       edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder:: First seed hit ..................\n" << (*itr_hits);
0058       first_seed = j;
0059       break;
0060     }
0061   }
0062
0063   temp.push_back(cluster);
0064
0065   std::vector<SimG4HcalHitCluster>::iterator itr_clus;
0066
0067   for (j = 0, itr_hits = input.begin(); itr_hits != input.end(); j++, itr_hits++) {
0068     int h_type = itr_hits->det();  // if desired HCAL hits (only) clusterfinding
0069     if ((((h_type == static_cast<int>(HcalBarrel) || h_type == static_cast<int>(HcalEndcap) ||
0070            h_type == static_cast<int>(HcalForward)) &&
0071           hcal_only) ||
0072          (!hcal_only)) &&
0073         (j != first_seed)) {
0074       edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder:: ........... Consider hit ..................\n"
0075                                     << (*itr_hits);
0076
0077       int incl = 0;  // if the hit is included in one of clusters
0078
0079       int iclus;
0080       for (itr_clus = temp.begin(), iclus = 0; itr_clus != temp.end(); itr_clus++, iclus++) {
0081         edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder::=======> Cluster " << iclus << "\n" << (*itr_clus);
0082
0083         double d = rDist(&(*itr_clus), &(*itr_hits));
0084         if (d < jetcone) {
0085           edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder:: -> associated ... ";
0086           temp[iclus] += *itr_hits;
0087           incl = 1;
0088           break;
0089         }
0090       }
0091
0092       // to here jumps "break"
0093       if (incl == 0) {
0094         SimG4HcalHitCluster cl;
0095         cl += *itr_hits;
0096         temp.push_back(cl);
0097         edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder:: ************ NEW CLUSTER !\n" << cl;
0098       }
0099     }
0100   }
0101
0102   clusvector = temp;
0103   return &clusvector;
0104 }
0105
0106 double SimG4HcalHitJetFinder::rDist(const SimG4HcalHitCluster *cluster, const CaloHit *hit) const {
0107   double etac = cluster->eta();
0108   double phic = cluster->phi();
0109
0110   double etah = hit->eta();
0111   double phih = hit->phi();
0112
0113   return rDist(etac, phic, etah, phih);
0114 }
0115
0116 double SimG4HcalHitJetFinder::rDist(const double etac, const double phic, const double etah, const double phih) const {
0117   double delta_eta = etac - etah;
0118   double delta_phi = phic - phih;
0119
0120   if (phic < phih)
0121     delta_phi = phih - phic;
0122   if (delta_phi > M_PI)
0123     delta_phi = 2 * M_PI - delta_phi;
0124
0125   double tmp = sqrt(delta_eta * delta_eta + delta_phi * delta_phi);
0126
0127   edm::LogVerbatim("ValidHcal") << "HcalHitJetFinder::rDist:\n Clus. eta, phi = " << etac << " " << phic
0128                                 << "\n hit   eta, phi = " << etah << " " << phih << " rDist = " << tmp;
0129
0130   return tmp;
0131 }