Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 14:31:46 

0001 // -*- C++ -*-
0002 //
0003 // Package:    TrackAssociator
0004 // Class:      TrackDetectorAssociator
0005 //
0006 /*
0007 
0008  Description: <one line class summary>
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Dmytro Kovalskyi
0015 //         Created:  Fri Apr 21 10:59:41 PDT 2006
0016 //
0017 //
0018 
0019 #include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
0020 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0021 #include "TrackingTools/TrackAssociator/interface/DetIdInfo.h"
0022 #include "TrackingTools/Records/interface/DetIdAssociatorRecord.h"
0023 
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Utilities/interface/isFinite.h"
0026 #include "DataFormats/Common/interface/Handle.h"
0027 
0028 #include "DataFormats/DetId/interface/DetId.h"
0029 #include "DataFormats/TrackReco/interface/Track.h"
0030 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0031 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0032 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0033 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0034 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0035 #include "DataFormats/DTRecHit/interface/DTRecHitCollection.h"
0036 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0037 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0038 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
0039 #include "DataFormats/CSCRecHit/interface/CSCSegmentCollection.h"
0040 #include "DataFormats/GEMRecHit/interface/GEMSegmentCollection.h"
0041 #include "DataFormats/GEMRecHit/interface/ME0SegmentCollection.h"
0042 
0043 // calorimeter and muon infos
0044 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0045 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0046 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0047 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0048 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0049 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0050 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0051 
0052 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0053 #include <stack>
0054 #include <set>
0055 
0056 #include "DataFormats/Math/interface/LorentzVector.h"
0057 #include "Math/VectorUtil.h"
0058 #include <algorithm>
0059 
0060 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
0061 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
0062 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0063 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0064 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0065 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0066 
0067 using namespace reco;
0068 
0069 TrackDetectorAssociator::TrackDetectorAssociator() {
0070   ivProp_ = nullptr;
0071   useDefaultPropagator_ = false;
0072 }
0073 
0074 TrackDetectorAssociator::~TrackDetectorAssociator() = default;
0075 
0076 void TrackDetectorAssociator::setPropagator(const Propagator* ptr) {
0077   ivProp_ = ptr;
0078   cachedTrajectory_.setPropagator(ivProp_);
0079 }
0080 
0081 void TrackDetectorAssociator::useDefaultPropagator() { useDefaultPropagator_ = true; }
0082 
0083 void TrackDetectorAssociator::init(const edm::EventSetup& iSetup, const AssociatorParameters& parameters) {
0084   // access the calorimeter geometry
0085   theCaloGeometry_ = &iSetup.getData(parameters.theCaloGeometryToken);
0086 
0087   // get the tracking Geometry
0088   theTrackingGeometry_ = &iSetup.getData(parameters.theTrackingGeometryToken);
0089 
0090   if (useDefaultPropagator_ && (!defProp_ || theMagneticFieldWatcher_.check(iSetup))) {
0091     // setup propagator
0092     const MagneticField* bField = &iSetup.getData(parameters.bFieldToken);
0093 
0094     auto prop = std::make_unique<SteppingHelixPropagator>(bField, anyDirection);
0095     prop->setMaterialMode(false);
0096     prop->applyRadX0Correction(true);
0097     // prop->setDebug(true); // tmp
0098     defProp_ = std::move(prop);
0099     setPropagator(defProp_.get());
0100   }
0101 
0102   ecalDetIdAssociator_ = &iSetup.getData(parameters.ecalDetIdAssociatorToken);
0103   hcalDetIdAssociator_ = &iSetup.getData(parameters.hcalDetIdAssociatorToken);
0104   hoDetIdAssociator_ = &iSetup.getData(parameters.hoDetIdAssociatorToken);
0105   caloDetIdAssociator_ = &iSetup.getData(parameters.caloDetIdAssociatorToken);
0106   muonDetIdAssociator_ = &iSetup.getData(parameters.muonDetIdAssociatorToken);
0107   preshowerDetIdAssociator_ = &iSetup.getData(parameters.preshowerDetIdAssociatorToken);
0108 }
0109 
0110 TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent,
0111                                                      const edm::EventSetup& iSetup,
0112                                                      const FreeTrajectoryState& fts,
0113                                                      const AssociatorParameters& parameters) {
0114   return associate(iEvent, iSetup, parameters, &fts);
0115 }
0116 
0117 TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent,
0118                                                      const edm::EventSetup& iSetup,
0119                                                      const AssociatorParameters& parameters,
0120                                                      const FreeTrajectoryState* innerState,
0121                                                      const FreeTrajectoryState* outerState) {
0122   TrackDetMatchInfo info;
0123   if (!parameters.useEcal && !parameters.useCalo && !parameters.useHcal && !parameters.useHO && !parameters.useMuon &&
0124       !parameters.usePreshower)
0125     throw cms::Exception("ConfigurationError")
0126         << "Configuration error! No subdetector was selected for the track association.";
0127 
0128   SteppingHelixStateInfo trackOrigin(*innerState);
0129   info.stateAtIP = *innerState;
0130   cachedTrajectory_.setStateAtIP(trackOrigin);
0131 
0132   init(iSetup, parameters);
0133   // get track trajectory
0134   // ECAL points (EB+EE)
0135   // If the phi angle between a track entrance and exit points is more
0136   // than 2 crystals, it is possible that the track will cross 3 crystals
0137   // and therefore one has to check at least 3 points along the track
0138   // trajectory inside ECAL. In order to have a chance to cross 4 crystalls
0139   // in the barrel, a track should have P_t as low as 3 GeV or smaller
0140   // If it's necessary, number of points along trajectory can be increased
0141 
0142   info.setCaloGeometry(theCaloGeometry_);
0143 
0144   cachedTrajectory_.reset_trajectory();
0145   // estimate propagation outer boundaries based on
0146   // requested sub-detector information. For now limit
0147   // propagation region only if muon matching is not
0148   // requested.
0149   double HOmaxR = hoDetIdAssociator_->volume().maxR();
0150   double HOmaxZ = hoDetIdAssociator_->volume().maxZ();
0151   double minR = ecalDetIdAssociator_->volume().minR();
0152   double minZ = preshowerDetIdAssociator_->volume().minZ();
0153   cachedTrajectory_.setMaxHORadius(HOmaxR);
0154   cachedTrajectory_.setMaxHOLength(HOmaxZ * 2.);
0155   cachedTrajectory_.setMinDetectorRadius(minR);
0156   cachedTrajectory_.setMinDetectorLength(minZ * 2.);
0157 
0158   if (parameters.useMuon) {
0159     double maxR = muonDetIdAssociator_->volume().maxR();
0160     double maxZ = muonDetIdAssociator_->volume().maxZ();
0161     cachedTrajectory_.setMaxDetectorRadius(maxR);
0162     cachedTrajectory_.setMaxDetectorLength(maxZ * 2.);
0163   } else {
0164     cachedTrajectory_.setMaxDetectorRadius(HOmaxR);
0165     cachedTrajectory_.setMaxDetectorLength(HOmaxZ * 2.);
0166   }
0167 
0168   // If track extras exist and outerState is before HO maximum, then use outerState
0169   if (outerState) {
0170     if (outerState->position().perp() < HOmaxR && std::abs(outerState->position().z()) < HOmaxZ) {
0171       LogTrace("TrackAssociator") << "Using outerState as trackOrigin at Rho=" << outerState->position().perp()
0172                                   << "  Z=" << outerState->position().z() << "\n";
0173       trackOrigin = SteppingHelixStateInfo(*outerState);
0174     } else if (innerState) {
0175       LogTrace("TrackAssociator") << "Using innerState as trackOrigin at Rho=" << innerState->position().perp()
0176                                   << "  Z=" << innerState->position().z() << "\n";
0177       trackOrigin = SteppingHelixStateInfo(*innerState);
0178     }
0179   }
0180 
0181   if (trackOrigin.momentum().mag() == 0)
0182     return info;
0183   if (edm::isNotFinite(trackOrigin.momentum().x()) or edm::isNotFinite(trackOrigin.momentum().y()) or
0184       edm::isNotFinite(trackOrigin.momentum().z()))
0185     return info;
0186   if (!cachedTrajectory_.propagateAll(trackOrigin))
0187     return info;
0188 
0189   // get trajectory in calorimeters
0190   cachedTrajectory_.findEcalTrajectory(ecalDetIdAssociator_->volume());
0191   cachedTrajectory_.findHcalTrajectory(hcalDetIdAssociator_->volume());
0192   cachedTrajectory_.findHOTrajectory(hoDetIdAssociator_->volume());
0193   cachedTrajectory_.findPreshowerTrajectory(preshowerDetIdAssociator_->volume());
0194 
0195   info.trkGlobPosAtEcal = getPoint(cachedTrajectory_.getStateAtEcal().position());
0196   info.trkGlobPosAtHcal = getPoint(cachedTrajectory_.getStateAtHcal().position());
0197   info.trkGlobPosAtHO = getPoint(cachedTrajectory_.getStateAtHO().position());
0198 
0199   info.trkMomAtEcal = cachedTrajectory_.getStateAtEcal().momentum();
0200   info.trkMomAtHcal = cachedTrajectory_.getStateAtHcal().momentum();
0201   info.trkMomAtHO = cachedTrajectory_.getStateAtHO().momentum();
0202 
0203   if (parameters.useEcal)
0204     fillEcal(iEvent, info, parameters);
0205   if (parameters.useCalo)
0206     fillCaloTowers(iEvent, info, parameters);
0207   if (parameters.useHcal)
0208     fillHcal(iEvent, info, parameters);
0209   if (parameters.useHO)
0210     fillHO(iEvent, info, parameters);
0211   if (parameters.usePreshower)
0212     fillPreshower(iEvent, info, parameters);
0213   if (parameters.useMuon)
0214     fillMuon(iEvent, info, parameters);
0215   if (parameters.truthMatch)
0216     fillCaloTruth(iEvent, info, parameters);
0217 
0218   return info;
0219 }
0220 
0221 void TrackDetectorAssociator::fillEcal(const edm::Event& iEvent,
0222                                        TrackDetMatchInfo& info,
0223                                        const AssociatorParameters& parameters) {
0224   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getEcalTrajectory();
0225 
0226   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
0227        itr != trajectoryStates.end();
0228        itr++)
0229     LogTrace("TrackAssociator") << "ECAL trajectory point (rho, z, phi): " << itr->position().perp() << ", "
0230                                 << itr->position().z() << ", " << itr->position().phi();
0231 
0232   std::vector<GlobalPoint> coreTrajectory;
0233   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
0234        itr != trajectoryStates.end();
0235        itr++)
0236     coreTrajectory.push_back(itr->position());
0237 
0238   if (coreTrajectory.empty()) {
0239     LogTrace("TrackAssociator") << "ECAL track trajectory is empty; moving on\n";
0240     info.isGoodEcal = false;
0241     return;
0242   }
0243   info.isGoodEcal = true;
0244 
0245   // Find ECAL crystals
0246   edm::Handle<EBRecHitCollection> EBRecHits;
0247   iEvent.getByToken(parameters.EBRecHitsToken, EBRecHits);
0248   if (!EBRecHits.isValid())
0249     throw cms::Exception("FatalError") << "Unable to find EBRecHitCollection in the event!\n";
0250 
0251   edm::Handle<EERecHitCollection> EERecHits;
0252   iEvent.getByToken(parameters.EERecHitsToken, EERecHits);
0253   if (!EERecHits.isValid())
0254     throw cms::Exception("FatalError") << "Unable to find EERecHitCollection in event!\n";
0255 
0256   std::set<DetId> ecalIdsInRegion;
0257   if (parameters.accountForTrajectoryChangeCalo) {
0258     // get trajectory change with respect to initial state
0259     DetIdAssociator::MapRange mapRange =
0260         getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToEcal), parameters.dREcalPreselection);
0261     ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
0262   } else
0263     ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dREcalPreselection);
0264   LogTrace("TrackAssociator") << "ECAL hits in the region: " << ecalIdsInRegion.size();
0265   if (parameters.dREcalPreselection > parameters.dREcal)
0266     ecalIdsInRegion = ecalDetIdAssociator_->getDetIdsInACone(ecalIdsInRegion, coreTrajectory, parameters.dREcal);
0267   LogTrace("TrackAssociator") << "ECAL hits in the cone: " << ecalIdsInRegion.size();
0268   info.crossedEcalIds = ecalDetIdAssociator_->getCrossedDetIds(ecalIdsInRegion, coreTrajectory);
0269   const std::vector<DetId>& crossedEcalIds = info.crossedEcalIds;
0270   LogTrace("TrackAssociator") << "ECAL crossed hits " << crossedEcalIds.size();
0271 
0272   // add EcalRecHits
0273   for (std::vector<DetId>::const_iterator itr = crossedEcalIds.begin(); itr != crossedEcalIds.end(); itr++) {
0274     std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
0275     std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
0276     if (ebHit != (*EBRecHits).end())
0277       info.crossedEcalRecHits.push_back(&*ebHit);
0278     else if (eeHit != (*EERecHits).end())
0279       info.crossedEcalRecHits.push_back(&*eeHit);
0280     else
0281       LogTrace("TrackAssociator") << "Crossed EcalRecHit is not found for DetId: " << itr->rawId();
0282   }
0283   for (std::set<DetId>::const_iterator itr = ecalIdsInRegion.begin(); itr != ecalIdsInRegion.end(); itr++) {
0284     std::vector<EcalRecHit>::const_iterator ebHit = (*EBRecHits).find(*itr);
0285     std::vector<EcalRecHit>::const_iterator eeHit = (*EERecHits).find(*itr);
0286     if (ebHit != (*EBRecHits).end())
0287       info.ecalRecHits.push_back(&*ebHit);
0288     else if (eeHit != (*EERecHits).end())
0289       info.ecalRecHits.push_back(&*eeHit);
0290     else
0291       LogTrace("TrackAssociator") << "EcalRecHit from the cone is not found for DetId: " << itr->rawId();
0292   }
0293 }
0294 
0295 void TrackDetectorAssociator::fillCaloTowers(const edm::Event& iEvent,
0296                                              TrackDetMatchInfo& info,
0297                                              const AssociatorParameters& parameters) {
0298   // use ECAL and HCAL trajectories to match a tower. (HO isn't used for matching).
0299   std::vector<GlobalPoint> trajectory;
0300   const std::vector<SteppingHelixStateInfo>& ecalTrajectoryStates = cachedTrajectory_.getEcalTrajectory();
0301   const std::vector<SteppingHelixStateInfo>& hcalTrajectoryStates = cachedTrajectory_.getHcalTrajectory();
0302   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = ecalTrajectoryStates.begin();
0303        itr != ecalTrajectoryStates.end();
0304        itr++)
0305     trajectory.push_back(itr->position());
0306   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = hcalTrajectoryStates.begin();
0307        itr != hcalTrajectoryStates.end();
0308        itr++)
0309     trajectory.push_back(itr->position());
0310 
0311   if (trajectory.empty()) {
0312     LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
0313     info.isGoodCalo = false;
0314     return;
0315   }
0316   info.isGoodCalo = true;
0317 
0318   // find crossed CaloTowers
0319   edm::Handle<CaloTowerCollection> caloTowers;
0320   iEvent.getByToken(parameters.caloTowersToken, caloTowers);
0321   if (!caloTowers.isValid())
0322     throw cms::Exception("FatalError") << "Unable to find CaloTowers in event!\n";
0323 
0324   std::set<DetId> caloTowerIdsInRegion;
0325   if (parameters.accountForTrajectoryChangeCalo) {
0326     // get trajectory change with respect to initial state
0327     DetIdAssociator::MapRange mapRange =
0328         getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal), parameters.dRHcalPreselection);
0329     caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], mapRange);
0330   } else
0331     caloTowerIdsInRegion = caloDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRHcalPreselection);
0332 
0333   LogTrace("TrackAssociator") << "Towers in the region: " << caloTowerIdsInRegion.size();
0334 
0335   auto caloTowerIdsInAConeBegin = caloTowerIdsInRegion.begin();
0336   auto caloTowerIdsInAConeEnd = caloTowerIdsInRegion.end();
0337   std::set<DetId> caloTowerIdsInAConeTmp;
0338   if (!caloDetIdAssociator_->selectAllInACone(parameters.dRHcal)) {
0339     caloTowerIdsInAConeTmp =
0340         caloDetIdAssociator_->getDetIdsInACone(caloTowerIdsInRegion, trajectory, parameters.dRHcal);
0341     caloTowerIdsInAConeBegin = caloTowerIdsInAConeTmp.begin();
0342     caloTowerIdsInAConeEnd = caloTowerIdsInAConeTmp.end();
0343   }
0344   LogTrace("TrackAssociator") << "Towers in the cone: "
0345                               << std::distance(caloTowerIdsInAConeBegin, caloTowerIdsInAConeEnd);
0346 
0347   info.crossedTowerIds = caloDetIdAssociator_->getCrossedDetIds(caloTowerIdsInRegion, trajectory);
0348   const std::vector<DetId>& crossedCaloTowerIds = info.crossedTowerIds;
0349   LogTrace("TrackAssociator") << "Towers crossed: " << crossedCaloTowerIds.size();
0350 
0351   // add CaloTowers
0352   for (std::vector<DetId>::const_iterator itr = crossedCaloTowerIds.begin(); itr != crossedCaloTowerIds.end(); itr++) {
0353     CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
0354     if (tower != (*caloTowers).end())
0355       info.crossedTowers.push_back(&*tower);
0356     else
0357       LogTrace("TrackAssociator") << "Crossed CaloTower is not found for DetId: " << (*itr).rawId();
0358   }
0359 
0360   for (std::set<DetId>::const_iterator itr = caloTowerIdsInAConeBegin; itr != caloTowerIdsInAConeEnd; itr++) {
0361     CaloTowerCollection::const_iterator tower = (*caloTowers).find(*itr);
0362     if (tower != (*caloTowers).end())
0363       info.towers.push_back(&*tower);
0364     else
0365       LogTrace("TrackAssociator") << "CaloTower from the cone is not found for DetId: " << (*itr).rawId();
0366   }
0367 }
0368 
0369 void TrackDetectorAssociator::fillPreshower(const edm::Event& iEvent,
0370                                             TrackDetMatchInfo& info,
0371                                             const AssociatorParameters& parameters) {
0372   std::vector<GlobalPoint> trajectory;
0373   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getPreshowerTrajectory();
0374   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
0375        itr != trajectoryStates.end();
0376        itr++)
0377     trajectory.push_back(itr->position());
0378 
0379   if (trajectory.empty()) {
0380     LogTrace("TrackAssociator") << "Preshower trajectory is empty; moving on\n";
0381     return;
0382   }
0383 
0384   std::set<DetId> idsInRegion =
0385       preshowerDetIdAssociator_->getDetIdsCloseToAPoint(trajectory[0], parameters.dRPreshowerPreselection);
0386 
0387   LogTrace("TrackAssociator") << "Number of Preshower Ids in the region: " << idsInRegion.size();
0388   info.crossedPreshowerIds = preshowerDetIdAssociator_->getCrossedDetIds(idsInRegion, trajectory);
0389   LogTrace("TrackAssociator") << "Number of Preshower Ids in crossed: " << info.crossedPreshowerIds.size();
0390 }
0391 
0392 void TrackDetectorAssociator::fillHcal(const edm::Event& iEvent,
0393                                        TrackDetMatchInfo& info,
0394                                        const AssociatorParameters& parameters) {
0395   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHcalTrajectory();
0396 
0397   std::vector<GlobalPoint> coreTrajectory;
0398   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
0399        itr != trajectoryStates.end();
0400        itr++)
0401     coreTrajectory.push_back(itr->position());
0402 
0403   if (coreTrajectory.empty()) {
0404     LogTrace("TrackAssociator") << "HCAL trajectory is empty; moving on\n";
0405     info.isGoodHcal = false;
0406     return;
0407   }
0408   info.isGoodHcal = true;
0409 
0410   // find crossed Hcals
0411   edm::Handle<HBHERecHitCollection> collection;
0412   iEvent.getByToken(parameters.HBHEcollToken, collection);
0413   if (!collection.isValid())
0414     throw cms::Exception("FatalError") << "Unable to find HBHERecHits in event!\n";
0415 
0416   std::set<DetId> idsInRegion;
0417   if (parameters.accountForTrajectoryChangeCalo) {
0418     // get trajectory change with respect to initial state
0419     DetIdAssociator::MapRange mapRange =
0420         getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHcal), parameters.dRHcalPreselection);
0421     idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
0422   } else
0423     idsInRegion = hcalDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
0424 
0425   LogTrace("TrackAssociator") << "HCAL hits in the region: " << idsInRegion.size() << "\n"
0426                               << DetIdInfo::info(idsInRegion, nullptr);
0427 
0428   auto idsInAConeBegin = idsInRegion.begin();
0429   auto idsInAConeEnd = idsInRegion.end();
0430   std::set<DetId> idsInAConeTmp;
0431   if (!hcalDetIdAssociator_->selectAllInACone(parameters.dRHcal)) {
0432     idsInAConeTmp = hcalDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
0433     idsInAConeBegin = idsInAConeTmp.begin();
0434     idsInAConeEnd = idsInAConeTmp.end();
0435   }
0436   LogTrace("TrackAssociator") << "HCAL hits in the cone: " << std::distance(idsInAConeBegin, idsInAConeEnd) << "\n"
0437                               << DetIdInfo::info(std::set<DetId>(idsInAConeBegin, idsInAConeEnd), nullptr);
0438   info.crossedHcalIds = hcalDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
0439   const std::vector<DetId>& crossedIds = info.crossedHcalIds;
0440   LogTrace("TrackAssociator") << "HCAL hits crossed: " << crossedIds.size() << "\n"
0441                               << DetIdInfo::info(crossedIds, nullptr);
0442 
0443   // add Hcal
0444   for (std::vector<DetId>::const_iterator itr = crossedIds.begin(); itr != crossedIds.end(); itr++) {
0445     HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
0446     if (hit != (*collection).end())
0447       info.crossedHcalRecHits.push_back(&*hit);
0448     else
0449       LogTrace("TrackAssociator") << "Crossed HBHERecHit is not found for DetId: " << itr->rawId();
0450   }
0451   for (std::set<DetId>::const_iterator itr = idsInAConeBegin; itr != idsInAConeEnd; itr++) {
0452     HBHERecHitCollection::const_iterator hit = (*collection).find(*itr);
0453     if (hit != (*collection).end())
0454       info.hcalRecHits.push_back(&*hit);
0455     else
0456       LogTrace("TrackAssociator") << "HBHERecHit from the cone is not found for DetId: " << itr->rawId();
0457   }
0458 }
0459 
0460 void TrackDetectorAssociator::fillHO(const edm::Event& iEvent,
0461                                      TrackDetMatchInfo& info,
0462                                      const AssociatorParameters& parameters) {
0463   const std::vector<SteppingHelixStateInfo>& trajectoryStates = cachedTrajectory_.getHOTrajectory();
0464 
0465   std::vector<GlobalPoint> coreTrajectory;
0466   for (std::vector<SteppingHelixStateInfo>::const_iterator itr = trajectoryStates.begin();
0467        itr != trajectoryStates.end();
0468        itr++)
0469     coreTrajectory.push_back(itr->position());
0470 
0471   if (coreTrajectory.empty()) {
0472     LogTrace("TrackAssociator") << "HO trajectory is empty; moving on\n";
0473     info.isGoodHO = false;
0474     return;
0475   }
0476   info.isGoodHO = true;
0477 
0478   // find crossed HOs
0479   edm::Handle<HORecHitCollection> collection;
0480   iEvent.getByToken(parameters.HOcollToken, collection);
0481   if (!collection.isValid())
0482     throw cms::Exception("FatalError") << "Unable to find HORecHits in event!\n";
0483 
0484   std::set<DetId> idsInRegion;
0485   if (parameters.accountForTrajectoryChangeCalo) {
0486     // get trajectory change with respect to initial state
0487     DetIdAssociator::MapRange mapRange =
0488         getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::IpToHO), parameters.dRHcalPreselection);
0489     idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], mapRange);
0490   } else
0491     idsInRegion = hoDetIdAssociator_->getDetIdsCloseToAPoint(coreTrajectory[0], parameters.dRHcalPreselection);
0492 
0493   LogTrace("TrackAssociator") << "idsInRegion.size(): " << idsInRegion.size();
0494 
0495   auto idsInAConeBegin = idsInRegion.begin();
0496   auto idsInAConeEnd = idsInRegion.end();
0497   std::set<DetId> idsInAConeTmp;
0498   if (!hoDetIdAssociator_->selectAllInACone(parameters.dRHcal)) {
0499     idsInAConeTmp = hoDetIdAssociator_->getDetIdsInACone(idsInRegion, coreTrajectory, parameters.dRHcal);
0500     idsInAConeBegin = idsInAConeTmp.begin();
0501     idsInAConeEnd = idsInAConeTmp.end();
0502   }
0503   LogTrace("TrackAssociator") << "idsInACone.size(): " << std::distance(idsInAConeBegin, idsInAConeEnd);
0504   info.crossedHOIds = hoDetIdAssociator_->getCrossedDetIds(idsInRegion, coreTrajectory);
0505   const std::vector<DetId>& crossedIds = info.crossedHOIds;
0506 
0507   // add HO
0508   for (std::vector<DetId>::const_iterator itr = crossedIds.begin(); itr != crossedIds.end(); itr++) {
0509     HORecHitCollection::const_iterator hit = (*collection).find(*itr);
0510     if (hit != (*collection).end())
0511       info.crossedHORecHits.push_back(&*hit);
0512     else
0513       LogTrace("TrackAssociator") << "Crossed HORecHit is not found for DetId: " << itr->rawId();
0514   }
0515 
0516   for (std::set<DetId>::const_iterator itr = idsInAConeBegin; itr != idsInAConeEnd; itr++) {
0517     HORecHitCollection::const_iterator hit = (*collection).find(*itr);
0518     if (hit != (*collection).end())
0519       info.hoRecHits.push_back(&*hit);
0520     else
0521       LogTrace("TrackAssociator") << "HORecHit from the cone is not found for DetId: " << itr->rawId();
0522   }
0523 }
0524 
0525 FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState(const MagneticField* bField,
0526                                                                     const SimTrack& track,
0527                                                                     const SimVertex& vertex) {
0528   GlobalVector vector(track.momentum().x(), track.momentum().y(), track.momentum().z());
0529   GlobalPoint point(vertex.position().x(), vertex.position().y(), vertex.position().z());
0530 
0531   int charge = track.type() > 0 ? -1 : 1;  // lepton convention
0532   if (abs(track.type()) == 211 ||          // pion
0533       abs(track.type()) == 321 ||          // kaon
0534       abs(track.type()) == 2212)
0535     charge = track.type() < 0 ? -1 : 1;
0536   return getFreeTrajectoryState(bField, vector, point, charge);
0537 }
0538 
0539 FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState(const MagneticField* bField,
0540                                                                     const GlobalVector& momentum,
0541                                                                     const GlobalPoint& vertex,
0542                                                                     const int charge) {
0543   GlobalTrajectoryParameters tPars(vertex, momentum, charge, bField);
0544 
0545   ROOT::Math::SMatrixIdentity id;
0546   AlgebraicSymMatrix66 covT(id);
0547   covT *= 1e-6;  // initialize to sigma=1e-3
0548   CartesianTrajectoryError tCov(covT);
0549 
0550   return FreeTrajectoryState(tPars, tCov);
0551 }
0552 
0553 FreeTrajectoryState TrackDetectorAssociator::getFreeTrajectoryState(const MagneticField* bField,
0554                                                                     const reco::Track& track) {
0555   GlobalVector vector(track.momentum().x(), track.momentum().y(), track.momentum().z());
0556 
0557   GlobalPoint point(track.vertex().x(), track.vertex().y(), track.vertex().z());
0558 
0559   GlobalTrajectoryParameters tPars(point, vector, track.charge(), bField);
0560 
0561   // FIX THIS !!!
0562   // need to convert from perigee to global or helix (curvilinear) frame
0563   // for now just an arbitrary matrix.
0564   ROOT::Math::SMatrixIdentity id;
0565   AlgebraicSymMatrix66 covT(id);
0566   covT *= 1e-6;  // initialize to sigma=1e-3
0567   CartesianTrajectoryError tCov(covT);
0568 
0569   return FreeTrajectoryState(tPars, tCov);
0570 }
0571 
0572 DetIdAssociator::MapRange TrackDetectorAssociator::getMapRange(const std::pair<float, float>& delta, const float dR) {
0573   DetIdAssociator::MapRange mapRange;
0574   mapRange.dThetaPlus = dR;
0575   mapRange.dThetaMinus = dR;
0576   mapRange.dPhiPlus = dR;
0577   mapRange.dPhiMinus = dR;
0578   if (delta.first > 0)
0579     mapRange.dThetaPlus += delta.first;
0580   else
0581     mapRange.dThetaMinus += std::abs(delta.first);
0582   if (delta.second > 0)
0583     mapRange.dPhiPlus += delta.second;
0584   else
0585     mapRange.dPhiMinus += std::abs(delta.second);
0586   LogTrace("TrackAssociator") << "Selection range: (dThetaPlus, dThetaMinus, dPhiPlus, dPhiMinus, dRPreselection): "
0587                               << mapRange.dThetaPlus << ", " << mapRange.dThetaMinus << ", " << mapRange.dPhiPlus
0588                               << ", " << mapRange.dPhiMinus << ", " << dR;
0589   return mapRange;
0590 }
0591 
0592 void TrackDetectorAssociator::getTAMuonChamberMatches(std::vector<TAMuonChamberMatch>& matches,
0593                                                       const AssociatorParameters& parameters) {
0594   // Strategy:
0595   //    Propagate through the whole detector, estimate change in eta and phi
0596   //    along the trajectory, add this to dRMuon and find DetIds around this
0597   //    direction using the map. Then propagate fast to each surface and apply
0598   //    final matching criteria.
0599 
0600   // get the direction first
0601   SteppingHelixStateInfo trajectoryPoint = cachedTrajectory_.getStateAtHcal();
0602   // If trajectory point at HCAL is not valid, try to use the outer most state of the
0603   // trajectory instead.
0604   if (!trajectoryPoint.isValid())
0605     trajectoryPoint = cachedTrajectory_.getOuterState();
0606   if (!trajectoryPoint.isValid()) {
0607     LogTrace("TrackAssociator")
0608         << "trajectory position at HCAL is not valid. Assume the track cannot reach muon detectors and skip it";
0609     return;
0610   }
0611 
0612   GlobalVector direction = trajectoryPoint.momentum().unit();
0613   LogTrace("TrackAssociator") << "muon direction: " << direction
0614                               << "\n\t and corresponding point: " << trajectoryPoint.position() << "\n";
0615 
0616   DetIdAssociator::MapRange mapRange =
0617       getMapRange(cachedTrajectory_.trajectoryDelta(CachedTrajectory::FullTrajectory), parameters.dRMuonPreselection);
0618 
0619   // and find chamber DetIds
0620 
0621   std::set<DetId> muonIdsInRegion = muonDetIdAssociator_->getDetIdsCloseToAPoint(trajectoryPoint.position(), mapRange);
0622   LogTrace("TrackAssociator") << "Number of chambers to check: " << muonIdsInRegion.size();
0623   for (std::set<DetId>::const_iterator detId = muonIdsInRegion.begin(); detId != muonIdsInRegion.end(); detId++) {
0624     const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet(*detId);
0625     TrajectoryStateOnSurface stateOnSurface = cachedTrajectory_.propagate(&geomDet->surface());
0626     if (!stateOnSurface.isValid()) {
0627       LogTrace("TrackAssociator") << "Failed to propagate the track; moving on\n\t"
0628                                   << "Element is not crosssed: " << DetIdInfo::info(*detId, nullptr) << "\n";
0629       continue;
0630     }
0631     LocalPoint localPoint = geomDet->surface().toLocal(stateOnSurface.freeState()->position());
0632     LocalError localError = stateOnSurface.localError().positionError();
0633     float distanceX = 0.f;
0634     float distanceY = 0.f;
0635     if (const CSCChamber* cscChamber = dynamic_cast<const CSCChamber*>(geomDet)) {
0636       const CSCChamberSpecs* chamberSpecs = cscChamber->specs();
0637       if (!chamberSpecs) {
0638         LogTrace("TrackAssociator") << "Failed to get CSCChamberSpecs from CSCChamber; moving on\n";
0639         continue;
0640       }
0641       const CSCLayerGeometry* layerGeometry = chamberSpecs->oddLayerGeometry(1);
0642       if (!layerGeometry) {
0643         LogTrace("TrackAssociator") << "Failed to get CSCLayerGeometry from CSCChamberSpecs; moving on\n";
0644         continue;
0645       }
0646       const CSCWireTopology* wireTopology = layerGeometry->wireTopology();
0647       if (!wireTopology) {
0648         LogTrace("TrackAssociator") << "Failed to get CSCWireTopology from CSCLayerGeometry; moving on\n";
0649         continue;
0650       }
0651 
0652       float wideWidth = wireTopology->wideWidthOfPlane();
0653       float narrowWidth = wireTopology->narrowWidthOfPlane();
0654       float length = wireTopology->lengthOfPlane();
0655       // If slanted, there is no y offset between local origin and symmetry center of wire plane
0656       float yOfFirstWire = std::abs(wireTopology->wireAngle()) > 1.E-06f ? -0.5 * length : wireTopology->yOfWire(1);
0657       // y offset between local origin and symmetry center of wire plane
0658       float yCOWPOffset = yOfFirstWire + 0.5f * length;
0659 
0660       // tangent of the incline angle from inside the trapezoid
0661       float tangent = (wideWidth - narrowWidth) / (2.f * length);
0662       // y position wrt bottom of trapezoid
0663       float yPrime = localPoint.y() + std::abs(yOfFirstWire);
0664       // half trapezoid width at y' is 0.5 * narrowWidth + x side of triangle with the above tangent and side y'
0665       float halfWidthAtYPrime = 0.5f * narrowWidth + yPrime * tangent;
0666       distanceX = std::abs(localPoint.x()) - halfWidthAtYPrime;
0667       distanceY = std::abs(localPoint.y() - yCOWPOffset) - 0.5f * length;
0668     } else {
0669       distanceX = std::abs(localPoint.x()) - 0.5f * geomDet->surface().bounds().width();
0670       distanceY = std::abs(localPoint.y()) - 0.5f * geomDet->surface().bounds().length();
0671     }
0672     if ((distanceX < parameters.muonMaxDistanceX && distanceY < parameters.muonMaxDistanceY) ||
0673         (distanceX * distanceX <
0674              localError.xx() * parameters.muonMaxDistanceSigmaX * parameters.muonMaxDistanceSigmaX &&
0675          distanceY * distanceY <
0676              localError.yy() * parameters.muonMaxDistanceSigmaY * parameters.muonMaxDistanceSigmaY)) {
0677       LogTrace("TrackAssociator") << "found a match: " << DetIdInfo::info(*detId, nullptr) << "\n";
0678       TAMuonChamberMatch match;
0679       match.tState = stateOnSurface;
0680       match.localDistanceX = distanceX;
0681       match.localDistanceY = distanceY;
0682       match.id = *detId;
0683       matches.push_back(match);
0684     } else {
0685       LogTrace("TrackAssociator") << "chamber is too far: " << DetIdInfo::info(*detId, nullptr)
0686                                   << "\n\tdistanceX: " << distanceX << "\t distanceY: " << distanceY
0687                                   << "\t sigmaX: " << distanceX / sqrt(localError.xx())
0688                                   << "\t sigmaY: " << distanceY / sqrt(localError.yy()) << "\n";
0689     }
0690   }
0691 }
0692 
0693 void TrackDetectorAssociator::fillMuon(const edm::Event& iEvent,
0694                                        TrackDetMatchInfo& info,
0695                                        const AssociatorParameters& parameters) {
0696   // Get the segments from the event
0697   edm::Handle<DTRecSegment4DCollection> dtSegments;
0698   iEvent.getByToken(parameters.dtSegmentsToken, dtSegments);
0699   if (!dtSegments.isValid())
0700     throw cms::Exception("FatalError") << "Unable to find DTRecSegment4DCollection in event!\n";
0701 
0702   edm::Handle<CSCSegmentCollection> cscSegments;
0703   iEvent.getByToken(parameters.cscSegmentsToken, cscSegments);
0704   if (!cscSegments.isValid())
0705     throw cms::Exception("FatalError") << "Unable to find CSCSegmentCollection in event!\n";
0706 
0707   edm::Handle<GEMSegmentCollection> gemSegments;
0708   if (parameters.useGEM)
0709     iEvent.getByToken(parameters.gemSegmentsToken, gemSegments);
0710   edm::Handle<ME0SegmentCollection> me0Segments;
0711   if (parameters.useME0)
0712     iEvent.getByToken(parameters.me0SegmentsToken, me0Segments);
0713 
0714   ///// get a set of DetId's in a given direction
0715 
0716   // check the map of available segments
0717   // if there is no segments in a given direction at all,
0718   // then there is no point to fly there.
0719   //
0720   // MISSING
0721   // Possible solution: quick search for presence of segments
0722   // for the set of DetIds
0723 
0724   // get a set of matches corresponding to muon chambers
0725   std::vector<TAMuonChamberMatch> matchedChambers;
0726   LogTrace("TrackAssociator") << "Trying to Get ChamberMatches" << std::endl;
0727   getTAMuonChamberMatches(matchedChambers, parameters);
0728   LogTrace("TrackAssociator") << "Chambers matched: " << matchedChambers.size() << "\n";
0729 
0730   // Iterate over all chamber matches and fill segment matching
0731   // info if it's available
0732   for (std::vector<TAMuonChamberMatch>::iterator matchedChamber = matchedChambers.begin();
0733        matchedChamber != matchedChambers.end();
0734        matchedChamber++) {
0735     const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet((*matchedChamber).id);
0736     // DT chamber
0737     if (const DTChamber* chamber = dynamic_cast<const DTChamber*>(geomDet)) {
0738       // Get the range for the corresponding segments
0739       DTRecSegment4DCollection::range range = dtSegments->get(chamber->id());
0740       // Loop over the segments of this chamber
0741       for (DTRecSegment4DCollection::const_iterator segment = range.first; segment != range.second; segment++) {
0742         if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
0743           matchedChamber->segments.back().dtSegmentRef = DTRecSegment4DRef(dtSegments, segment - dtSegments->begin());
0744         }
0745       }
0746     }
0747     // CSC Chamber
0748     else if (const CSCChamber* chamber = dynamic_cast<const CSCChamber*>(geomDet)) {
0749       // Get the range for the corresponding segments
0750       CSCSegmentCollection::range range = cscSegments->get(chamber->id());
0751       // Loop over the segments
0752       for (CSCSegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) {
0753         if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
0754           matchedChamber->segments.back().cscSegmentRef = CSCSegmentRef(cscSegments, segment - cscSegments->begin());
0755         }
0756       }
0757     } else {
0758       // GEM Chamber
0759       if (parameters.useGEM) {
0760         if (const GEMSuperChamber* chamber = dynamic_cast<const GEMSuperChamber*>(geomDet)) {
0761           // Get the range for the corresponding segments
0762           GEMSegmentCollection::range range = gemSegments->get(chamber->id());
0763           // Loop over the segments
0764           for (GEMSegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) {
0765             if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
0766               matchedChamber->segments.back().gemSegmentRef =
0767                   GEMSegmentRef(gemSegments, segment - gemSegments->begin());
0768             }
0769           }
0770         }
0771       }
0772       // ME0 Chamber
0773       if (parameters.useME0) {
0774         if (const ME0Chamber* chamber = dynamic_cast<const ME0Chamber*>(geomDet)) {
0775           // Get the range for the corresponding segments
0776           ME0SegmentCollection::range range = me0Segments->get(chamber->id());
0777           // Loop over the segments
0778           for (ME0SegmentCollection::const_iterator segment = range.first; segment != range.second; segment++) {
0779             if (addTAMuonSegmentMatch(*matchedChamber, &(*segment), parameters)) {
0780               matchedChamber->segments.back().me0SegmentRef =
0781                   ME0SegmentRef(me0Segments, segment - me0Segments->begin());
0782             }
0783           }
0784         }
0785       }
0786     }
0787     info.chambers.push_back(*matchedChamber);
0788   }
0789 }
0790 
0791 bool TrackDetectorAssociator::addTAMuonSegmentMatch(TAMuonChamberMatch& matchedChamber,
0792                                                     const RecSegment* segment,
0793                                                     const AssociatorParameters& parameters) {
0794   LogTrace("TrackAssociator") << "Segment local position: " << segment->localPosition() << "\n"
0795                               << std::hex << segment->geographicalId().rawId() << "\n";
0796 
0797   const GeomDet* chamber = muonDetIdAssociator_->getGeomDet(matchedChamber.id);
0798   TrajectoryStateOnSurface trajectoryStateOnSurface = matchedChamber.tState;
0799   GlobalPoint segmentGlobalPosition = chamber->toGlobal(segment->localPosition());
0800 
0801   LogTrace("TrackAssociator") << "Segment global position: " << segmentGlobalPosition
0802                               << " \t (R_xy,eta,phi): " << segmentGlobalPosition.perp() << ","
0803                               << segmentGlobalPosition.eta() << "," << segmentGlobalPosition.phi() << "\n";
0804 
0805   LogTrace("TrackAssociator") << "\teta hit: " << segmentGlobalPosition.eta()
0806                               << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().eta() << "\n"
0807                               << "\tphi hit: " << segmentGlobalPosition.phi()
0808                               << " \tpropagator: " << trajectoryStateOnSurface.freeState()->position().phi()
0809                               << std::endl;
0810 
0811   bool isGood = false;
0812   bool isDTWithoutY = false;
0813   const DTRecSegment4D* dtseg = dynamic_cast<const DTRecSegment4D*>(segment);
0814   if (dtseg && (!dtseg->hasZed()))
0815     isDTWithoutY = true;
0816 
0817   float deltaPhi(std::abs(segmentGlobalPosition.phi() - trajectoryStateOnSurface.freeState()->position().phi()));
0818   if (deltaPhi > float(M_PI))
0819     deltaPhi = std::abs(deltaPhi - float(M_PI) * 2.f);
0820   float deltaEta = std::abs(segmentGlobalPosition.eta() - trajectoryStateOnSurface.freeState()->position().eta());
0821 
0822   if (isDTWithoutY) {
0823     isGood = deltaPhi < parameters.dRMuon;
0824     // Be in chamber
0825     isGood &= deltaEta < .3f;
0826   } else
0827     isGood = deltaEta * deltaEta + deltaPhi * deltaPhi < parameters.dRMuon * parameters.dRMuon;
0828 
0829   if (isGood) {
0830     TAMuonSegmentMatch muonSegment;
0831     muonSegment.segmentGlobalPosition = getPoint(segmentGlobalPosition);
0832     muonSegment.segmentLocalPosition = getPoint(segment->localPosition());
0833     muonSegment.segmentLocalDirection = getVector(segment->localDirection());
0834     muonSegment.segmentLocalErrorXX = segment->localPositionError().xx();
0835     muonSegment.segmentLocalErrorYY = segment->localPositionError().yy();
0836     muonSegment.segmentLocalErrorXY = segment->localPositionError().xy();
0837     muonSegment.segmentLocalErrorDxDz = segment->localDirectionError().xx();
0838     muonSegment.segmentLocalErrorDyDz = segment->localDirectionError().yy();
0839 
0840     // DANGEROUS - compiler cannot guaranty parameters ordering
0841     // AlgebraicSymMatrix segmentCovMatrix = segment->parametersError();
0842     // muonSegment.segmentLocalErrorXDxDz = segmentCovMatrix[2][0];
0843     // muonSegment.segmentLocalErrorYDyDz = segmentCovMatrix[3][1];
0844     muonSegment.segmentLocalErrorXDxDz = -999.f;
0845     muonSegment.segmentLocalErrorYDyDz = -999.f;
0846     muonSegment.hasZed = true;
0847     muonSegment.hasPhi = true;
0848 
0849     // timing information
0850     muonSegment.t0 = 0.f;
0851     if (dtseg) {
0852       if ((dtseg->hasPhi()) && (!isDTWithoutY)) {
0853         int phiHits = dtseg->phiSegment()->specificRecHits().size();
0854         //    int zHits = dtseg->zSegment()->specificRecHits().size();
0855         int hits = 0;
0856         double t0 = 0.;
0857         // TODO: cuts on hit numbers not optimized in any way yet...
0858         if (phiHits > 5 && dtseg->phiSegment()->ist0Valid()) {
0859           t0 += dtseg->phiSegment()->t0() * phiHits;
0860           hits += phiHits;
0861           LogTrace("TrackAssociator") << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
0862         }
0863         // the z segments seem to contain little useful information...
0864         //    if (zHits>3) {
0865         //      t0+=s->zSegment()->t0()*zHits;
0866         //      hits+=zHits;
0867         //      LogTrace("TrackAssociator") << "   Z t0: " << s->zSegment()->t0() << " hits: " << zHits << std::endl;
0868         //    }
0869         if (hits)
0870           muonSegment.t0 = t0 / hits;
0871         //    LogTrace("TrackAssociator") << " --- t0: " << muonSegment.t0 << std::endl;
0872       } else {
0873         // check and set dimensionality
0874         if (isDTWithoutY)
0875           muonSegment.hasZed = false;
0876         if (!dtseg->hasPhi())
0877           muonSegment.hasPhi = false;
0878       }
0879     }
0880     matchedChamber.segments.push_back(muonSegment);
0881   }
0882 
0883   return isGood;
0884 }
0885 
0886 //********************** NON-CORE CODE ******************************//
0887 
0888 void TrackDetectorAssociator::fillCaloTruth(const edm::Event& iEvent,
0889                                             TrackDetMatchInfo& info,
0890                                             const AssociatorParameters& parameters) {
0891   // get list of simulated tracks and their vertices
0892   using namespace edm;
0893   Handle<SimTrackContainer> simTracks;
0894   iEvent.getByToken(parameters.simTracksToken, simTracks);
0895   if (!simTracks.isValid())
0896     throw cms::Exception("FatalError") << "No simulated tracks found\n";
0897 
0898   Handle<SimVertexContainer> simVertices;
0899   iEvent.getByToken(parameters.simVerticesToken, simVertices);
0900   if (!simVertices.isValid())
0901     throw cms::Exception("FatalError") << "No simulated vertices found\n";
0902 
0903   // get sim calo hits
0904   Handle<PCaloHitContainer> simEcalHitsEB;
0905   iEvent.getByToken(parameters.simEcalHitsEBToken, simEcalHitsEB);
0906   if (!simEcalHitsEB.isValid())
0907     throw cms::Exception("FatalError") << "No simulated ECAL EB hits found\n";
0908 
0909   Handle<PCaloHitContainer> simEcalHitsEE;
0910   iEvent.getByToken(parameters.simEcalHitsEEToken, simEcalHitsEE);
0911   if (!simEcalHitsEE.isValid())
0912     throw cms::Exception("FatalError") << "No simulated ECAL EE hits found\n";
0913 
0914   Handle<PCaloHitContainer> simHcalHits;
0915   iEvent.getByToken(parameters.simHcalHitsToken, simHcalHits);
0916   if (!simHcalHits.isValid())
0917     throw cms::Exception("FatalError") << "No simulated HCAL hits found\n";
0918 
0919   // find truth partner
0920   SimTrackContainer::const_iterator simTrack = simTracks->begin();
0921   for (; simTrack != simTracks->end(); ++simTrack) {
0922     math::XYZVector simP3(simTrack->momentum().x(), simTrack->momentum().y(), simTrack->momentum().z());
0923     math::XYZVector recoP3(info.stateAtIP.momentum().x(), info.stateAtIP.momentum().y(), info.stateAtIP.momentum().z());
0924     if (ROOT::Math::VectorUtil::DeltaR(recoP3, simP3) < 0.1)
0925       break;
0926   }
0927   if (simTrack != simTracks->end()) {
0928     info.simTrack = &(*simTrack);
0929     float ecalTrueEnergy(0);
0930     float hcalTrueEnergy(0);
0931 
0932     // loop over calo hits
0933     for (PCaloHitContainer::const_iterator hit = simEcalHitsEB->begin(); hit != simEcalHitsEB->end(); ++hit)
0934       if (hit->geantTrackId() == info.simTrack->genpartIndex())
0935         ecalTrueEnergy += hit->energy();
0936 
0937     for (PCaloHitContainer::const_iterator hit = simEcalHitsEE->begin(); hit != simEcalHitsEE->end(); ++hit)
0938       if (hit->geantTrackId() == info.simTrack->genpartIndex())
0939         ecalTrueEnergy += hit->energy();
0940 
0941     for (PCaloHitContainer::const_iterator hit = simHcalHits->begin(); hit != simHcalHits->end(); ++hit)
0942       if (hit->geantTrackId() == info.simTrack->genpartIndex())
0943         hcalTrueEnergy += hit->energy();
0944 
0945     info.ecalTrueEnergy = ecalTrueEnergy;
0946     info.hcalTrueEnergy = hcalTrueEnergy;
0947     info.hcalTrueEnergyCorrected = hcalTrueEnergy;
0948     if (std::abs(info.trkGlobPosAtHcal.eta()) < 1.3f)
0949       info.hcalTrueEnergyCorrected = hcalTrueEnergy * 113.2f;
0950     else if (std::abs(info.trkGlobPosAtHcal.eta()) < 3.0f)
0951       info.hcalTrueEnergyCorrected = hcalTrueEnergy * 167.2f;
0952   }
0953 }
0954 
0955 TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent,
0956                                                      const edm::EventSetup& iSetup,
0957                                                      const reco::Track& track,
0958                                                      const AssociatorParameters& parameters,
0959                                                      Direction direction /*= Any*/) {
0960   double currentStepSize = cachedTrajectory_.getPropagationStep();
0961 
0962   const MagneticField* bField = &iSetup.getData(parameters.bFieldToken);
0963 
0964   if (track.extra().isNull()) {
0965     if (direction != InsideOut)
0966       throw cms::Exception("FatalError") << "No TrackExtra information is available and association is done with "
0967                                             "something else than InsideOut track.\n"
0968                                          << "Either change the parameter or provide needed data!\n";
0969     LogTrace("TrackAssociator") << "Track Extras not found\n";
0970     FreeTrajectoryState initialState = trajectoryStateTransform::initialFreeState(track, bField);
0971     return associate(iEvent, iSetup, parameters, &initialState);  // 5th argument is null pointer
0972   }
0973 
0974   LogTrace("TrackAssociator") << "Track Extras found\n";
0975   FreeTrajectoryState innerState = trajectoryStateTransform::innerFreeState(track, bField);
0976   FreeTrajectoryState outerState = trajectoryStateTransform::outerFreeState(track, bField);
0977   FreeTrajectoryState referenceState = trajectoryStateTransform::initialFreeState(track, bField);
0978 
0979   LogTrace("TrackAssociator") << "inner track state (rho, z, phi):" << track.innerPosition().Rho() << ", "
0980                               << track.innerPosition().z() << ", " << track.innerPosition().phi() << "\n";
0981   LogTrace("TrackAssociator") << "innerFreeState (rho, z, phi):" << innerState.position().perp() << ", "
0982                               << innerState.position().z() << ", " << innerState.position().phi() << "\n";
0983 
0984   LogTrace("TrackAssociator") << "outer track state (rho, z, phi):" << track.outerPosition().Rho() << ", "
0985                               << track.outerPosition().z() << ", " << track.outerPosition().phi() << "\n";
0986   LogTrace("TrackAssociator") << "outerFreeState (rho, z, phi):" << outerState.position().perp() << ", "
0987                               << outerState.position().z() << ", " << outerState.position().phi() << "\n";
0988 
0989   // InsideOut first
0990   if (crossedIP(track)) {
0991     switch (direction) {
0992       case InsideOut:
0993       case Any:
0994         return associate(iEvent, iSetup, parameters, &referenceState, &outerState);
0995         break;
0996       case OutsideIn: {
0997         cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize));
0998         TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &innerState, &referenceState);
0999         cachedTrajectory_.setPropagationStep(currentStepSize);
1000         return result;
1001         break;
1002       }
1003     }
1004   } else {
1005     switch (direction) {
1006       case InsideOut:
1007         return associate(iEvent, iSetup, parameters, &innerState, &outerState);
1008         break;
1009       case OutsideIn: {
1010         cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize));
1011         TrackDetMatchInfo result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
1012         cachedTrajectory_.setPropagationStep(currentStepSize);
1013         return result;
1014         break;
1015       }
1016       case Any: {
1017         // check if we deal with clear outside-in case
1018         if (track.innerPosition().Dot(track.innerMomentum()) < 0 &&
1019             track.outerPosition().Dot(track.outerMomentum()) < 0) {
1020           cachedTrajectory_.setPropagationStep(-std::abs(currentStepSize));
1021           TrackDetMatchInfo result;
1022           if (track.innerPosition().Mag2() < track.outerPosition().Mag2())
1023             result = associate(iEvent, iSetup, parameters, &innerState, &outerState);
1024           else
1025             result = associate(iEvent, iSetup, parameters, &outerState, &innerState);
1026           cachedTrajectory_.setPropagationStep(currentStepSize);
1027           return result;
1028         }
1029       }
1030     }
1031   }
1032 
1033   // all other cases
1034   return associate(iEvent, iSetup, parameters, &innerState, &outerState);
1035 }
1036 
1037 TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent,
1038                                                      const edm::EventSetup& iSetup,
1039                                                      const SimTrack& track,
1040                                                      const SimVertex& vertex,
1041                                                      const AssociatorParameters& parameters) {
1042   auto const* bField = &iSetup.getData(parameters.bFieldToken);
1043   return associate(iEvent, iSetup, getFreeTrajectoryState(bField, track, vertex), parameters);
1044 }
1045 
1046 TrackDetMatchInfo TrackDetectorAssociator::associate(const edm::Event& iEvent,
1047                                                      const edm::EventSetup& iSetup,
1048                                                      const GlobalVector& momentum,
1049                                                      const GlobalPoint& vertex,
1050                                                      const int charge,
1051                                                      const AssociatorParameters& parameters) {
1052   auto const* bField = &iSetup.getData(parameters.bFieldToken);
1053   return associate(iEvent, iSetup, getFreeTrajectoryState(bField, momentum, vertex, charge), parameters);
1054 }
1055 
1056 bool TrackDetectorAssociator::crossedIP(const reco::Track& track) {
1057   bool crossed = true;
1058   crossed &= (track.innerPosition().rho() > 3);  // something close to active volume
1059   crossed &= (track.outerPosition().rho() > 3);  // something close to active volume
1060   crossed &=
1061       ((track.innerPosition().x() * track.innerMomentum().x() + track.innerPosition().y() * track.innerMomentum().y() <
1062         0) !=
1063        (track.outerPosition().x() * track.outerMomentum().x() + track.outerPosition().y() * track.outerMomentum().y() <
1064         0));
1065   return crossed;
1066 }
[ 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!