Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-08-10 22:49:13 

0001 #include "Calibration/IsolatedParticles/interface/MatchingSimTrack.h"
0002 #include "Calibration/IsolatedParticles/interface/eHCALMatrix.h"
0003 #include "Calibration/IsolatedParticles/interface/eECALMatrix.h"
0004 #include "Calibration/IsolatedParticles/interface/FindCaloHitCone.h"
0005 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 
0008 #include <algorithm>
0009 #include <sstream>
0010 
0011 namespace spr {
0012 
0013   template <typename T>
0014   std::map<std::string, double> eECALSimInfo(const edm::Event& iEvent,
0015                                              const DetId& det,
0016                                              const CaloGeometry* geo,
0017                                              const CaloTopology* caloTopology,
0018                                              edm::Handle<T>& hitsEB,
0019                                              edm::Handle<T>& hitsEE,
0020                                              edm::Handle<edm::SimTrackContainer>& SimTk,
0021                                              edm::Handle<edm::SimVertexContainer>& SimVtx,
0022                                              const reco::Track* pTrack,
0023                                              TrackerHitAssociator& associate,
0024                                              int ieta,
0025                                              int iphi,
0026                                              double timeCut,
0027                                              bool debug) {
0028     spr::caloSimInfo info;
0029     spr::eECALSimInfo(iEvent,
0030                       det,
0031                       geo,
0032                       caloTopology,
0033                       hitsEB,
0034                       hitsEE,
0035                       SimTk,
0036                       SimVtx,
0037                       pTrack,
0038                       associate,
0039                       ieta,
0040                       iphi,
0041                       info,
0042                       timeCut,
0043                       debug);
0044     // return a map of matching type and energy of SimHits
0045     return spr::eCaloSimInfo(info);
0046   }
0047 
0048   template <typename T>
0049   void eECALSimInfo(const edm::Event& iEvent,
0050                     const DetId& det,
0051                     const CaloGeometry* geo,
0052                     const CaloTopology* caloTopology,
0053                     edm::Handle<T>& hitsEB,
0054                     edm::Handle<T>& hitsEE,
0055                     edm::Handle<edm::SimTrackContainer>& SimTk,
0056                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0057                     const reco::Track* pTrack,
0058                     TrackerHitAssociator& associate,
0059                     int ieta,
0060                     int iphi,
0061                     spr::caloSimInfo& info,
0062                     double timeCut,
0063                     bool debug) {
0064     if (debug)
0065       edm::LogVerbatim("IsoTrack") << "Processing eECALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0066                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0067                                    << pTrack->recHitsSize();
0068 
0069     //matching SimTrack
0070     edm::SimTrackContainer::const_iterator trkInfo =
0071         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0072 
0073     //vector of Ecal cells in NxN
0074     std::vector<DetId> vdets = spr::matrixECALIds(det, ieta, iphi, geo, caloTopology, debug);
0075 
0076     //Construct a struct with matching type and energy of SimHits
0077     spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, info, timeCut, debug);
0078   }
0079 
0080   template <typename T>
0081   void eECALSimInfo(const edm::Event& iEvent,
0082                     const DetId& det,
0083                     const CaloGeometry* geo,
0084                     const CaloTopology* caloTopology,
0085                     edm::Handle<T>& hitsEB,
0086                     edm::Handle<T>& hitsEE,
0087                     edm::Handle<edm::SimTrackContainer>& SimTk,
0088                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0089                     const reco::Track* pTrack,
0090                     TrackerHitAssociator& associate,
0091                     int ietaE,
0092                     int ietaW,
0093                     int iphiN,
0094                     int iphiS,
0095                     spr::caloSimInfo& info,
0096                     double timeCut,
0097                     bool debug) {
0098     if (debug)
0099       edm::LogVerbatim("IsoTrack") << "Processing eECALSimInfo " << ietaE + ietaW + 1 << "x" << iphiN + iphiS + 1
0100                                    << "\ntrkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0101                                    << pTrack->recHitsSize();
0102 
0103     //matching SimTrack
0104     edm::SimTrackContainer::const_iterator trkInfo =
0105         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0106 
0107     //vector of Ecal cells in NxN
0108     std::vector<DetId> vdets = spr::matrixECALIds(det, ietaE, ietaW, iphiN, iphiS, geo, caloTopology, debug);
0109 
0110     // return a map of matching type and energy of SimHits
0111     spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, info, timeCut, debug);
0112   }
0113 
0114   template <typename T>
0115   std::map<std::string, double> eHCALSimInfo(const edm::Event& iEvent,
0116                                              const HcalTopology* topology,
0117                                              const DetId& det,
0118                                              const CaloGeometry* geo,
0119                                              edm::Handle<T>& hits,
0120                                              edm::Handle<edm::SimTrackContainer>& SimTk,
0121                                              edm::Handle<edm::SimVertexContainer>& SimVtx,
0122                                              const reco::Track* pTrack,
0123                                              TrackerHitAssociator& associate,
0124                                              int ieta,
0125                                              int iphi,
0126                                              double timeCut,
0127                                              bool includeHO,
0128                                              bool debug) {
0129     spr::caloSimInfo info;
0130     spr::eHCALSimInfo(
0131         iEvent, topology, det, geo, hits, SimTk, SimVtx, pTrack, associate, ieta, iphi, info, timeCut, includeHO, debug);
0132     // return a map of matching type and energy of SimHits
0133     return spr::eCaloSimInfo(info);
0134   }
0135 
0136   template <typename T>
0137   void eHCALSimInfo(const edm::Event& iEvent,
0138                     const HcalTopology* topology,
0139                     const DetId& det,
0140                     const CaloGeometry* geo,
0141                     edm::Handle<T>& hits,
0142                     edm::Handle<edm::SimTrackContainer>& SimTk,
0143                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0144                     const reco::Track* pTrack,
0145                     TrackerHitAssociator& associate,
0146                     int ieta,
0147                     int iphi,
0148                     spr::caloSimInfo& info,
0149                     double timeCut,
0150                     bool includeHO,
0151                     bool debug) {
0152     if (debug)
0153       edm::LogVerbatim("IsoTrack") << "Processing eHCALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0154                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0155                                    << pTrack->recHitsSize();
0156 
0157     // get the matching SimTrack pointer
0158     edm::SimTrackContainer::const_iterator trkInfo =
0159         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0160 
0161     // get the hits in Hcal in NxN around det
0162     std::vector<typename T::const_iterator> hit;
0163     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0164 
0165     spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, info, timeCut, includeHO, debug);
0166   }
0167 
0168   template <typename T>
0169   void eHCALSimInfo(const edm::Event& iEvent,
0170                     const HcalTopology* topology,
0171                     const DetId& det,
0172                     const CaloGeometry* geo,
0173                     edm::Handle<T>& hits,
0174                     edm::Handle<edm::SimTrackContainer>& SimTk,
0175                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0176                     const reco::Track* pTrack,
0177                     TrackerHitAssociator& associate,
0178                     int ietaE,
0179                     int ietaW,
0180                     int iphiN,
0181                     int iphiS,
0182                     spr::caloSimInfo& info,
0183                     double timeCut,
0184                     bool includeHO,
0185                     bool debug) {
0186     if (debug)
0187       edm::LogVerbatim("IsoTrack") << "Processing eHCALSimInfo " << ietaE + ietaW + 1 << "x" << iphiN + iphiS + 1
0188                                    << "\ntrkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0189                                    << pTrack->recHitsSize();
0190 
0191     // get the matching SimTrack pointer
0192     edm::SimTrackContainer::const_iterator trkInfo =
0193         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0194 
0195     // get the hits in Hcal in NxN around det
0196     std::vector<typename T::const_iterator> hit;
0197     spr::hitHCALmatrixTotal(topology, det, hits, ietaE, ietaW, iphiN, iphiS, hit, includeHO, debug);
0198 
0199     spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, info, timeCut, includeHO, debug);
0200   }
0201 
0202   // NxN method modified to return vector of
0203   // simhit multiplicities by type of particle
0204   template <typename T>
0205   std::map<std::string, double> eHCALSimInfo(const edm::Event& iEvent,
0206                                              const HcalTopology* topology,
0207                                              const DetId& det,
0208                                              edm::Handle<T>& hits,
0209                                              edm::Handle<edm::SimTrackContainer>& SimTk,
0210                                              edm::Handle<edm::SimVertexContainer>& SimVtx,
0211                                              const reco::Track* pTrack,
0212                                              TrackerHitAssociator& associate,
0213                                              int ieta,
0214                                              int iphi,
0215                                              std::vector<int>& multiplicityVector,
0216                                              bool debug) {
0217     if (debug) {
0218       edm::LogVerbatim("IsoTrack") << "Processing HcalSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1;
0219       edm::LogVerbatim("IsoTrack") << "trkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0220                                    << pTrack->recHitsSize();
0221     }
0222 
0223     // Get SimTrack matched to RecoTrack
0224     edm::SimTrackContainer::const_iterator trkInfo =
0225         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0226 
0227     // Get list of simhits in matrix ("hits" is list of pcalo simhits)
0228     std::vector<typename T::const_iterator> hit;
0229     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, false, debug);
0230 
0231     return eCaloSimInfo(hits, SimTk, SimVtx, hit, trkInfo, multiplicityVector, debug);
0232   }
0233 
0234   // eHCALSimInfo for Cone
0235   template <typename T>
0236   std::map<std::string, double> eHCALSimInfoCone(const edm::Event& iEvent,
0237                                                  edm::Handle<T>& hits,
0238                                                  edm::Handle<edm::SimTrackContainer>& SimTk,
0239                                                  edm::Handle<edm::SimVertexContainer>& SimVtx,
0240                                                  const reco::Track* pTrack,
0241                                                  TrackerHitAssociator& associate,
0242                                                  const CaloGeometry* geo,
0243                                                  const GlobalPoint& hpoint1,
0244                                                  const GlobalPoint& point1,
0245                                                  double dR,
0246                                                  const GlobalVector& trackMom,
0247                                                  std::vector<int>& multiplicityVector) {
0248     bool debug_dummy = false;
0249 
0250     // Get SimTrack matched to RecoTrack
0251     edm::SimTrackContainer::const_iterator trkInfo =
0252         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug_dummy);
0253 
0254     //    spr::hitHCALmatrix(topology,det,hits,ieta,iphi,debug);
0255 
0256     // Get list of simhits in matrix ("hits" is list of pcalo simhits)
0257     std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom);
0258 
0259     return eCaloSimInfo(hits, SimTk, SimVtx, hit, trkInfo, multiplicityVector, debug_dummy);
0260   }
0261 
0262   // eHCALSimInfo for Cone that only clusters simhits in Hcal cells
0263   // that contained a rechit
0264   template <typename T>
0265   std::map<std::string, double> eHCALSimInfoCone(const edm::Event& iEvent,
0266                                                  edm::Handle<T>& hits,
0267                                                  edm::Handle<edm::SimTrackContainer>& SimTk,
0268                                                  edm::Handle<edm::SimVertexContainer>& SimVtx,
0269                                                  const reco::Track* pTrack,
0270                                                  TrackerHitAssociator& associate,
0271                                                  const CaloGeometry* geo,
0272                                                  const GlobalPoint& hpoint1,
0273                                                  const GlobalPoint& point1,
0274                                                  double dR,
0275                                                  const GlobalVector& trackMom,
0276                                                  std::vector<int>& multiplicityVector,
0277                                                  std::vector<DetId>& coneRecHitDetIds) {
0278     bool debug_dummy = false;
0279 
0280     // Get SimTrack matched to RecoTrack
0281     edm::SimTrackContainer::const_iterator trkInfo =
0282         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug_dummy);
0283 
0284     //    spr::hitHCALmatrix(topology,det,hits,ieta,iphi,debug);
0285 
0286     // Get list of simhits in matrix ("hits" is list of pcalo simhits)
0287     std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom);
0288 
0289     // Loop over hits and only keep those whose detid matches a detid in
0290     // coneRecHitDetIds:
0291 
0292     std::vector<typename T::const_iterator> matchedhit;
0293     for (int ihit = 0; ihit < hit.size(); ihit++) {
0294       bool keephit = false;
0295       for (size_t idetid = 0; idetid < coneRecHitDetIds.size(); idetid++) {
0296         if (hit.at(ihit)->id() == coneRecHitDetIds.at(idetid)) {
0297           keephit = true;
0298           break;
0299         }
0300       }
0301       if (keephit)
0302         matchedhit.push_back(hit.at(ihit));
0303     }
0304 
0305     return eCaloSimInfo(hits, SimTk, SimVtx, matchedhit, trkInfo, multiplicityVector, debug_dummy);
0306   }
0307 
0308   template <typename T>
0309   void eCaloSimInfo(std::vector<DetId> vdets,
0310                     const CaloGeometry* geo,
0311                     edm::Handle<T>& hitsEB,
0312                     edm::Handle<T>& hitsEE,
0313                     edm::Handle<edm::SimTrackContainer>& SimTk,
0314                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0315                     edm::SimTrackContainer::const_iterator trkInfo,
0316                     spr::caloSimInfo& info,
0317                     double timeCut,
0318                     bool debug) {
0319     std::vector<typename T::const_iterator> hitEB, hitEE, hit;
0320     for (unsigned int i1 = 0; i1 < vdets.size(); i1++) {
0321       if (vdets[i1] != DetId(0)) {
0322         if (vdets[i1].subdetId() == EcalBarrel) {
0323           hit = spr::findHit(hitsEB, vdets[i1]);
0324           for (unsigned int i2 = 0; i2 < hit.size(); i2++) {
0325             if (std::count(hitEB.begin(), hitEB.end(), hit[i2]) == 0)
0326               hitEB.push_back(hit[i2]);
0327           }
0328         } else if (vdets[i1].subdetId() == EcalEndcap) {
0329           hit = spr::findHit(hitsEE, vdets[i1]);
0330           for (unsigned int i2 = 0; i2 < hit.size(); i2++) {
0331             if (std::count(hitEE.begin(), hitEE.end(), hit[i2]) == 0)
0332               hitEE.push_back(hit[i2]);
0333           }
0334         }
0335       }
0336     }
0337 
0338     spr::caloSimInfo eeInfo;
0339     spr::eCaloSimInfo(geo, hitsEB, SimTk, SimVtx, hitEB, trkInfo, info, timeCut, false, debug);
0340     spr::eCaloSimInfo(geo, hitsEE, SimTk, SimVtx, hitEE, trkInfo, eeInfo, timeCut, false, debug);
0341 
0342     if (info.pdgMatched == 0)
0343       info.pdgMatched = eeInfo.pdgMatched;
0344     info.eMatched += eeInfo.eMatched;
0345     info.eGamma += eeInfo.eGamma;
0346     info.eNeutralHad += eeInfo.eNeutralHad;
0347     info.eChargedHad += eeInfo.eChargedHad;
0348     info.eRest += eeInfo.eRest;
0349     info.eTotal += eeInfo.eTotal;
0350     if (debug) {
0351       edm::LogVerbatim("IsoTrack") << " energyMatched " << info.eMatched << " energyGamma   " << info.eGamma
0352                                    << " energyNeutral " << info.eNeutralHad << " energyCharged " << info.eChargedHad
0353                                    << " energyRest    " << info.eRest << " energyTotal   " << info.eTotal;
0354     }
0355   }
0356 
0357   template <typename T>
0358   void eCaloSimInfo(const CaloGeometry* geo,
0359                     edm::Handle<T>& hits,
0360                     edm::Handle<edm::SimTrackContainer>& SimTk,
0361                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0362                     std::vector<typename T::const_iterator> hit,
0363                     edm::SimTrackContainer::const_iterator trkInfo,
0364                     spr::caloSimInfo& info,
0365                     double timeCut,
0366                     bool includeHO,
0367                     bool debug) {
0368     int matchedId = 0;  //pdgid
0369     if (debug) {
0370       if (trkInfo != SimTk->end())
0371         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  matchSimTrk:" << trkInfo->trackId() << " matchedId "
0372                                      << trkInfo->type();
0373       else
0374         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  not valid track pointer";
0375     }
0376 
0377     if (trkInfo != SimTk->end()) {
0378       unsigned int matchSimTrk = trkInfo->trackId();
0379       matchedId = trkInfo->type();  //pdgid
0380 
0381       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0382         DetId id_ = (DetId)(hit[ihit]->id());
0383         double tof = timeOfFlight(id_, geo, debug);
0384         bool ok = true;
0385         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
0386           ok = false;
0387         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
0388           ok = false;
0389         if ((hit[ihit]->time() <= (tof + timeCut)) && ok) {
0390           // if the hitId matches with matching trackId
0391           if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0392             info.eMatched += hit[ihit]->energy();
0393           } else {
0394             // trace back the history and check the pdgId of origin SimTrack of SimHit
0395             bool found = false;
0396             for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0397               if (hit[ihit]->geantTrackId() == (int)simTrkItr->trackId()) {
0398                 found = true;
0399                 bool match = spr::validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, debug);
0400                 if (match) {
0401                   info.eMatched += hit[ihit]->energy();
0402                 } else {
0403                   edm::SimTrackContainer::const_iterator parentItr =
0404                       spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0405                   if (debug) {
0406                     if (parentItr != SimTk->end())
0407                       edm::LogVerbatim("IsoTrack")
0408                           << "original parent of " << simTrkItr->trackId() << " " << parentItr->trackId() << ", "
0409                           << parentItr->type() << " Energy " << hit[ihit]->energy();
0410                     else
0411                       edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId()
0412                                                    << " not found; Energy " << hit[ihit]->energy();
0413                   }
0414                   if (parentItr == SimTk->end()) {
0415                     info.eRest += hit[ihit]->energy();
0416                   } else if (parentItr->type() == 22) {
0417                     info.eGamma += hit[ihit]->energy();
0418                   } else if ((int)parentItr->charge() == 0) {
0419                     info.eNeutralHad += hit[ihit]->energy();
0420                   } else
0421                     info.eChargedHad += hit[ihit]->energy();
0422                 }
0423                 break;
0424               }
0425             }
0426             if (!found)
0427               info.eRest += hit[ihit]->energy();
0428             if (debug)
0429               edm::LogVerbatim("IsoTrack") << "Hit " << ihit << ": " << *hit[ihit];
0430           }
0431         }
0432       }  // loop over hits
0433     }
0434 
0435     info.eTotal = info.eMatched + info.eGamma + info.eNeutralHad + info.eChargedHad + info.eRest;
0436     info.pdgMatched = matchedId;
0437     if (debug) {
0438       edm::LogVerbatim("IsoTrack") << " Energy (matched) " << info.eMatched << " (gamma) " << info.eGamma
0439                                    << " (neutral) " << info.eNeutralHad << " (charged) " << info.eChargedHad
0440                                    << " (rest) " << info.eRest << " (total) " << info.eTotal;
0441     }
0442   }
0443 
0444   inline std::map<std::string, double> eCaloSimInfo(spr::caloSimInfo& info) {
0445     std::map<std::string, double> simInfo;
0446     simInfo.insert(std::pair<std::string, double>("eMatched", info.eMatched));
0447     simInfo.insert(std::pair<std::string, double>("pdgMatched", info.pdgMatched));
0448     simInfo.insert(std::pair<std::string, double>("eGamma", info.eGamma));
0449     simInfo.insert(std::pair<std::string, double>("eNeutralHad", info.eNeutralHad));
0450     simInfo.insert(std::pair<std::string, double>("eChargedHad", info.eChargedHad));
0451     simInfo.insert(std::pair<std::string, double>("eRest", info.eRest));
0452     simInfo.insert(std::pair<std::string, double>("eTotal", info.eTotal));
0453     return simInfo;
0454   }
0455 
0456   // Returns total energy of CaloSimHits which originate from the matching SimTrack
0457   template <typename T>
0458   double eCaloSimInfo(const edm::Event& iEvent,
0459                       const CaloGeometry* geo,
0460                       edm::Handle<T>& hits,
0461                       edm::Handle<edm::SimTrackContainer>& SimTk,
0462                       edm::Handle<edm::SimVertexContainer>& SimVtx,
0463                       const reco::Track* pTrack,
0464                       TrackerHitAssociator& associate,
0465                       double timeCut,
0466                       bool includeHO,
0467                       bool debug) {
0468     //
0469     std::vector<int> matchedId = spr::matchedSimTrackId(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0470     double energySum = 0.0;
0471     if (debug) {
0472       std::ostringstream st1;
0473       for (unsigned int it = 0; it < matchedId.size(); ++it)
0474         st1 << " " << matchedId[it];
0475       edm::LogVerbatim("IsoTrack") << "eCaloSimInfo:: Found " << matchedId.size()
0476                                    << " track IDs originating from the current track" << st1.str();
0477     }
0478     if (!matchedId.empty()) {
0479       typename T::const_iterator ihit;
0480       for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
0481         //cut on time
0482         DetId id_ = (DetId)(ihit->id());
0483         double tof = timeOfFlight(id_, geo, debug);
0484         bool ok = true;
0485         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
0486           ok = false;
0487         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
0488           ok = false;
0489         if ((ihit->time() <= (tof + timeCut)) && ok) {
0490           int id = ihit->geantTrackId();
0491           bool found = false;
0492           for (unsigned int it = 0; it < matchedId.size(); ++it) {
0493             if (id == matchedId[it]) {
0494               found = true;
0495               break;
0496             }
0497           }
0498           if (found) {
0499             energySum += ihit->energy();
0500             if (debug)
0501               edm::LogVerbatim("IsoTrack") << "Hit " << *ihit;
0502           }
0503         }
0504       }
0505     }
0506     return energySum;
0507   }
0508 
0509   template <typename T>
0510   double eCaloSimInfo(const edm::Event& iEvent,
0511                       const CaloGeometry* geo,
0512                       edm::Handle<T>& hitsEB,
0513                       edm::Handle<T>& hitsEE,
0514                       edm::Handle<edm::SimTrackContainer>& SimTk,
0515                       edm::Handle<edm::SimVertexContainer>& SimVtx,
0516                       const reco::Track* pTrack,
0517                       TrackerHitAssociator& associate,
0518                       double timeCut,
0519                       bool debug) {
0520     //
0521     std::vector<int> matchedId = spr::matchedSimTrackId(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0522     double energySum = 0.0;
0523     if (debug) {
0524       std::ostringstream st1;
0525       for (unsigned int it = 0; it < matchedId.size(); ++it)
0526         st1 << " " << matchedId[it];
0527       edm::LogVerbatim("IsoTrack") << "eCaloSimInfo:: Found " << matchedId.size()
0528                                    << " track IDs originating from the current track" << st1.str();
0529     }
0530     if (!matchedId.empty()) {
0531       typename T::const_iterator ihit;
0532       for (ihit = hitsEB->begin(); ihit != hitsEB->end(); ihit++) {
0533         //cut on time
0534         double tof = timeOfFlight((DetId)(ihit->id()), geo, debug);
0535         if (ihit->time() <= (tof + timeCut)) {
0536           int id = ihit->geantTrackId();
0537           bool found = false;
0538           for (unsigned int it = 0; it < matchedId.size(); ++it) {
0539             if (id == matchedId[it]) {
0540               found = true;
0541               break;
0542             }
0543           }
0544           if (found) {
0545             energySum += ihit->energy();
0546             if (debug)
0547               edm::LogVerbatim("IsoTrack") << "Hit " << *ihit;
0548           }
0549         }
0550       }
0551       for (ihit = hitsEE->begin(); ihit != hitsEE->end(); ihit++) {
0552         //cut on time
0553         double tof = timeOfFlight((DetId)(ihit->id()), geo, debug);
0554         if (ihit->time() <= (tof + timeCut)) {
0555           int id = ihit->geantTrackId();
0556           bool found = false;
0557           for (unsigned int it = 0; it < matchedId.size(); ++it) {
0558             if (id == matchedId[it]) {
0559               found = true;
0560               break;
0561             }
0562           }
0563           if (found) {
0564             energySum += ihit->energy();
0565             if (debug)
0566               edm::LogVerbatim("IsoTrack") << "Hit " << *ihit;
0567           }
0568         }
0569       }
0570     }
0571     return energySum;
0572   }
0573 
0574   template <typename T>
0575   std::map<std::string, double> eCaloSimInfo(edm::Handle<T>& hits,
0576                                              edm::Handle<edm::SimTrackContainer>& SimTk,
0577                                              edm::Handle<edm::SimVertexContainer>& SimVtx,
0578                                              std::vector<typename T::const_iterator> hit,
0579                                              edm::SimTrackContainer::const_iterator trkInfo,
0580                                              std::vector<int>& multiplicityVector,
0581                                              bool debug) {
0582     int matchedId = 0;  //pdgid
0583 
0584     if (debug) {
0585       if (trkInfo != SimTk->end()) {
0586         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  matchSimTrk:" << trkInfo->trackId() << " matchedId "
0587                                      << trkInfo->type();
0588       } else {
0589         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  not valid track pointer";
0590       }
0591     }
0592 
0593     // Find sets of hcal cells with energy
0594     // Use std::set to avoid duplicates
0595 
0596     std::set<int> uniqueIds_;
0597     std::set<int> uniqueIds_matched;
0598     std::set<int> uniqueIds_total;
0599     std::set<int> uniqueIds_neut;
0600     std::set<int> uniqueIds_char;
0601     std::set<int> uniqueIds_gamma;
0602     std::set<int> uniqueIds_rest;
0603 
0604     // Counters for SimHits
0605     int nMatched = 0, nRest = 0, nGamma = 0, nNeutralHad = 0, nChargedHad = 0;
0606 
0607     double energySum = 0.0;
0608     double energyGamma = 0.0, energyNeutral = 0.0, energyCharged = 0.0, energyRest = 0.0;
0609 
0610     // What is this checking?
0611     if (trkInfo != SimTk->end()) {
0612       unsigned int matchSimTrk = trkInfo->trackId();
0613       matchedId = trkInfo->type();  //pdgid
0614 
0615       // Loop over hits
0616       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0617         // For simhit tower multiplicity, we need to avoid double
0618         // counting towers with multiple simhits:
0619         // This will break if called for pcalohits from the ecal!!
0620         HcalDetId detId(hit[ihit]->id());
0621         int ieta = detId.ieta();
0622         int iphi = detId.iphi();
0623         int uniqueId = 100 * ieta + iphi;
0624 
0625         if (hit[ihit]->time() > 150.0)
0626           continue;
0627 
0628         if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0629           energySum += spr::getEnergy(hit[ihit]);
0630           uniqueIds_matched.insert(uniqueId);
0631           uniqueIds_total.insert(uniqueId);
0632           nMatched++;
0633         } else {
0634           bool found = false;
0635           for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0636             if (hit[ihit]->geantTrackId() == (int)simTrkItr->trackId()) {
0637               found = true;
0638               bool match = spr::validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, debug);
0639               if (match) {
0640                 energySum += spr::getEnergy(hit[ihit]);
0641                 uniqueIds_matched.insert(uniqueId);
0642                 uniqueIds_total.insert(uniqueId);
0643                 nMatched++;
0644               } else {
0645                 edm::SimTrackContainer::const_iterator parentItr = spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0646                 if (debug) {
0647                   if (parentItr != SimTk->end()) {
0648                     edm::LogVerbatim("IsoTrack")
0649                         << "original parent of " << simTrkItr->trackId() << " " << parentItr->trackId() << ", "
0650                         << parentItr->type() << " Energy " << spr::getEnergy(hit[ihit]);
0651                   } else {
0652                     edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId()
0653                                                  << " not found; Energy " << spr::getEnergy(hit[ihit]);
0654                   }
0655                 }
0656                 if (parentItr == SimTk->end()) {
0657                   energyRest += spr::getEnergy(hit[ihit]);
0658                   uniqueIds_rest.insert(uniqueId);
0659                   uniqueIds_total.insert(uniqueId);
0660                   nRest++;
0661                 } else if (parentItr->type() == 22) {
0662                   energyGamma += spr::getEnergy(hit[ihit]);
0663                   uniqueIds_gamma.insert(uniqueId);
0664                   uniqueIds_total.insert(uniqueId);
0665                   nGamma++;
0666                 } else if ((int)parentItr->charge() == 0) {
0667                   energyNeutral += spr::getEnergy(hit[ihit]);
0668                   uniqueIds_neut.insert(uniqueId);
0669                   uniqueIds_total.insert(uniqueId);
0670                   nNeutralHad++;
0671                 } else {
0672                   energyCharged += spr::getEnergy(hit[ihit]);
0673                   uniqueIds_char.insert(uniqueId);
0674                   uniqueIds_total.insert(uniqueId);
0675                   nChargedHad++;
0676                 }
0677               }
0678               break;
0679             }
0680           }  // for (simTrkItr = ...)
0681 
0682           if (!found) {
0683             energyRest += spr::getEnergy(hit[ihit]);
0684             uniqueIds_rest.insert(uniqueId);
0685             uniqueIds_total.insert(uniqueId);
0686             nRest++;
0687           }
0688 
0689           if (debug)
0690             edm::LogVerbatim("IsoTrack") << "Hit " << ihit << ": " << *hit[ihit];
0691         }  // else condition, i.e. (hit[ihit]->geantTrackId() != (int)matchSimTrk )
0692       }    // loop over hits
0693     }      // if (trkInfo != SimTk->end())
0694 
0695     double energyTot = energySum + energyGamma + energyNeutral + energyCharged + energyRest;
0696     multiplicityVector.push_back(uniqueIds_matched.size());
0697     multiplicityVector.push_back(uniqueIds_total.size());
0698     multiplicityVector.push_back(uniqueIds_neut.size());
0699     multiplicityVector.push_back(uniqueIds_char.size());
0700     multiplicityVector.push_back(uniqueIds_gamma.size());
0701     multiplicityVector.push_back(uniqueIds_rest.size());
0702 
0703     std::map<std::string, double> simInfo;
0704     simInfo.insert(std::pair<std::string, double>("eMatched", energySum));
0705     simInfo.insert(std::pair<std::string, double>("pdgMatched", matchedId));
0706     simInfo.insert(std::pair<std::string, double>("eGamma", energyGamma));
0707     simInfo.insert(std::pair<std::string, double>("eNeutralHad", energyNeutral));
0708     simInfo.insert(std::pair<std::string, double>("eChargedHad", energyCharged));
0709     simInfo.insert(std::pair<std::string, double>("eRest", energyRest));
0710     simInfo.insert(std::pair<std::string, double>("eTotal", energyTot));
0711     if (debug) {
0712       edm::LogVerbatim("IsoTrack") << " energySum " << energySum << "  energyGamma " << energyGamma
0713                                    << "  energyNeutral " << energyNeutral << "  energyCharged " << energyCharged
0714                                    << "  energyRest " << energyRest << "  energyTot " << energyTot;
0715     }
0716     return simInfo;
0717   }  // eCaloSimInfo
0718 
0719   template <typename T>
0720   std::map<std::string, double> eECALSimInfo(const edm::Event& iEvent,
0721                                              CaloNavigator<DetId>& navigator,
0722                                              const CaloGeometry* geo,
0723                                              edm::Handle<T>& hits,
0724                                              edm::Handle<edm::SimTrackContainer>& SimTk,
0725                                              edm::Handle<edm::SimVertexContainer>& SimVtx,
0726                                              const reco::Track* pTrack,
0727                                              TrackerHitAssociator& associate,
0728                                              int ieta,
0729                                              int iphi,
0730                                              double timeCut,
0731                                              bool debug) {
0732     if (debug)
0733       edm::LogVerbatim("IsoTrack") << "Processing eECALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0734                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0735                                    << pTrack->recHitsSize();
0736 
0737     //matching SimTrack
0738     edm::SimTrackContainer::const_iterator trkInfo =
0739         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0740     //vector of Ecal hits in NxN
0741     std::vector<typename T::const_iterator> ecalHits;
0742     spr::hitECALmatrix(navigator, hits, ieta, iphi, ecalHits, debug);
0743 
0744     // return a map of matching type and energy of SimHits
0745     spr::caloSimInfo info;
0746     spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, ecalHits, trkInfo, info, timeCut, false, debug);
0747     return spr::eCaloSimInfo(info);
0748   }
0749 
0750   template <typename T>
0751   std::map<std::string, double> eECALSimInfoTotal(const edm::Event& iEvent,
0752                                                   const DetId& det,
0753                                                   const CaloGeometry* geo,
0754                                                   const CaloTopology* caloTopology,
0755                                                   edm::Handle<T>& hitsEB,
0756                                                   edm::Handle<T>& hitsEE,
0757                                                   edm::Handle<edm::SimTrackContainer>& SimTk,
0758                                                   edm::Handle<edm::SimVertexContainer>& SimVtx,
0759                                                   const reco::Track* pTrack,
0760                                                   TrackerHitAssociator& associate,
0761                                                   int ieta,
0762                                                   int iphi,
0763                                                   int itry,
0764                                                   double timeCut,
0765                                                   bool debug) {
0766     if (debug)
0767       edm::LogVerbatim("IsoTrack") << "Processing eECALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0768                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0769                                    << pTrack->recHitsSize();
0770 
0771     //matching SimTrack
0772     edm::SimTrackContainer::const_iterator trkInfo =
0773         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0774 
0775     spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);
0776     std::map<std::string, double> simInfo;
0777 
0778     if (itry >= 0) {
0779       //vector of Ecal cells in NxN
0780       std::vector<DetId> vdets = spr::matrixECALIds(
0781           det, etaphi.ietaE[itry], etaphi.ietaW[itry], etaphi.iphiN[itry], etaphi.iphiS[itry], geo, caloTopology, debug);
0782       // get a map of matching type and energy of SimHits
0783       spr::caloSimInfo info;
0784       spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, info, timeCut, debug);
0785       simInfo = spr::eCaloSimInfo(info);
0786     } else {
0787       int itrym = 0;
0788       std::vector<DetId> vdets = spr::matrixECALIds(det,
0789                                                     etaphi.ietaE[itrym],
0790                                                     etaphi.ietaW[itrym],
0791                                                     etaphi.iphiN[itrym],
0792                                                     etaphi.iphiS[itrym],
0793                                                     geo,
0794                                                     caloTopology,
0795                                                     debug);
0796       // get a map of matching type and energy of SimHits
0797       spr::caloSimInfo info;
0798       spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, info, timeCut, debug);
0799       simInfo = spr::eCaloSimInfo(info);
0800       for (int itrys = 1; itrys < etaphi.ntrys; ++itrys) {
0801         vdets = spr::matrixECALIds(det,
0802                                    etaphi.ietaE[itrys],
0803                                    etaphi.ietaW[itrys],
0804                                    etaphi.iphiN[itrys],
0805                                    etaphi.iphiS[itrys],
0806                                    geo,
0807                                    caloTopology,
0808                                    debug);
0809         // get a map of matching type and energy of SimHits
0810         spr::caloSimInfo infox;
0811         spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, infox, timeCut, debug);
0812         std::map<std::string, double> simInfoX = spr::eCaloSimInfo(infox);
0813         if (simInfoX["eTotal"] > simInfo["eTotal"]) {
0814           simInfo = simInfoX;
0815           itrym = itrys;
0816         }
0817       }
0818     }
0819     return simInfo;
0820   }
0821 
0822   template <typename T>
0823   std::map<std::string, double> eHCALSimInfoTotal(const edm::Event& iEvent,
0824                                                   const HcalTopology* topology,
0825                                                   const DetId& det,
0826                                                   const CaloGeometry* geo,
0827                                                   edm::Handle<T>& hits,
0828                                                   edm::Handle<edm::SimTrackContainer>& SimTk,
0829                                                   edm::Handle<edm::SimVertexContainer>& SimVtx,
0830                                                   const reco::Track* pTrack,
0831                                                   TrackerHitAssociator& associate,
0832                                                   int ieta,
0833                                                   int iphi,
0834                                                   int itry,
0835                                                   double timeCut,
0836                                                   bool includeHO,
0837                                                   bool debug) {
0838     if (debug)
0839       edm::LogVerbatim("IsoTrack") << "Processing eHCALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0840                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0841                                    << pTrack->recHitsSize();
0842 
0843     // get the matching SimTrack pointer
0844     edm::SimTrackContainer::const_iterator trkInfo =
0845         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0846 
0847     spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);
0848     std::map<std::string, double> simInfo;
0849 
0850     if (itry >= 0) {
0851       // get the hits in Hcal in NxN around det
0852       std::vector<typename T::const_iterator> hit;
0853       spr::hitHCALmatrixTotal(topology,
0854                               det,
0855                               hits,
0856                               etaphi.ietaE[itry],
0857                               etaphi.ietaW[itry],
0858                               etaphi.iphiN[itry],
0859                               etaphi.iphiS[itry],
0860                               hit,
0861                               false,
0862                               debug);
0863       spr::caloSimInfo info;
0864       spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, info, timeCut, includeHO, debug);
0865       simInfo = spr::eCaloSimInfo(info);
0866     } else {
0867       int itrym = 0;
0868       std::vector<typename T::const_iterator> hit;
0869       spr::hitHCALmatrixTotal(topology,
0870                               det,
0871                               hits,
0872                               etaphi.ietaE[itrym],
0873                               etaphi.ietaW[itrym],
0874                               etaphi.iphiN[itrym],
0875                               etaphi.iphiS[itrym],
0876                               hit,
0877                               includeHO,
0878                               debug);
0879       spr::caloSimInfo info;
0880       spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, info, timeCut, includeHO, debug);
0881       simInfo = spr::eCaloSimInfo(info);
0882       for (int itrys = 1; itrys < etaphi.ntrys; ++itrys) {
0883         hit.clear();
0884         spr::hitHCALmatrixTotal(topology,
0885                                 det,
0886                                 hits,
0887                                 etaphi.ietaE[itrys],
0888                                 etaphi.ietaW[itrys],
0889                                 etaphi.iphiN[itrys],
0890                                 etaphi.iphiS[itrys],
0891                                 hit,
0892                                 includeHO,
0893                                 debug);
0894         spr::caloSimInfo infox;
0895         spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, infox, timeCut, includeHO, debug);
0896         std::map<std::string, double> simInfoX = spr::eCaloSimInfo(infox);
0897         if (simInfoX["eTotal"] > simInfo["eTotal"]) {
0898           simInfo = simInfoX;
0899           itrym = itrys;
0900         }
0901       }
0902     }
0903     return simInfo;
0904   }
0905 
0906   template <typename T>
0907   spr::energyMap eHCALSimInfoMatrix(const edm::Event& iEvent,
0908                                     const HcalTopology* topology,
0909                                     const DetId& det,
0910                                     const CaloGeometry* geo,
0911                                     edm::Handle<T>& hits,
0912                                     edm::Handle<edm::SimTrackContainer>& SimTk,
0913                                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0914                                     const reco::Track* pTrack,
0915                                     TrackerHitAssociator& associate,
0916                                     int ieta,
0917                                     int iphi,
0918                                     double timeCut,
0919                                     bool includeHO,
0920                                     bool debug) {
0921     if (debug)
0922       edm::LogVerbatim("IsoTrack") << "Processing eHCALSimInfoMatrix " << 2 * ieta + 1 << "x" << 2 * iphi + 1
0923                                    << "\ntrkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0924                                    << pTrack->recHitsSize();
0925 
0926     // get the matching SimTrack pointer
0927     edm::SimTrackContainer::const_iterator trkInfo =
0928         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0929 
0930     // get the hits in Hcal in NxN around det
0931     std::vector<typename T::const_iterator> hit;
0932     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0933 
0934     return spr::caloSimInfoMatrix(geo, hits, SimTk, SimVtx, hit, trkInfo, timeCut, includeHO, debug);
0935   }
0936 
0937   template <typename T>
0938   spr::energyMap caloSimInfoMatrix(const CaloGeometry* geo,
0939                                    edm::Handle<T>& hits,
0940                                    edm::Handle<edm::SimTrackContainer>& SimTk,
0941                                    edm::Handle<edm::SimVertexContainer>& SimVtx,
0942                                    std::vector<typename T::const_iterator> hit,
0943                                    edm::SimTrackContainer::const_iterator trkInfo,
0944                                    double timeCut,
0945                                    bool includeHO,
0946                                    bool debug) {
0947     int matchedId = 0;  //pdgid
0948 
0949     if (debug) {
0950       if (trkInfo != SimTk->end())
0951         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  matchSimTrk:" << trkInfo->trackId() << " matchedId "
0952                                      << trkInfo->type();
0953       else
0954         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  not valid track pointer";
0955     }
0956 
0957     std::vector<std::pair<DetId, double> > detSum, detGamma, detCharged, detNeutral, detRest, detAll;
0958 
0959     if (trkInfo != SimTk->end()) {
0960       unsigned int matchSimTrk = trkInfo->trackId();
0961       matchedId = trkInfo->type();  //pdgid
0962 
0963       edm::SimTrackContainer::const_iterator simTrkItr;
0964 
0965       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0966         DetId id_ = (DetId)(hit[ihit]->id());
0967         double tof = timeOfFlight(id_, geo, debug);
0968         double energySum = 0, energyRest = 0;
0969         double energyGamma = 0, energyNeutral = 0, energyCharged = 0;
0970         bool ok = true;
0971         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
0972           ok = false;
0973         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
0974           ok = false;
0975         if ((hit[ihit]->time() <= (tof + timeCut)) && ok) {
0976           // if the hitId matches with matching trackId
0977           if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0978             energySum = hit[ihit]->energy();
0979           } else {
0980             // trace back the history and check the pdgId of origin SimTrack of SimHit
0981             bool found = false;
0982             for (simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0983               if (hit[ihit]->geantTrackId() == (int)simTrkItr->trackId()) {
0984                 found = true;
0985                 bool match = spr::validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, debug);
0986                 if (match) {
0987                   energySum = hit[ihit]->energy();
0988                 } else {
0989                   edm::SimTrackContainer::const_iterator parentItr =
0990                       spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0991                   if (debug) {
0992                     if (parentItr != SimTk->end())
0993                       edm::LogVerbatim("IsoTrack")
0994                           << "original parent of " << simTrkItr->trackId() << " " << parentItr->trackId() << ", "
0995                           << parentItr->type() << " Energy " << hit[ihit]->energy();
0996                     else
0997                       edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId()
0998                                                    << " not found; Energy " << hit[ihit]->energy();
0999                   }
1000                   if (parentItr == SimTk->end()) {
1001                     energyRest = hit[ihit]->energy();
1002                   } else if (parentItr->type() == 22) {
1003                     energyGamma = hit[ihit]->energy();
1004                   } else if ((int)parentItr->charge() == 0) {
1005                     energyNeutral = hit[ihit]->energy();
1006                   } else
1007                     energyCharged = hit[ihit]->energy();
1008                 }
1009                 break;
1010               }
1011             }
1012             if (!found)
1013               energyRest = hit[ihit]->energy();
1014             if (debug)
1015               edm::LogVerbatim("IsoTrack") << "Hit " << ihit << ": " << *hit[ihit];
1016             found = false;
1017             if (energySum > 0) {
1018               for (unsigned int k = 0; k < detSum.size(); k++) {
1019                 if (id_ == detSum[k].first) {
1020                   found = true;
1021                   detSum[k].second += energySum;
1022                   break;
1023                 }
1024               }
1025               if (!found)
1026                 detSum.push_back(std::pair<DetId, double>(id_, energySum));
1027             } else if (energyRest > 0) {
1028               for (unsigned int k = 0; k < detRest.size(); k++) {
1029                 if (id_ == detRest[k].first) {
1030                   found = true;
1031                   detRest[k].second += energyRest;
1032                   break;
1033                 }
1034               }
1035               if (!found)
1036                 detRest.push_back(std::pair<DetId, double>(id_, energyRest));
1037             } else if (energyGamma > 0) {
1038               for (unsigned int k = 0; k < detGamma.size(); k++) {
1039                 if (id_ == detGamma[k].first) {
1040                   found = true;
1041                   detGamma[k].second += energyGamma;
1042                   break;
1043                 }
1044               }
1045               if (!found)
1046                 detGamma.push_back(std::pair<DetId, double>(id_, energyGamma));
1047             } else if (energyCharged > 0) {
1048               for (unsigned int k = 0; k < detCharged.size(); k++) {
1049                 if (id_ == detCharged[k].first) {
1050                   found = true;
1051                   detCharged[k].second += energyCharged;
1052                   break;
1053                 }
1054               }
1055               if (!found)
1056                 detCharged.push_back(std::pair<DetId, double>(id_, energyCharged));
1057             } else if (energyNeutral > 0) {
1058               for (unsigned int k = 0; k < detNeutral.size(); k++) {
1059                 if (id_ == detNeutral[k].first) {
1060                   found = true;
1061                   detNeutral[k].second += energyNeutral;
1062                   break;
1063                 }
1064               }
1065               if (!found)
1066                 detNeutral.push_back(std::pair<DetId, double>(id_, energyNeutral));
1067             }
1068             found = false;
1069             double energyTot = energySum + energyGamma + energyNeutral + energyCharged + energyRest;
1070             for (unsigned int k = 0; k < detAll.size(); k++) {
1071               if (id_ == detAll[k].first) {
1072                 found = true;
1073                 detAll[k].second += energyTot;
1074                 break;
1075               }
1076             }
1077             if (!found)
1078               detAll.push_back(std::pair<DetId, double>(id_, energyTot));
1079           }
1080         }
1081       }  // loop over hits
1082     }
1083 
1084     spr::energyMap simInfo;
1085     simInfo.pdgId = matchedId;
1086     simInfo.matched = detSum;
1087     simInfo.gamma = detGamma;
1088     simInfo.charged = detCharged;
1089     simInfo.neutral = detNeutral;
1090     simInfo.rest = detRest;
1091     simInfo.all = detAll;
1092     if (debug) {
1093       edm::LogVerbatim("IsoTrack") << "CaloSimInfo:: for particle " << simInfo.pdgId << "\n"
1094                                    << "All detIds " << detAll.size();
1095       for (unsigned int k = 0; k < detAll.size(); ++k)
1096         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detAll[k].first.rawId() << std::dec
1097                                      << detAll[k].second;
1098       edm::LogVerbatim("IsoTrack") << "Matched detIds" << detSum.size();
1099       for (unsigned int k = 0; k < detSum.size(); ++k)
1100         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detSum[k].first.rawId() << std::dec
1101                                      << detSum[k].second;
1102       edm::LogVerbatim("IsoTrack") << "Gamma detIds" << detGamma.size();
1103       for (unsigned int k = 0; k < detGamma.size(); ++k)
1104         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detGamma[k].first.rawId() << std::dec
1105                                      << detGamma[k].second;
1106       edm::LogVerbatim("IsoTrack") << "Charged detIds" << detSum.size();
1107       for (unsigned int k = 0; k < detCharged.size(); ++k)
1108         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detCharged[k].first.rawId() << std::dec
1109                                      << detCharged[k].second;
1110       edm::LogVerbatim("IsoTrack") << "Neutral detIds" << detNeutral.size();
1111       for (unsigned int k = 0; k < detNeutral.size(); ++k)
1112         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detNeutral[k].first.rawId() << std::dec
1113                                      << detNeutral[k].second;
1114       edm::LogVerbatim("IsoTrack") << "Rest detIds" << detRest.size();
1115       for (unsigned int k = 0; k < detRest.size(); ++k)
1116         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detRest[k].first.rawId() << std::dec
1117                                      << detRest[k].second;
1118     }
1119     return simInfo;
1120   }
1121 
1122   template <typename T>
1123   std::vector<typename T::const_iterator> missedECALHits(const edm::Event& iEvent,
1124                                                          CaloNavigator<DetId>& navigator,
1125                                                          edm::Handle<T>& hits,
1126                                                          edm::Handle<edm::SimTrackContainer>& SimTk,
1127                                                          edm::Handle<edm::SimVertexContainer>& SimVtx,
1128                                                          const reco::Track* pTrack,
1129                                                          TrackerHitAssociator& associate,
1130                                                          int ieta,
1131                                                          int iphi,
1132                                                          bool flag,
1133                                                          bool debug) {
1134     if (debug) {
1135       edm::LogVerbatim("IsoTrack") << "Processing missedECALHits " << 2 * ieta + 1 << "x" << 2 * iphi + 1
1136                                    << " and Flag " << flag;
1137       edm::LogVerbatim("IsoTrack") << "trkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
1138                                    << pTrack->recHitsSize();
1139     }
1140 
1141     std::vector<typename T::const_iterator> ecalHits;
1142     spr::hitECALmatrix(navigator, hits, ieta, iphi, ecalHits, debug);
1143 
1144     std::vector<int> matchedId = spr::matchedSimTrackId(iEvent, SimTk, SimVtx, pTrack, associate, debug);
1145 
1146     return spr::missedCaloHits(hits, matchedId, ecalHits, flag, false, debug);
1147   }
1148 
1149   template <typename T>
1150   std::vector<typename T::const_iterator> missedHCALHits(const edm::Event& iEvent,
1151                                                          const HcalTopology* topology,
1152                                                          const DetId& det,
1153                                                          edm::Handle<T>& hits,
1154                                                          edm::Handle<edm::SimTrackContainer>& SimTk,
1155                                                          edm::Handle<edm::SimVertexContainer>& SimVtx,
1156                                                          const reco::Track* pTrack,
1157                                                          TrackerHitAssociator& associate,
1158                                                          int ieta,
1159                                                          int iphi,
1160                                                          bool flag,
1161                                                          bool includeHO,
1162                                                          bool debug) {
1163     if (debug) {
1164       edm::LogVerbatim("IsoTrack") << "Processing missedHCALHits " << 2 * ieta + 1 << "x" << 2 * iphi + 1
1165                                    << " and Flag " << flag;
1166       edm::LogVerbatim("IsoTrack") << "trkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
1167                                    << pTrack->recHitsSize();
1168     }
1169 
1170     std::vector<typename T::const_iterator> hit;
1171     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
1172 
1173     std::vector<int> matchedId = spr::matchedSimTrackId(iEvent, SimTk, SimVtx, pTrack, associate, debug);
1174 
1175     return spr::missedCaloHits(hits, matchedId, hit, flag, includeHO, debug);
1176   }
1177 
1178   template <typename T>
1179   std::vector<typename T::const_iterator> missedCaloHits(edm::Handle<T>& hits,
1180                                                          std::vector<int> matchedId,
1181                                                          std::vector<typename T::const_iterator> caloHits,
1182                                                          bool flag,
1183                                                          bool includeHO,
1184                                                          bool debug) {
1185     std::vector<typename T::const_iterator> missedHits;
1186     std::vector<typename T::const_iterator> usedHits;
1187     if (!matchedId.empty()) {
1188       typename T::const_iterator ihit;
1189       for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
1190         int id = ihit->geantTrackId();
1191         bool found = false;
1192         for (unsigned int it = 0; it < matchedId.size(); ++it) {
1193           if (id == matchedId[it]) {
1194             found = true;
1195             break;
1196           }
1197         }
1198         DetId id_ = (DetId)(ihit->id());
1199         bool ok = true;
1200         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
1201           ok = false;
1202         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
1203           ok = false;
1204         if (found && ok) {
1205           if (flag) {
1206             if (count(caloHits.begin(), caloHits.end(), ihit) == 0)
1207               missedHits.push_back(ihit);
1208           } else {
1209             usedHits.push_back(ihit);
1210           }
1211         }
1212       }
1213     }
1214     if (!flag) {
1215       for (unsigned int ii = 0; ii < caloHits.size(); ii++) {
1216         DetId id_ = (DetId)(caloHits[ii]->id());
1217         bool ok = true;
1218         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
1219           ok = false;
1220         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
1221           ok = false;
1222         if (count(usedHits.begin(), usedHits.end(), caloHits[ii]) == 0 && ok)
1223           missedHits.push_back(caloHits[ii]);
1224       }
1225     }
1226     if (debug) {
1227       edm::LogVerbatim("IsoTrack") << "missedCaloHits finds " << missedHits.size() << " missed hits";
1228       for (unsigned int i = 0; i < missedHits.size(); i++)
1229         edm::LogVerbatim("IsoTrack") << "Hit " << i << " " << *missedHits[i];
1230     }
1231     return missedHits;
1232   }
1233 
1234 }  // namespace spr
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!