Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:35:01 

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     double energySum = 0.0;
0606     double energyGamma = 0.0, energyNeutral = 0.0, energyCharged = 0.0, energyRest = 0.0;
0607 
0608     // What is this checking?
0609     if (trkInfo != SimTk->end()) {
0610       unsigned int matchSimTrk = trkInfo->trackId();
0611       matchedId = trkInfo->type();  //pdgid
0612 
0613       // Loop over hits
0614       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0615         // For simhit tower multiplicity, we need to avoid double
0616         // counting towers with multiple simhits:
0617         // This will break if called for pcalohits from the ecal!!
0618         HcalDetId detId(hit[ihit]->id());
0619         int ieta = detId.ieta();
0620         int iphi = detId.iphi();
0621         int uniqueId = 100 * ieta + iphi;
0622 
0623         if (hit[ihit]->time() > 150.0)
0624           continue;
0625 
0626         if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0627           energySum += spr::getEnergy(hit[ihit]);
0628           uniqueIds_matched.insert(uniqueId);
0629           uniqueIds_total.insert(uniqueId);
0630         } else {
0631           bool found = false;
0632           for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0633             if (hit[ihit]->geantTrackId() == (int)simTrkItr->trackId()) {
0634               found = true;
0635               bool match = spr::validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, debug);
0636               if (match) {
0637                 energySum += spr::getEnergy(hit[ihit]);
0638                 uniqueIds_matched.insert(uniqueId);
0639                 uniqueIds_total.insert(uniqueId);
0640               } else {
0641                 edm::SimTrackContainer::const_iterator parentItr = spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0642                 if (debug) {
0643                   if (parentItr != SimTk->end()) {
0644                     edm::LogVerbatim("IsoTrack")
0645                         << "original parent of " << simTrkItr->trackId() << " " << parentItr->trackId() << ", "
0646                         << parentItr->type() << " Energy " << spr::getEnergy(hit[ihit]);
0647                   } else {
0648                     edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId()
0649                                                  << " not found; Energy " << spr::getEnergy(hit[ihit]);
0650                   }
0651                 }
0652                 if (parentItr == SimTk->end()) {
0653                   energyRest += spr::getEnergy(hit[ihit]);
0654                   uniqueIds_rest.insert(uniqueId);
0655                   uniqueIds_total.insert(uniqueId);
0656                 } else if (parentItr->type() == 22) {
0657                   energyGamma += spr::getEnergy(hit[ihit]);
0658                   uniqueIds_gamma.insert(uniqueId);
0659                   uniqueIds_total.insert(uniqueId);
0660                 } else if ((int)parentItr->charge() == 0) {
0661                   energyNeutral += spr::getEnergy(hit[ihit]);
0662                   uniqueIds_neut.insert(uniqueId);
0663                   uniqueIds_total.insert(uniqueId);
0664                 } else {
0665                   energyCharged += spr::getEnergy(hit[ihit]);
0666                   uniqueIds_char.insert(uniqueId);
0667                   uniqueIds_total.insert(uniqueId);
0668                 }
0669               }
0670               break;
0671             }
0672           }  // for (simTrkItr = ...)
0673 
0674           if (!found) {
0675             energyRest += spr::getEnergy(hit[ihit]);
0676             uniqueIds_rest.insert(uniqueId);
0677             uniqueIds_total.insert(uniqueId);
0678           }
0679 
0680           if (debug)
0681             edm::LogVerbatim("IsoTrack") << "Hit " << ihit << ": " << *hit[ihit];
0682         }  // else condition, i.e. (hit[ihit]->geantTrackId() != (int)matchSimTrk )
0683       }  // loop over hits
0684     }  // if (trkInfo != SimTk->end())
0685 
0686     double energyTot = energySum + energyGamma + energyNeutral + energyCharged + energyRest;
0687     multiplicityVector.push_back(uniqueIds_matched.size());
0688     multiplicityVector.push_back(uniqueIds_total.size());
0689     multiplicityVector.push_back(uniqueIds_neut.size());
0690     multiplicityVector.push_back(uniqueIds_char.size());
0691     multiplicityVector.push_back(uniqueIds_gamma.size());
0692     multiplicityVector.push_back(uniqueIds_rest.size());
0693 
0694     std::map<std::string, double> simInfo;
0695     simInfo.insert(std::pair<std::string, double>("eMatched", energySum));
0696     simInfo.insert(std::pair<std::string, double>("pdgMatched", matchedId));
0697     simInfo.insert(std::pair<std::string, double>("eGamma", energyGamma));
0698     simInfo.insert(std::pair<std::string, double>("eNeutralHad", energyNeutral));
0699     simInfo.insert(std::pair<std::string, double>("eChargedHad", energyCharged));
0700     simInfo.insert(std::pair<std::string, double>("eRest", energyRest));
0701     simInfo.insert(std::pair<std::string, double>("eTotal", energyTot));
0702     if (debug) {
0703       edm::LogVerbatim("IsoTrack") << " energySum " << energySum << "  energyGamma " << energyGamma
0704                                    << "  energyNeutral " << energyNeutral << "  energyCharged " << energyCharged
0705                                    << "  energyRest " << energyRest << "  energyTot " << energyTot;
0706     }
0707     return simInfo;
0708   }  // eCaloSimInfo
0709 
0710   template <typename T>
0711   std::map<std::string, double> eECALSimInfo(const edm::Event& iEvent,
0712                                              CaloNavigator<DetId>& navigator,
0713                                              const CaloGeometry* geo,
0714                                              edm::Handle<T>& hits,
0715                                              edm::Handle<edm::SimTrackContainer>& SimTk,
0716                                              edm::Handle<edm::SimVertexContainer>& SimVtx,
0717                                              const reco::Track* pTrack,
0718                                              TrackerHitAssociator& associate,
0719                                              int ieta,
0720                                              int iphi,
0721                                              double timeCut,
0722                                              bool debug) {
0723     if (debug)
0724       edm::LogVerbatim("IsoTrack") << "Processing eECALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0725                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0726                                    << pTrack->recHitsSize();
0727 
0728     //matching SimTrack
0729     edm::SimTrackContainer::const_iterator trkInfo =
0730         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0731     //vector of Ecal hits in NxN
0732     std::vector<typename T::const_iterator> ecalHits;
0733     spr::hitECALmatrix(navigator, hits, ieta, iphi, ecalHits, debug);
0734 
0735     // return a map of matching type and energy of SimHits
0736     spr::caloSimInfo info;
0737     spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, ecalHits, trkInfo, info, timeCut, false, debug);
0738     return spr::eCaloSimInfo(info);
0739   }
0740 
0741   template <typename T>
0742   std::map<std::string, double> eECALSimInfoTotal(const edm::Event& iEvent,
0743                                                   const DetId& det,
0744                                                   const CaloGeometry* geo,
0745                                                   const CaloTopology* caloTopology,
0746                                                   edm::Handle<T>& hitsEB,
0747                                                   edm::Handle<T>& hitsEE,
0748                                                   edm::Handle<edm::SimTrackContainer>& SimTk,
0749                                                   edm::Handle<edm::SimVertexContainer>& SimVtx,
0750                                                   const reco::Track* pTrack,
0751                                                   TrackerHitAssociator& associate,
0752                                                   int ieta,
0753                                                   int iphi,
0754                                                   int itry,
0755                                                   double timeCut,
0756                                                   bool debug) {
0757     if (debug)
0758       edm::LogVerbatim("IsoTrack") << "Processing eECALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0759                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0760                                    << pTrack->recHitsSize();
0761 
0762     //matching SimTrack
0763     edm::SimTrackContainer::const_iterator trkInfo =
0764         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0765 
0766     spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);
0767     std::map<std::string, double> simInfo;
0768 
0769     if (itry >= 0) {
0770       //vector of Ecal cells in NxN
0771       std::vector<DetId> vdets = spr::matrixECALIds(
0772           det, etaphi.ietaE[itry], etaphi.ietaW[itry], etaphi.iphiN[itry], etaphi.iphiS[itry], geo, caloTopology, debug);
0773       // get a map of matching type and energy of SimHits
0774       spr::caloSimInfo info;
0775       spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, info, timeCut, debug);
0776       simInfo = spr::eCaloSimInfo(info);
0777     } else {
0778       int itrym = 0;
0779       std::vector<DetId> vdets = spr::matrixECALIds(det,
0780                                                     etaphi.ietaE[itrym],
0781                                                     etaphi.ietaW[itrym],
0782                                                     etaphi.iphiN[itrym],
0783                                                     etaphi.iphiS[itrym],
0784                                                     geo,
0785                                                     caloTopology,
0786                                                     debug);
0787       // get a map of matching type and energy of SimHits
0788       spr::caloSimInfo info;
0789       spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, info, timeCut, debug);
0790       simInfo = spr::eCaloSimInfo(info);
0791       for (int itrys = 1; itrys < etaphi.ntrys; ++itrys) {
0792         vdets = spr::matrixECALIds(det,
0793                                    etaphi.ietaE[itrys],
0794                                    etaphi.ietaW[itrys],
0795                                    etaphi.iphiN[itrys],
0796                                    etaphi.iphiS[itrys],
0797                                    geo,
0798                                    caloTopology,
0799                                    debug);
0800         // get a map of matching type and energy of SimHits
0801         spr::caloSimInfo infox;
0802         spr::eCaloSimInfo(vdets, geo, hitsEB, hitsEE, SimTk, SimVtx, trkInfo, infox, timeCut, debug);
0803         std::map<std::string, double> simInfoX = spr::eCaloSimInfo(infox);
0804         if (simInfoX["eTotal"] > simInfo["eTotal"]) {
0805           simInfo = simInfoX;
0806           itrym = itrys;
0807         }
0808       }
0809     }
0810     return simInfo;
0811   }
0812 
0813   template <typename T>
0814   std::map<std::string, double> eHCALSimInfoTotal(const edm::Event& iEvent,
0815                                                   const HcalTopology* topology,
0816                                                   const DetId& det,
0817                                                   const CaloGeometry* geo,
0818                                                   edm::Handle<T>& hits,
0819                                                   edm::Handle<edm::SimTrackContainer>& SimTk,
0820                                                   edm::Handle<edm::SimVertexContainer>& SimVtx,
0821                                                   const reco::Track* pTrack,
0822                                                   TrackerHitAssociator& associate,
0823                                                   int ieta,
0824                                                   int iphi,
0825                                                   int itry,
0826                                                   double timeCut,
0827                                                   bool includeHO,
0828                                                   bool debug) {
0829     if (debug)
0830       edm::LogVerbatim("IsoTrack") << "Processing eHCALSimInfo " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << "\ntrkMom "
0831                                    << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0832                                    << pTrack->recHitsSize();
0833 
0834     // get the matching SimTrack pointer
0835     edm::SimTrackContainer::const_iterator trkInfo =
0836         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0837 
0838     spr::EtaPhi etaphi = spr::getEtaPhi(ieta, iphi, debug);
0839     std::map<std::string, double> simInfo;
0840 
0841     if (itry >= 0) {
0842       // get the hits in Hcal in NxN around det
0843       std::vector<typename T::const_iterator> hit;
0844       spr::hitHCALmatrixTotal(topology,
0845                               det,
0846                               hits,
0847                               etaphi.ietaE[itry],
0848                               etaphi.ietaW[itry],
0849                               etaphi.iphiN[itry],
0850                               etaphi.iphiS[itry],
0851                               hit,
0852                               false,
0853                               debug);
0854       spr::caloSimInfo info;
0855       spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, info, timeCut, includeHO, debug);
0856       simInfo = spr::eCaloSimInfo(info);
0857     } else {
0858       int itrym = 0;
0859       std::vector<typename T::const_iterator> hit;
0860       spr::hitHCALmatrixTotal(topology,
0861                               det,
0862                               hits,
0863                               etaphi.ietaE[itrym],
0864                               etaphi.ietaW[itrym],
0865                               etaphi.iphiN[itrym],
0866                               etaphi.iphiS[itrym],
0867                               hit,
0868                               includeHO,
0869                               debug);
0870       spr::caloSimInfo info;
0871       spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, info, timeCut, includeHO, debug);
0872       simInfo = spr::eCaloSimInfo(info);
0873       for (int itrys = 1; itrys < etaphi.ntrys; ++itrys) {
0874         hit.clear();
0875         spr::hitHCALmatrixTotal(topology,
0876                                 det,
0877                                 hits,
0878                                 etaphi.ietaE[itrys],
0879                                 etaphi.ietaW[itrys],
0880                                 etaphi.iphiN[itrys],
0881                                 etaphi.iphiS[itrys],
0882                                 hit,
0883                                 includeHO,
0884                                 debug);
0885         spr::caloSimInfo infox;
0886         spr::eCaloSimInfo(geo, hits, SimTk, SimVtx, hit, trkInfo, infox, timeCut, includeHO, debug);
0887         std::map<std::string, double> simInfoX = spr::eCaloSimInfo(infox);
0888         if (simInfoX["eTotal"] > simInfo["eTotal"]) {
0889           simInfo = simInfoX;
0890           itrym = itrys;
0891         }
0892       }
0893     }
0894     return simInfo;
0895   }
0896 
0897   template <typename T>
0898   spr::energyMap eHCALSimInfoMatrix(const edm::Event& iEvent,
0899                                     const HcalTopology* topology,
0900                                     const DetId& det,
0901                                     const CaloGeometry* geo,
0902                                     edm::Handle<T>& hits,
0903                                     edm::Handle<edm::SimTrackContainer>& SimTk,
0904                                     edm::Handle<edm::SimVertexContainer>& SimVtx,
0905                                     const reco::Track* pTrack,
0906                                     TrackerHitAssociator& associate,
0907                                     int ieta,
0908                                     int iphi,
0909                                     double timeCut,
0910                                     bool includeHO,
0911                                     bool debug) {
0912     if (debug)
0913       edm::LogVerbatim("IsoTrack") << "Processing eHCALSimInfoMatrix " << 2 * ieta + 1 << "x" << 2 * iphi + 1
0914                                    << "\ntrkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
0915                                    << pTrack->recHitsSize();
0916 
0917     // get the matching SimTrack pointer
0918     edm::SimTrackContainer::const_iterator trkInfo =
0919         spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0920 
0921     // get the hits in Hcal in NxN around det
0922     std::vector<typename T::const_iterator> hit;
0923     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
0924 
0925     return spr::caloSimInfoMatrix(geo, hits, SimTk, SimVtx, hit, trkInfo, timeCut, includeHO, debug);
0926   }
0927 
0928   template <typename T>
0929   spr::energyMap caloSimInfoMatrix(const CaloGeometry* geo,
0930                                    edm::Handle<T>& hits,
0931                                    edm::Handle<edm::SimTrackContainer>& SimTk,
0932                                    edm::Handle<edm::SimVertexContainer>& SimVtx,
0933                                    std::vector<typename T::const_iterator> hit,
0934                                    edm::SimTrackContainer::const_iterator trkInfo,
0935                                    double timeCut,
0936                                    bool includeHO,
0937                                    bool debug) {
0938     int matchedId = 0;  //pdgid
0939 
0940     if (debug) {
0941       if (trkInfo != SimTk->end())
0942         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  matchSimTrk:" << trkInfo->trackId() << " matchedId "
0943                                      << trkInfo->type();
0944       else
0945         edm::LogVerbatim("IsoTrack") << "In eCaloSimInfo::  not valid track pointer";
0946     }
0947 
0948     std::vector<std::pair<DetId, double> > detSum, detGamma, detCharged, detNeutral, detRest, detAll;
0949 
0950     if (trkInfo != SimTk->end()) {
0951       unsigned int matchSimTrk = trkInfo->trackId();
0952       matchedId = trkInfo->type();  //pdgid
0953 
0954       edm::SimTrackContainer::const_iterator simTrkItr;
0955 
0956       for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0957         DetId id_ = (DetId)(hit[ihit]->id());
0958         double tof = timeOfFlight(id_, geo, debug);
0959         double energySum = 0, energyRest = 0;
0960         double energyGamma = 0, energyNeutral = 0, energyCharged = 0;
0961         bool ok = true;
0962         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
0963           ok = false;
0964         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
0965           ok = false;
0966         if ((hit[ihit]->time() <= (tof + timeCut)) && ok) {
0967           // if the hitId matches with matching trackId
0968           if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0969             energySum = hit[ihit]->energy();
0970           } else {
0971             // trace back the history and check the pdgId of origin SimTrack of SimHit
0972             bool found = false;
0973             for (simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
0974               if (hit[ihit]->geantTrackId() == (int)simTrkItr->trackId()) {
0975                 found = true;
0976                 bool match = spr::validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, debug);
0977                 if (match) {
0978                   energySum = hit[ihit]->energy();
0979                 } else {
0980                   edm::SimTrackContainer::const_iterator parentItr =
0981                       spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
0982                   if (debug) {
0983                     if (parentItr != SimTk->end())
0984                       edm::LogVerbatim("IsoTrack")
0985                           << "original parent of " << simTrkItr->trackId() << " " << parentItr->trackId() << ", "
0986                           << parentItr->type() << " Energy " << hit[ihit]->energy();
0987                     else
0988                       edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId()
0989                                                    << " not found; Energy " << hit[ihit]->energy();
0990                   }
0991                   if (parentItr == SimTk->end()) {
0992                     energyRest = hit[ihit]->energy();
0993                   } else if (parentItr->type() == 22) {
0994                     energyGamma = hit[ihit]->energy();
0995                   } else if ((int)parentItr->charge() == 0) {
0996                     energyNeutral = hit[ihit]->energy();
0997                   } else
0998                     energyCharged = hit[ihit]->energy();
0999                 }
1000                 break;
1001               }
1002             }
1003             if (!found)
1004               energyRest = hit[ihit]->energy();
1005             if (debug)
1006               edm::LogVerbatim("IsoTrack") << "Hit " << ihit << ": " << *hit[ihit];
1007             found = false;
1008             if (energySum > 0) {
1009               for (unsigned int k = 0; k < detSum.size(); k++) {
1010                 if (id_ == detSum[k].first) {
1011                   found = true;
1012                   detSum[k].second += energySum;
1013                   break;
1014                 }
1015               }
1016               if (!found)
1017                 detSum.push_back(std::pair<DetId, double>(id_, energySum));
1018             } else if (energyRest > 0) {
1019               for (unsigned int k = 0; k < detRest.size(); k++) {
1020                 if (id_ == detRest[k].first) {
1021                   found = true;
1022                   detRest[k].second += energyRest;
1023                   break;
1024                 }
1025               }
1026               if (!found)
1027                 detRest.push_back(std::pair<DetId, double>(id_, energyRest));
1028             } else if (energyGamma > 0) {
1029               for (unsigned int k = 0; k < detGamma.size(); k++) {
1030                 if (id_ == detGamma[k].first) {
1031                   found = true;
1032                   detGamma[k].second += energyGamma;
1033                   break;
1034                 }
1035               }
1036               if (!found)
1037                 detGamma.push_back(std::pair<DetId, double>(id_, energyGamma));
1038             } else if (energyCharged > 0) {
1039               for (unsigned int k = 0; k < detCharged.size(); k++) {
1040                 if (id_ == detCharged[k].first) {
1041                   found = true;
1042                   detCharged[k].second += energyCharged;
1043                   break;
1044                 }
1045               }
1046               if (!found)
1047                 detCharged.push_back(std::pair<DetId, double>(id_, energyCharged));
1048             } else if (energyNeutral > 0) {
1049               for (unsigned int k = 0; k < detNeutral.size(); k++) {
1050                 if (id_ == detNeutral[k].first) {
1051                   found = true;
1052                   detNeutral[k].second += energyNeutral;
1053                   break;
1054                 }
1055               }
1056               if (!found)
1057                 detNeutral.push_back(std::pair<DetId, double>(id_, energyNeutral));
1058             }
1059             found = false;
1060             double energyTot = energySum + energyGamma + energyNeutral + energyCharged + energyRest;
1061             for (unsigned int k = 0; k < detAll.size(); k++) {
1062               if (id_ == detAll[k].first) {
1063                 found = true;
1064                 detAll[k].second += energyTot;
1065                 break;
1066               }
1067             }
1068             if (!found)
1069               detAll.push_back(std::pair<DetId, double>(id_, energyTot));
1070           }
1071         }
1072       }  // loop over hits
1073     }
1074 
1075     spr::energyMap simInfo;
1076     simInfo.pdgId = matchedId;
1077     simInfo.matched = detSum;
1078     simInfo.gamma = detGamma;
1079     simInfo.charged = detCharged;
1080     simInfo.neutral = detNeutral;
1081     simInfo.rest = detRest;
1082     simInfo.all = detAll;
1083     if (debug) {
1084       edm::LogVerbatim("IsoTrack") << "CaloSimInfo:: for particle " << simInfo.pdgId << "\n"
1085                                    << "All detIds " << detAll.size();
1086       for (unsigned int k = 0; k < detAll.size(); ++k)
1087         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detAll[k].first.rawId() << std::dec
1088                                      << detAll[k].second;
1089       edm::LogVerbatim("IsoTrack") << "Matched detIds" << detSum.size();
1090       for (unsigned int k = 0; k < detSum.size(); ++k)
1091         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detSum[k].first.rawId() << std::dec
1092                                      << detSum[k].second;
1093       edm::LogVerbatim("IsoTrack") << "Gamma detIds" << detGamma.size();
1094       for (unsigned int k = 0; k < detGamma.size(); ++k)
1095         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detGamma[k].first.rawId() << std::dec
1096                                      << detGamma[k].second;
1097       edm::LogVerbatim("IsoTrack") << "Charged detIds" << detSum.size();
1098       for (unsigned int k = 0; k < detCharged.size(); ++k)
1099         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detCharged[k].first.rawId() << std::dec
1100                                      << detCharged[k].second;
1101       edm::LogVerbatim("IsoTrack") << "Neutral detIds" << detNeutral.size();
1102       for (unsigned int k = 0; k < detNeutral.size(); ++k)
1103         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detNeutral[k].first.rawId() << std::dec
1104                                      << detNeutral[k].second;
1105       edm::LogVerbatim("IsoTrack") << "Rest detIds" << detRest.size();
1106       for (unsigned int k = 0; k < detRest.size(); ++k)
1107         edm::LogVerbatim("IsoTrack") << k << " detId 0x" << std::hex << detRest[k].first.rawId() << std::dec
1108                                      << detRest[k].second;
1109     }
1110     return simInfo;
1111   }
1112 
1113   template <typename T>
1114   std::vector<typename T::const_iterator> missedECALHits(const edm::Event& iEvent,
1115                                                          CaloNavigator<DetId>& navigator,
1116                                                          edm::Handle<T>& hits,
1117                                                          edm::Handle<edm::SimTrackContainer>& SimTk,
1118                                                          edm::Handle<edm::SimVertexContainer>& SimVtx,
1119                                                          const reco::Track* pTrack,
1120                                                          TrackerHitAssociator& associate,
1121                                                          int ieta,
1122                                                          int iphi,
1123                                                          bool flag,
1124                                                          bool debug) {
1125     if (debug) {
1126       edm::LogVerbatim("IsoTrack") << "Processing missedECALHits " << 2 * ieta + 1 << "x" << 2 * iphi + 1
1127                                    << " and Flag " << flag;
1128       edm::LogVerbatim("IsoTrack") << "trkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
1129                                    << pTrack->recHitsSize();
1130     }
1131 
1132     std::vector<typename T::const_iterator> ecalHits;
1133     spr::hitECALmatrix(navigator, hits, ieta, iphi, ecalHits, debug);
1134 
1135     std::vector<int> matchedId = spr::matchedSimTrackId(iEvent, SimTk, SimVtx, pTrack, associate, debug);
1136 
1137     return spr::missedCaloHits(hits, matchedId, ecalHits, flag, false, debug);
1138   }
1139 
1140   template <typename T>
1141   std::vector<typename T::const_iterator> missedHCALHits(const edm::Event& iEvent,
1142                                                          const HcalTopology* topology,
1143                                                          const DetId& det,
1144                                                          edm::Handle<T>& hits,
1145                                                          edm::Handle<edm::SimTrackContainer>& SimTk,
1146                                                          edm::Handle<edm::SimVertexContainer>& SimVtx,
1147                                                          const reco::Track* pTrack,
1148                                                          TrackerHitAssociator& associate,
1149                                                          int ieta,
1150                                                          int iphi,
1151                                                          bool flag,
1152                                                          bool includeHO,
1153                                                          bool debug) {
1154     if (debug) {
1155       edm::LogVerbatim("IsoTrack") << "Processing missedHCALHits " << 2 * ieta + 1 << "x" << 2 * iphi + 1
1156                                    << " and Flag " << flag;
1157       edm::LogVerbatim("IsoTrack") << "trkMom " << pTrack->p() << " eta " << pTrack->eta() << " trkRecHits "
1158                                    << pTrack->recHitsSize();
1159     }
1160 
1161     std::vector<typename T::const_iterator> hit;
1162     spr::hitHCALmatrix(topology, det, hits, ieta, iphi, hit, includeHO, debug);
1163 
1164     std::vector<int> matchedId = spr::matchedSimTrackId(iEvent, SimTk, SimVtx, pTrack, associate, debug);
1165 
1166     return spr::missedCaloHits(hits, matchedId, hit, flag, includeHO, debug);
1167   }
1168 
1169   template <typename T>
1170   std::vector<typename T::const_iterator> missedCaloHits(edm::Handle<T>& hits,
1171                                                          std::vector<int> matchedId,
1172                                                          std::vector<typename T::const_iterator> caloHits,
1173                                                          bool flag,
1174                                                          bool includeHO,
1175                                                          bool debug) {
1176     std::vector<typename T::const_iterator> missedHits;
1177     std::vector<typename T::const_iterator> usedHits;
1178     if (!matchedId.empty()) {
1179       typename T::const_iterator ihit;
1180       for (ihit = hits->begin(); ihit != hits->end(); ihit++) {
1181         int id = ihit->geantTrackId();
1182         bool found = false;
1183         for (unsigned int it = 0; it < matchedId.size(); ++it) {
1184           if (id == matchedId[it]) {
1185             found = true;
1186             break;
1187           }
1188         }
1189         DetId id_ = (DetId)(ihit->id());
1190         bool ok = true;
1191         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
1192           ok = false;
1193         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
1194           ok = false;
1195         if (found && ok) {
1196           if (flag) {
1197             if (count(caloHits.begin(), caloHits.end(), ihit) == 0)
1198               missedHits.push_back(ihit);
1199           } else {
1200             usedHits.push_back(ihit);
1201           }
1202         }
1203       }
1204     }
1205     if (!flag) {
1206       for (unsigned int ii = 0; ii < caloHits.size(); ii++) {
1207         DetId id_ = (DetId)(caloHits[ii]->id());
1208         bool ok = true;
1209         if (((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalForward)))
1210           ok = false;
1211         if ((!includeHO) && ((int)(id_.det()) == 4) && (id_.subdetId() == (int)(HcalOuter)))
1212           ok = false;
1213         if (count(usedHits.begin(), usedHits.end(), caloHits[ii]) == 0 && ok)
1214           missedHits.push_back(caloHits[ii]);
1215       }
1216     }
1217     if (debug) {
1218       edm::LogVerbatim("IsoTrack") << "missedCaloHits finds " << missedHits.size() << " missed hits";
1219       for (unsigned int i = 0; i < missedHits.size(); i++)
1220         edm::LogVerbatim("IsoTrack") << "Hit " << i << " " << *missedHits[i];
1221     }
1222     return missedHits;
1223   }
1224 
1225 }  // 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!