CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Validation/HcalHits/src/SimG4HcalHitJetFinder.cc

Go to the documentation of this file.
00001 
00002 // File: SimG4HcalHitJetFinder.cc
00003 // Description: Jet finder class for SimG4HcalValidation
00005 #include "Validation/HcalHits/interface/SimG4HcalHitJetFinder.h"
00006 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
00007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00008 
00009 #include <algorithm>
00010 #include <iostream>
00011 #include <cmath>
00012 
00013 SimG4HcalHitJetFinder::SimG4HcalHitJetFinder(double cone): jetcone(cone) {}
00014 
00015 SimG4HcalHitJetFinder::~SimG4HcalHitJetFinder() {
00016   edm::LogInfo("ValidHcal") << "SimG4HcalHitJetFinder:: Deleting";
00017 }
00018 
00019 void SimG4HcalHitJetFinder::setCone(double cone) { 
00020   jetcone = cone;
00021 }   
00022 
00023 void SimG4HcalHitJetFinder::setInput(std::vector<CaloHit>* hhit) { 
00024   input   = * hhit;
00025 }
00026 
00027 std::vector<SimG4HcalHitCluster>* SimG4HcalHitJetFinder::getClusters(bool hcal_only){
00028   clusvector.erase(clusvector.begin(),clusvector.end()); 
00029   if (input.size() == 0) { 
00030     return &clusvector;
00031   }
00032 
00033   std::vector<CaloHit>::iterator itr;
00034   for (itr = input.begin(); itr != input.end(); itr++) {
00035     LogDebug("ValidHcal") << "HcalHitJetFinder::getClusters_1 - input :  e " 
00036                           << itr->e() << "  eta " << itr->eta() << "  phi " 
00037                           << itr->phi() << "  subdet " << itr->det();
00038   }
00039 
00040   sort(input.begin(),input.end()); // sort input in descending order
00041  
00042   for (itr = input.begin(); itr != input.end(); itr++) {
00043     LogDebug("ValidHcal") << "HcalHitJetFinder::getClusters_2 - input :  e " 
00044                           << itr->e() << "  eta " << itr->eta() << "  phi " 
00045                           << itr->phi() << "  subdet " << itr->det();
00046   }
00047 
00048   std::vector<SimG4HcalHitCluster> temp;   // dummy container for clusters
00049 
00050   //  first input hit -> first cluster
00051 
00052   CaloHit             hit;
00053   SimG4HcalHitCluster cluster;
00054 
00055   std::vector<CaloHit>::iterator itr_hits;
00056 
00057   int j, first_seed = 0;
00058   for (j=0, itr_hits = input.begin(); itr_hits != input.end();
00059        j++, itr_hits++) {
00060     int h_type = itr_hits->det(); //if desired HCAL hits (only) clusterfinding
00061     if (((h_type == static_cast<int>(HcalBarrel) ||
00062           h_type == static_cast<int>(HcalEndcap) ||
00063           h_type == static_cast<int>(HcalForward)) && hcal_only) || 
00064         (!hcal_only)) {
00065       cluster += input[j];
00066       LogDebug("ValidHcal") << "HcalHitJetFinder:: First seed hit "
00067                             << "..................\n" << (*itr_hits);
00068       first_seed = j;
00069       break;
00070     }
00071   }
00072   
00073   temp.push_back(cluster);
00074   
00075   std::vector<SimG4HcalHitCluster>::iterator itr_clus;
00076 
00077   for (j=0, itr_hits = input.begin(); itr_hits != input.end();
00078        j++, itr_hits++) {
00079     int h_type = itr_hits->det(); //if desired HCAL hits (only) clusterfinding
00080     if ((((h_type == static_cast<int>(HcalBarrel) ||
00081            h_type == static_cast<int>(HcalEndcap) ||
00082            h_type == static_cast<int>(HcalForward)) && hcal_only) || 
00083          (!hcal_only)) && (j != first_seed)) {
00084       LogDebug("ValidHcal") << "HcalHitJetFinder:: ........... Consider hit"
00085                             << " ..................\n" << (*itr_hits);
00086       
00087       int incl = 0; // if the hit is included in one of clusters
00088       
00089       int iclus;  
00090       for (itr_clus = temp.begin(), iclus = 0; itr_clus != temp.end();
00091            itr_clus++, iclus++) { 
00092         
00093         LogDebug("ValidHcal") << "HcalHitJetFinder::=======> Cluster " 
00094                               << iclus << "\n" << (*itr_clus);
00095         
00096         double d = rDist(&(*itr_clus), &(*itr_hits));
00097         if (d < jetcone) {
00098           LogDebug("ValidHcal") << "HcalHitJetFinder:: -> associated ... ";
00099           temp[iclus] += *itr_hits;
00100           incl = 1;
00101           break;  
00102         }
00103       }
00104       
00105       // to here jumps "break"
00106       if (incl == 0) {
00107         SimG4HcalHitCluster cl;
00108         cl += *itr_hits;
00109         temp.push_back(cl);
00110         LogDebug("ValidHcal") << "HcalHitJetFinder:: ************ NEW CLUSTER"
00111                               << " !\n" << cl;
00112       }
00113     }
00114   }
00115 
00116   clusvector = temp;
00117   return &clusvector;
00118 }
00119 
00120 double SimG4HcalHitJetFinder::rDist(const SimG4HcalHitCluster* cluster, 
00121                                     const CaloHit* hit) const {
00122 
00123   double etac = cluster->eta();
00124   double phic = cluster->phi();
00125   
00126   double etah = hit->eta();
00127   double phih = hit->phi();
00128 
00129   return rDist(etac, phic, etah, phih);
00130 }
00131 
00132 
00133 double SimG4HcalHitJetFinder::rDist(const double etac,const double phic,
00134                                     const double etah,const double phih) const{
00135 
00136   double delta_eta = etac - etah;
00137   double delta_phi = phic - phih;
00138 
00139   if (phic < phih)      delta_phi = phih - phic;
00140   if (delta_phi > M_PI) delta_phi = 2*M_PI - delta_phi;
00141 
00142   double tmp = sqrt(delta_eta * delta_eta + delta_phi * delta_phi);
00143 
00144   LogDebug("ValidHcal") << "HcalHitJetFinder::rDist:\n"
00145                         << " Clus. eta, phi = " << etac << " " << phic << "\n"
00146                         << " hit   eta, phi = " << etah << " " << phih 
00147                         << " rDist = " << tmp;
00148 
00149   return tmp;
00150 }
00151