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
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
0070 edm::SimTrackContainer::const_iterator trkInfo =
0071 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0072
0073
0074 std::vector<DetId> vdets = spr::matrixECALIds(det, ieta, iphi, geo, caloTopology, debug);
0075
0076
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
0104 edm::SimTrackContainer::const_iterator trkInfo =
0105 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0106
0107
0108 std::vector<DetId> vdets = spr::matrixECALIds(det, ietaE, ietaW, iphiN, iphiS, geo, caloTopology, debug);
0109
0110
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
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
0158 edm::SimTrackContainer::const_iterator trkInfo =
0159 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0160
0161
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
0192 edm::SimTrackContainer::const_iterator trkInfo =
0193 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0194
0195
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
0203
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
0224 edm::SimTrackContainer::const_iterator trkInfo =
0225 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0226
0227
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
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
0251 edm::SimTrackContainer::const_iterator trkInfo =
0252 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug_dummy);
0253
0254
0255
0256
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
0263
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
0281 edm::SimTrackContainer::const_iterator trkInfo =
0282 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug_dummy);
0283
0284
0285
0286
0287 std::vector<typename T::const_iterator> hit = spr::findHitCone(geo, hits, hpoint1, point1, dR, trackMom);
0288
0289
0290
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;
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();
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
0391 if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0392 info.eMatched += hit[ihit]->energy();
0393 } else {
0394
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 }
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
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
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
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
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;
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
0594
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
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
0611 if (trkInfo != SimTk->end()) {
0612 unsigned int matchSimTrk = trkInfo->trackId();
0613 matchedId = trkInfo->type();
0614
0615
0616 for (unsigned int ihit = 0; ihit < hit.size(); ihit++) {
0617
0618
0619
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 }
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 }
0692 }
0693 }
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 }
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
0738 edm::SimTrackContainer::const_iterator trkInfo =
0739 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0740
0741 std::vector<typename T::const_iterator> ecalHits;
0742 spr::hitECALmatrix(navigator, hits, ieta, iphi, ecalHits, debug);
0743
0744
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
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
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
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
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
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
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
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
0927 edm::SimTrackContainer::const_iterator trkInfo =
0928 spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
0929
0930
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;
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();
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
0977 if (hit[ihit]->geantTrackId() == (int)matchSimTrk) {
0978 energySum = hit[ihit]->energy();
0979 } else {
0980
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 }
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 }