File indexing completed on 2023-01-21 00:19:56
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
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
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
0086 theCaloGeometry_ = &iSetup.getData(parameters.theCaloGeometryToken);
0087
0088
0089 theTrackingGeometry_ = &iSetup.getData(parameters.theTrackingGeometryToken);
0090
0091 if (useDefaultPropagator_ && (!defProp_ || theMagneticFieldWatcher_.check(iSetup))) {
0092
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
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
0135
0136
0137
0138
0139
0140
0141
0142
0143 info.setCaloGeometry(theCaloGeometry_);
0144
0145 cachedTrajectory_.reset_trajectory();
0146
0147
0148
0149
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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
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;
0533 if (abs(track.type()) == 211 ||
0534 abs(track.type()) == 321 ||
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;
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
0563
0564
0565 ROOT::Math::SMatrixIdentity id;
0566 AlgebraicSymMatrix66 covT(id);
0567 covT *= 1e-6;
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
0596
0597
0598
0599
0600
0601
0602 SteppingHelixStateInfo trajectoryPoint = cachedTrajectory_.getStateAtHcal();
0603
0604
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
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
0657 float yOfFirstWire = std::abs(wireTopology->wireAngle()) > 1.E-06f ? -0.5 * length : wireTopology->yOfWire(1);
0658
0659 float yCOWPOffset = yOfFirstWire + 0.5f * length;
0660
0661
0662 float tangent = (wideWidth - narrowWidth) / (2.f * length);
0663
0664 float yPrime = localPoint.y() + std::abs(yOfFirstWire);
0665
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
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
0727
0728
0729
0730
0731
0732
0733
0734
0735
0736
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
0743
0744 for (std::vector<TAMuonChamberMatch>::iterator matchedChamber = matchedChambers.begin();
0745 matchedChamber != matchedChambers.end();
0746 matchedChamber++) {
0747 const GeomDet* geomDet = muonDetIdAssociator_->getGeomDet((*matchedChamber).id);
0748
0749 if (const DTChamber* chamber = dynamic_cast<const DTChamber*>(geomDet)) {
0750
0751 DTRecSegment4DCollection::range range = dtSegments->get(chamber->id());
0752
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
0760 else if (const CSCChamber* chamber = dynamic_cast<const CSCChamber*>(geomDet)) {
0761
0762 CSCSegmentCollection::range range = cscSegments->get(chamber->id());
0763
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
0771 if (parameters.useGEM) {
0772 if (const GEMSuperChamber* chamber = dynamic_cast<const GEMSuperChamber*>(geomDet)) {
0773
0774 GEMSegmentCollection::range range = gemSegments->get(chamber->id());
0775
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
0785 if (parameters.useME0) {
0786 if (const ME0Chamber* chamber = dynamic_cast<const ME0Chamber*>(geomDet)) {
0787
0788 ME0SegmentCollection::range range = me0Segments->get(chamber->id());
0789
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
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
0853
0854
0855
0856 muonSegment.segmentLocalErrorXDxDz = -999.f;
0857 muonSegment.segmentLocalErrorYDyDz = -999.f;
0858 muonSegment.hasZed = true;
0859 muonSegment.hasPhi = true;
0860
0861
0862 muonSegment.t0 = 0.f;
0863 if (dtseg) {
0864 if ((dtseg->hasPhi()) && (!isDTWithoutY)) {
0865 int phiHits = dtseg->phiSegment()->specificRecHits().size();
0866
0867 int hits = 0;
0868 double t0 = 0.;
0869
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
0876
0877
0878
0879
0880
0881 if (hits)
0882 muonSegment.t0 = t0 / hits;
0883
0884 } else {
0885
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
0899
0900 void TrackDetectorAssociator::fillCaloTruth(const edm::Event& iEvent,
0901 TrackDetMatchInfo& info,
0902 const AssociatorParameters& parameters) {
0903
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
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
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
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 ) {
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);
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
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
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
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);
1071 crossed &= (track.outerPosition().rho() > 3);
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 }