File indexing completed on 2023-03-17 10:43:41
0001 #include "Calibration/IsolatedParticles/interface/CaloConstants.h"
0002 #include "Calibration/IsolatedParticles/interface/ChargeIsolation.h"
0003 #include "Calibration/IsolatedParticles/interface/CaloPropagateTrack.h"
0004 #include "Calibration/IsolatedParticles/interface/FindDistCone.h"
0005 #include "Calibration/IsolatedParticles/interface/MatrixECALDetIds.h"
0006 #include "Calibration/IsolatedParticles/interface/MatrixHCALDetIds.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0009
0010 namespace spr {
0011
0012 double chargeIsolationEcal(unsigned int trkIndex,
0013 std::vector<spr::propagatedTrackID>& vdetIds,
0014 const CaloGeometry* geo,
0015 const CaloTopology* caloTopology,
0016 int ieta,
0017 int iphi,
0018 bool debug) {
0019 const DetId coreDet = vdetIds[trkIndex].detIdECAL;
0020 if (debug) {
0021 if (coreDet.subdetId() == EcalBarrel)
0022 edm::LogVerbatim("IsoTrack") << "DetId " << (EBDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
0023 else
0024 edm::LogVerbatim("IsoTrack") << "DetId " << (EEDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
0025 }
0026 double maxNearP = -1.0;
0027 if (vdetIds[trkIndex].okECAL) {
0028 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
0029 if (debug)
0030 edm::LogVerbatim("IsoTrack") << "chargeIsolationEcal:: eta/phi/dets " << ieta << " " << iphi << " "
0031 << vdets.size();
0032
0033 for (unsigned int indx = 0; indx < vdetIds.size(); ++indx) {
0034 if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okECAL) {
0035 const DetId anyCell = vdetIds[indx].detIdECAL;
0036 if (!spr::chargeIsolation(anyCell, vdets)) {
0037 const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
0038 if (maxNearP < pTrack->p())
0039 maxNearP = pTrack->p();
0040 }
0041 }
0042 }
0043 }
0044 return maxNearP;
0045 }
0046
0047 double chargeIsolationGenEcal(unsigned int trkIndex,
0048 std::vector<spr::propagatedGenParticleID>& vdetIds,
0049 const CaloGeometry* geo,
0050 const CaloTopology* caloTopology,
0051 int ieta,
0052 int iphi,
0053 bool debug) {
0054 const DetId coreDet = vdetIds[trkIndex].detIdECAL;
0055 if (debug) {
0056 if (coreDet.subdetId() == EcalBarrel)
0057 edm::LogVerbatim("IsoTrack") << "DetId " << (EBDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
0058 else
0059 edm::LogVerbatim("IsoTrack") << "DetId " << (EEDetId)(coreDet) << " Flag " << vdetIds[trkIndex].okECAL;
0060 }
0061 double maxNearP = -1.0;
0062 if (vdetIds[trkIndex].okECAL) {
0063 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
0064 if (debug)
0065 edm::LogVerbatim("IsoTrack") << "chargeIsolationGenEcal:: eta/phi/dets " << ieta << " " << iphi << " "
0066 << vdets.size();
0067
0068 for (unsigned int indx = 0; indx < vdetIds.size(); ++indx) {
0069 if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okECAL) {
0070 const DetId anyCell = vdetIds[indx].detIdECAL;
0071 if (!spr::chargeIsolation(anyCell, vdets)) {
0072 const reco::GenParticle* pTrack = &(*(vdetIds[indx].trkItr));
0073 if (maxNearP < pTrack->p())
0074 maxNearP = pTrack->p();
0075 }
0076 }
0077 }
0078 }
0079 return maxNearP;
0080 }
0081
0082 double chargeIsolationEcal(const DetId& coreDet,
0083 reco::TrackCollection::const_iterator trkItr,
0084 edm::Handle<reco::TrackCollection> trkCollection,
0085 const CaloGeometry* geo,
0086 const CaloTopology* caloTopology,
0087 const MagneticField* bField,
0088 int ieta,
0089 int iphi,
0090 const std::string& theTrackQuality,
0091 bool debug) {
0092 const CaloSubdetectorGeometry* barrelGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0093 const CaloSubdetectorGeometry* endcapGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
0094
0095 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
0096 if (debug)
0097 edm::LogVerbatim("IsoTrack") << "chargeIsolation:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size();
0098
0099 double maxNearP = -1.0;
0100 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0101
0102
0103 reco::TrackCollection::const_iterator trkItr2;
0104 for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0105 const reco::Track* pTrack2 = &(*trkItr2);
0106
0107 bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
0108 if ((trkItr2 != trkItr) && trkQuality) {
0109 std::pair<math::XYZPoint, bool> info = spr::propagateECAL(pTrack2, bField);
0110 const GlobalPoint point2(info.first.x(), info.first.y(), info.first.z());
0111
0112 if (info.second) {
0113 if (std::abs(point2.eta()) < spr::etaBEEcal) {
0114 const DetId anyCell = barrelGeom->getClosestCell(point2);
0115 if (!spr::chargeIsolation(anyCell, vdets)) {
0116 if (debug)
0117 edm::LogVerbatim("IsoTrack")
0118 << "chargeIsolationEcal Cell " << (EBDetId)(anyCell) << " pt " << pTrack2->p();
0119 if (maxNearP < pTrack2->p())
0120 maxNearP = pTrack2->p();
0121 }
0122 } else {
0123 if (endcapGeom) {
0124 const DetId anyCell = endcapGeom->getClosestCell(point2);
0125 if (!spr::chargeIsolation(anyCell, vdets)) {
0126 if (debug)
0127 edm::LogVerbatim("IsoTrack")
0128 << "chargeIsolationEcal Cell " << (EEDetId)(anyCell) << " pt " << pTrack2->p();
0129 if (maxNearP < pTrack2->p())
0130 maxNearP = pTrack2->p();
0131 }
0132 }
0133 }
0134 }
0135 }
0136 }
0137 return maxNearP;
0138 }
0139
0140 double chargeIsolationHcal(unsigned int trkIndex,
0141 std::vector<spr::propagatedTrackID>& vdetIds,
0142 const HcalTopology* topology,
0143 int ieta,
0144 int iphi,
0145 bool debug) {
0146 std::vector<DetId> dets(1, vdetIds[trkIndex].detIdHCAL);
0147 if (debug) {
0148 edm::LogVerbatim("IsoTrack") << "DetId " << (HcalDetId)(dets[0]) << " Flag " << vdetIds[trkIndex].okHCAL;
0149 }
0150
0151 double maxNearP = -1.0;
0152 if (vdetIds[trkIndex].okHCAL) {
0153 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
0154 if (debug)
0155 edm::LogVerbatim("IsoTrack") << "chargeIsolationHcal:: eta/phi/dets " << ieta << " " << iphi << " "
0156 << vdets.size();
0157
0158 for (unsigned indx = 0; indx < vdetIds.size(); ++indx) {
0159 if (indx != trkIndex && vdetIds[indx].ok && vdetIds[indx].okHCAL) {
0160 const DetId anyCell = vdetIds[indx].detIdHCAL;
0161 if (!spr::chargeIsolation(anyCell, vdets)) {
0162 const reco::Track* pTrack = &(*(vdetIds[indx].trkItr));
0163 if (debug)
0164 edm::LogVerbatim("IsoTrack")
0165 << "chargeIsolationHcal Cell " << (HcalDetId)(anyCell) << " pt " << pTrack->p();
0166 if (maxNearP < pTrack->p())
0167 maxNearP = pTrack->p();
0168 }
0169 }
0170 }
0171 }
0172 return maxNearP;
0173 }
0174
0175
0176
0177 double chargeIsolationHcal(reco::TrackCollection::const_iterator trkItr,
0178 edm::Handle<reco::TrackCollection> trkCollection,
0179 const DetId ClosestCell,
0180 const HcalTopology* topology,
0181 const CaloSubdetectorGeometry* gHB,
0182 const MagneticField* bField,
0183 int ieta,
0184 int iphi,
0185 const std::string& theTrackQuality,
0186 bool debug) {
0187 std::vector<DetId> dets(1, ClosestCell);
0188 if (debug)
0189 edm::LogVerbatim("IsoTrack") << (HcalDetId)ClosestCell;
0190 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
0191
0192 if (debug) {
0193 for (unsigned int i = 0; i < vdets.size(); i++) {
0194 edm::LogVerbatim("IsoTrack") << "HcalDetId in " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << " "
0195 << static_cast<HcalDetId>(vdets[i]);
0196 }
0197 }
0198
0199 double maxNearP = -1.0;
0200 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0201
0202 reco::TrackCollection::const_iterator trkItr2;
0203 for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0204 const reco::Track* pTrack2 = &(*trkItr2);
0205
0206 bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
0207 if ((trkItr2 != trkItr) && trkQuality) {
0208 std::pair<math::XYZPoint, bool> info = spr::propagateHCAL(pTrack2, bField);
0209 const GlobalPoint point2(info.first.x(), info.first.y(), info.first.z());
0210
0211 if (debug)
0212 edm::LogVerbatim("IsoTrack") << "Track2 (p,eta,phi) " << pTrack2->p() << " " << pTrack2->eta() << " "
0213 << pTrack2->phi();
0214 if (info.second) {
0215 const DetId anyCell = gHB->getClosestCell(point2);
0216 if (!spr::chargeIsolation(anyCell, vdets)) {
0217 if (maxNearP < pTrack2->p())
0218 maxNearP = pTrack2->p();
0219 }
0220
0221 if (debug)
0222 edm::LogVerbatim("IsoTrack") << "maxNearP " << maxNearP << " thisCell " << static_cast<HcalDetId>(anyCell)
0223 << " (" << info.first.x() << "," << info.first.y() << "," << info.first.z()
0224 << ")";
0225 }
0226 }
0227 }
0228 return maxNearP;
0229 }
0230
0231 bool chargeIsolation(const DetId anyCell, std::vector<DetId>& vdets) {
0232 bool isIsolated = true;
0233 for (unsigned int i = 0; i < vdets.size(); i++) {
0234 if (anyCell == vdets[i]) {
0235 isIsolated = false;
0236 break;
0237 }
0238 }
0239 return isIsolated;
0240 }
0241
0242 double coneChargeIsolation(const edm::Event& iEvent,
0243 const edm::EventSetup& iSetup,
0244 reco::TrackCollection::const_iterator trkItr,
0245 edm::Handle<reco::TrackCollection> trkCollection,
0246 TrackDetectorAssociator& associator,
0247 TrackAssociatorParameters& parameters_,
0248 const std::string& theTrackQuality,
0249 int& nNearTRKs,
0250 int& nLayers_maxNearP,
0251 int& trkQual_maxNearP,
0252 double& maxNearP_goodTrk,
0253 const GlobalPoint& hpoint1,
0254 const GlobalVector& trackMom,
0255 double dR) {
0256 nNearTRKs = 0;
0257 nLayers_maxNearP = 0;
0258 trkQual_maxNearP = -1;
0259 maxNearP_goodTrk = -999.0;
0260 double maxNearP = -999.0;
0261 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0262
0263
0264 reco::TrackCollection::const_iterator trkItr2;
0265 for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0266
0267 const reco::Track* pTrack2 = &(*trkItr2);
0268
0269
0270 bool trkQuality = (trackQuality_ != reco::TrackBase::undefQuality) ? (pTrack2->quality(trackQuality_)) : true;
0271 if (trkQuality)
0272 trkQual_maxNearP = 1;
0273 const reco::HitPattern& hitp = pTrack2->hitPattern();
0274 nLayers_maxNearP = hitp.trackerLayersWithMeasurement();
0275
0276
0277
0278 if (trkItr2 != trkItr) {
0279
0280 const FreeTrajectoryState fts2 =
0281 associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
0282 TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
0283
0284
0285 if (info2.isGoodHcal) {
0286 const GlobalPoint point2(info2.trkGlobPosAtHcal.x(), info2.trkGlobPosAtHcal.y(), info2.trkGlobPosAtHcal.z());
0287
0288 int isConeChargedIso = spr::coneChargeIsolation(hpoint1, point2, trackMom, dR);
0289
0290 if (isConeChargedIso == 0) {
0291 nNearTRKs++;
0292 if (maxNearP < pTrack2->p()) {
0293 maxNearP = pTrack2->p();
0294 if (trkQual_maxNearP > 0 && nLayers_maxNearP > 7 && maxNearP_goodTrk < pTrack2->p()) {
0295 maxNearP_goodTrk = pTrack2->p();
0296 }
0297 }
0298 }
0299 }
0300 }
0301 }
0302
0303 return maxNearP;
0304 }
0305
0306 double chargeIsolationCone(unsigned int trkIndex,
0307 std::vector<spr::propagatedTrackDirection>& trkDirs,
0308 double dR,
0309 int& nNearTRKs,
0310 bool debug) {
0311 double maxNearP = -1.0;
0312 nNearTRKs = 0;
0313 if (trkDirs[trkIndex].okHCAL) {
0314 if (debug)
0315 edm::LogVerbatim("IsoTrack") << "chargeIsolationCone with " << trkDirs.size() << " tracks ";
0316 for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
0317 if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
0318 int isConeChargedIso = spr::coneChargeIsolation(
0319 trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
0320 if (isConeChargedIso == 0) {
0321 nNearTRKs++;
0322 const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
0323 if (maxNearP < pTrack->p())
0324 maxNearP = pTrack->p();
0325 }
0326 }
0327 }
0328 }
0329
0330 if (debug)
0331 edm::LogVerbatim("IsoTrack") << "chargeIsolationCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP "
0332 << maxNearP;
0333 return maxNearP;
0334 }
0335
0336 double chargeIsolationGenCone(unsigned int trkIndex,
0337 std::vector<spr::propagatedGenParticleID>& trkDirs,
0338 double dR,
0339 int& nNearTRKs,
0340 bool debug) {
0341 double maxNearP = -1.0;
0342 nNearTRKs = 0;
0343 if (trkDirs[trkIndex].okHCAL) {
0344 if (debug)
0345 edm::LogVerbatim("IsoTrack") << "chargeIsolationCone with " << trkDirs.size() << " tracks ";
0346 for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
0347 if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
0348 int isConeChargedIso = spr::coneChargeIsolation(
0349 trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
0350 if (isConeChargedIso == 0) {
0351 nNearTRKs++;
0352 const reco::GenParticle* pTrack = &(*(trkDirs[indx].trkItr));
0353 if (maxNearP < pTrack->p())
0354 maxNearP = pTrack->p();
0355 }
0356 }
0357 }
0358 }
0359
0360 if (debug)
0361 edm::LogVerbatim("IsoTrack") << "chargeIsolationGenCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP "
0362 << maxNearP;
0363 return maxNearP;
0364 }
0365
0366 std::pair<double, double> chargeIsolationCone(unsigned int trkIndex,
0367 std::vector<spr::propagatedTrackDirection>& trkDirs,
0368 double dR,
0369 bool debug) {
0370 double maxNearP = -1.0;
0371 double sumP = 0;
0372 if (trkDirs[trkIndex].okHCAL) {
0373 if (debug)
0374 edm::LogVerbatim("IsoTrack") << "chargeIsolationCone with " << trkDirs.size() << " tracks ";
0375 for (unsigned int indx = 0; indx < trkDirs.size(); ++indx) {
0376 if (indx != trkIndex && trkDirs[indx].ok && trkDirs[indx].okHCAL) {
0377 int isConeChargedIso = spr::coneChargeIsolation(
0378 trkDirs[trkIndex].pointHCAL, trkDirs[indx].pointHCAL, trkDirs[trkIndex].directionHCAL, dR);
0379 if (isConeChargedIso == 0) {
0380 const reco::Track* pTrack = &(*(trkDirs[indx].trkItr));
0381 sumP += (pTrack->p());
0382 if (maxNearP < pTrack->p())
0383 maxNearP = pTrack->p();
0384 }
0385 }
0386 }
0387 }
0388
0389 if (debug)
0390 edm::LogVerbatim("IsoTrack") << "chargeIsolationCone Track " << trkDirs[trkIndex].okHCAL << " maxNearP "
0391 << maxNearP << ":" << sumP;
0392 return std::pair<double, double>(maxNearP, sumP);
0393 }
0394
0395 int coneChargeIsolation(const GlobalPoint& hpoint1,
0396 const GlobalPoint& point2,
0397 const GlobalVector& trackMom,
0398 double dR) {
0399 int isIsolated = 1;
0400 if (spr::getDistInPlaneTrackDir(hpoint1, trackMom, point2) > dR)
0401 isIsolated = 1;
0402 else
0403 isIsolated = 0;
0404 return isIsolated;
0405 }
0406
0407 double chargeIsolation(const edm::Event& iEvent,
0408 const edm::EventSetup& iSetup,
0409 CaloNavigator<DetId>& theNavigator,
0410 reco::TrackCollection::const_iterator trkItr,
0411 edm::Handle<reco::TrackCollection> trkCollection,
0412 const CaloSubdetectorGeometry* gEB,
0413 const CaloSubdetectorGeometry* gEE,
0414 TrackDetectorAssociator& associator,
0415 TrackAssociatorParameters& parameters_,
0416 int ieta,
0417 int iphi,
0418 const std::string& theTrackQuality,
0419 bool debug) {
0420 double maxNearP = -1.0;
0421 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0422
0423
0424 reco::TrackCollection::const_iterator trkItr2;
0425 for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0426 const reco::Track* pTrack2 = &(*trkItr2);
0427
0428 bool trkQuality = pTrack2->quality(trackQuality_);
0429 if ((trkItr2 != trkItr) && trkQuality) {
0430 const FreeTrajectoryState fts2 =
0431 associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
0432 TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
0433 const GlobalPoint point2(info2.trkGlobPosAtEcal.x(), info2.trkGlobPosAtEcal.y(), info2.trkGlobPosAtEcal.z());
0434
0435 if (info2.isGoodEcal) {
0436 if (std::abs(point2.eta()) < spr::etaBEEcal) {
0437 const DetId anyCell = gEB->getClosestCell(point2);
0438 if (debug)
0439 edm::LogVerbatim("IsoTrack")
0440 << "chargeIsolation:: EB cell " << (EBDetId)(anyCell) << " for pt " << pTrack2->p();
0441 if (!spr::chargeIsolation(anyCell, theNavigator, ieta, iphi)) {
0442 if (maxNearP < pTrack2->p())
0443 maxNearP = pTrack2->p();
0444 }
0445 } else {
0446 const DetId anyCell = gEE->getClosestCell(point2);
0447 if (debug)
0448 edm::LogVerbatim("IsoTrack")
0449 << "chargeIsolation:: EE cell " << (EEDetId)(anyCell) << " for pt " << pTrack2->p();
0450 if (!spr::chargeIsolation(anyCell, theNavigator, ieta, iphi)) {
0451 if (maxNearP < pTrack2->p())
0452 maxNearP = pTrack2->p();
0453 }
0454 }
0455 }
0456 }
0457 }
0458 return maxNearP;
0459 }
0460
0461
0462
0463 bool chargeIsolation(const DetId anyCell, CaloNavigator<DetId>& navigator, int ieta, int iphi) {
0464 bool isIsolated = false;
0465
0466 DetId thisDet;
0467
0468 for (int dx = -ieta; dx < ieta + 1; ++dx) {
0469 for (int dy = -iphi; dy < iphi + 1; ++dy) {
0470 thisDet = navigator.offsetBy(dx, dy);
0471 navigator.home();
0472
0473 if (thisDet != DetId(0)) {
0474 if (thisDet == anyCell) {
0475 isIsolated = false;
0476 return isIsolated;
0477 }
0478 }
0479 }
0480 }
0481 return isIsolated;
0482 }
0483
0484
0485
0486 double chargeIsolationEcal(const edm::Event& iEvent,
0487 const edm::EventSetup& iSetup,
0488 const DetId& coreDet,
0489 reco::TrackCollection::const_iterator trkItr,
0490 edm::Handle<reco::TrackCollection> trkCollection,
0491 const CaloGeometry* geo,
0492 const CaloTopology* caloTopology,
0493 TrackDetectorAssociator& associator,
0494 TrackAssociatorParameters& parameters_,
0495 int ieta,
0496 int iphi,
0497 const std::string& theTrackQuality,
0498 bool debug) {
0499 const CaloSubdetectorGeometry* barrelGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0500 const CaloSubdetectorGeometry* endcapGeom = (geo->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
0501
0502 std::vector<DetId> vdets = spr::matrixECALIds(coreDet, ieta, iphi, geo, caloTopology, debug);
0503 if (debug)
0504 edm::LogVerbatim("IsoTrack") << "chargeIsolation:: eta/phi/dets " << ieta << " " << iphi << " " << vdets.size();
0505
0506 double maxNearP = -1.0;
0507 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0508
0509
0510 reco::TrackCollection::const_iterator trkItr2;
0511 for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0512 const reco::Track* pTrack2 = &(*trkItr2);
0513
0514 bool trkQuality = pTrack2->quality(trackQuality_);
0515 if ((trkItr2 != trkItr) && trkQuality) {
0516 const FreeTrajectoryState fts2 =
0517 associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
0518 TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
0519 const GlobalPoint point2(info2.trkGlobPosAtEcal.x(), info2.trkGlobPosAtEcal.y(), info2.trkGlobPosAtEcal.z());
0520
0521 if (info2.isGoodEcal) {
0522 if (std::abs(point2.eta()) < spr::etaBEEcal) {
0523 const DetId anyCell = barrelGeom->getClosestCell(point2);
0524 if (debug)
0525 edm::LogVerbatim("IsoTrack")
0526 << "chargeIsolation:: EB cell " << (EBDetId)(anyCell) << " for pt " << pTrack2->p();
0527 if (!spr::chargeIsolation(anyCell, vdets)) {
0528 if (maxNearP < pTrack2->p())
0529 maxNearP = pTrack2->p();
0530 }
0531 } else {
0532 const DetId anyCell = endcapGeom->getClosestCell(point2);
0533 if (debug)
0534 edm::LogVerbatim("IsoTrack")
0535 << "chargeIsolation:: EE cell " << (EEDetId)(anyCell) << " for pt " << pTrack2->p();
0536 if (!spr::chargeIsolation(anyCell, vdets)) {
0537 if (maxNearP < pTrack2->p())
0538 maxNearP = pTrack2->p();
0539 }
0540 }
0541 }
0542 }
0543 }
0544 return maxNearP;
0545 }
0546
0547
0548
0549 double chargeIsolationHcal(const edm::Event& iEvent,
0550 const edm::EventSetup& iSetup,
0551 reco::TrackCollection::const_iterator trkItr,
0552 edm::Handle<reco::TrackCollection> trkCollection,
0553 const DetId ClosestCell,
0554 const HcalTopology* topology,
0555 const CaloSubdetectorGeometry* gHB,
0556 TrackDetectorAssociator& associator,
0557 TrackAssociatorParameters& parameters_,
0558 int ieta,
0559 int iphi,
0560 const std::string& theTrackQuality,
0561 bool debug) {
0562 std::vector<DetId> dets(1, ClosestCell);
0563
0564 if (debug)
0565 edm::LogVerbatim("IsoTrack") << static_cast<HcalDetId>(ClosestCell);
0566 std::vector<DetId> vdets = spr::matrixHCALIds(dets, topology, ieta, iphi, false, debug);
0567
0568 if (debug) {
0569 for (unsigned int i = 0; i < vdets.size(); i++) {
0570 edm::LogVerbatim("IsoTrack") << "HcalDetId in " << 2 * ieta + 1 << "x" << 2 * iphi + 1 << " "
0571 << static_cast<HcalDetId>(vdets[i]);
0572 }
0573 }
0574
0575 double maxNearP = -1.0;
0576 reco::TrackBase::TrackQuality trackQuality_ = reco::TrackBase::qualityByName(theTrackQuality);
0577
0578 reco::TrackCollection::const_iterator trkItr2;
0579 for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0580 const reco::Track* pTrack2 = &(*trkItr2);
0581
0582 bool trkQuality = pTrack2->quality(trackQuality_);
0583 if ((trkItr2 != trkItr) && trkQuality) {
0584 const FreeTrajectoryState fts2 =
0585 associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
0586 TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
0587 const GlobalPoint point2(info2.trkGlobPosAtHcal.x(), info2.trkGlobPosAtHcal.y(), info2.trkGlobPosAtHcal.z());
0588
0589 if (debug)
0590 edm::LogVerbatim("IsoTrack") << "Track2 (p,eta,phi) " << pTrack2->p() << " " << pTrack2->eta() << " "
0591 << pTrack2->phi();
0592
0593 if (info2.isGoodHcal) {
0594 const DetId anyCell = gHB->getClosestCell(point2);
0595 if (debug)
0596 edm::LogVerbatim("IsoTrack") << "chargeIsolation:: HCAL cell " << static_cast<HcalDetId>(anyCell)
0597 << " for pt " << pTrack2->p();
0598 if (!spr::chargeIsolation(anyCell, vdets)) {
0599 if (maxNearP < pTrack2->p())
0600 maxNearP = pTrack2->p();
0601 }
0602 if (debug) {
0603 edm::LogVerbatim("IsoTrack") << "maxNearP " << maxNearP << " thisCell " << static_cast<HcalDetId>(anyCell)
0604 << " (" << info2.trkGlobPosAtHcal.x() << "," << info2.trkGlobPosAtHcal.y()
0605 << "," << info2.trkGlobPosAtHcal.z() << ")";
0606 }
0607 }
0608 }
0609 }
0610 return maxNearP;
0611 }
0612
0613 }