Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-01-21 00:19:56 

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