Back to home page

Project CMSSW displayed by LXR
 
 

    

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     // const DetId anyCell,
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         }  //info.second
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     // Iterate over tracks
0264     reco::TrackCollection::const_iterator trkItr2;
0265     for (trkItr2 = trkCollection->begin(); trkItr2 != trkCollection->end(); ++trkItr2) {
0266       // Get track
0267       const reco::Track* pTrack2 = &(*trkItr2);
0268 
0269       // Get track qual, nlayers, and hit pattern
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       // Skip if the neighboring track candidate is the iso-track
0277       // candidate
0278       if (trkItr2 != trkItr) {
0279         // Get propagator
0280         const FreeTrajectoryState fts2 =
0281             associator.getFreeTrajectoryState(&iSetup.getData(parameters_.bFieldToken), *pTrack2);
0282         TrackDetMatchInfo info2 = associator.associate(iEvent, iSetup, fts2, parameters_);
0283 
0284         // Make sure it reaches Hcal
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     }  // Iterate over track loop
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     // const DetId anyCell,
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         }  //info2.isGoodEcal
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     // const DetId anyCell,
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         }  //info2.isGoodEcal
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 }  // 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!