Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:31:37 

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