File indexing completed on 2025-01-18 03:42:23
0001 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/TrackerRecHit2D/interface/OmniClusterRef.h"
0004 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h"
0005 #include "DataFormats/TrackerRecHit2D/interface/TkCloner.h"
0006 #include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
0007 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0010 #include "MagneticField/Engine/interface/MagneticField.h"
0011 #include "RecoTracker/SiTrackerMRHTools/interface/MultiRecHitCollector.h"
0012 #include "RecoTracker/SiTrackerMRHTools/interface/SiTrackerMultiRecHitUpdator.h"
0013 #include "RecoTracker/TrackProducer/interface/DAFTrackProducerAlgorithm.h"
0014 #include "RecoTracker/TransientTrackingRecHit/interface/TkClonerImpl.h"
0015 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0016 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0017 #include "TrackingTools/PatternTools/interface/TrajAnnealing.h"
0018 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0019 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0020 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0021 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
0022 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0023 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0024 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0025 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0026
0027 DAFTrackProducerAlgorithm::DAFTrackProducerAlgorithm(const edm::ParameterSet& conf)
0028 : conf_(conf), minHits_(conf.getParameter<int>("MinHits")) {}
0029
0030 void DAFTrackProducerAlgorithm::fillPSetDescription(edm::ParameterSetDescription& desc) { desc.add<int>("MinHits", 3); }
0031
0032 void DAFTrackProducerAlgorithm::runWithCandidate(const TrackingGeometry* theG,
0033 const MagneticField* theMF,
0034 const TrajTrackAssociationCollection& TTmap,
0035 const MeasurementTrackerEvent* measTk,
0036 const TrajectoryFitter* theFitter,
0037 const TransientTrackingRecHitBuilder* builder,
0038 const MultiRecHitCollector* measurementCollector,
0039 const SiTrackerMultiRecHitUpdator* updator,
0040 const reco::BeamSpot& bs,
0041 AlgoProductCollection& algoResults,
0042 TrajAnnealingCollection& trajann,
0043 bool TrajAnnSaving_,
0044 AlgoProductCollection& algoResultsBeforeDAF,
0045 AlgoProductCollection& algoResultsAfterDAF) const {
0046 LogDebug("DAFTrackProducerAlgorithm") << "Size of map: " << TTmap.size() << "\n";
0047
0048 int cont = 0;
0049 int nTracksChanged = 0;
0050
0051 for (TrajTrackAssociationCollection::const_iterator itTTmap = TTmap.begin(); itTTmap != TTmap.end(); itTTmap++) {
0052 const edm::Ref<std::vector<Trajectory> > BeforeDAFTraj = itTTmap->key;
0053 std::vector<TrajectoryMeasurement> BeforeDAFTrajMeas = BeforeDAFTraj->measurements();
0054 const reco::TrackRef BeforeDAFTrack = itTTmap->val;
0055
0056 float ndof = 0;
0057 Trajectory CurrentTraj;
0058
0059 if (BeforeDAFTraj->isValid()) {
0060 LogDebug("DAFTrackProducerAlgorithm") << "The trajectory #" << cont + 1 << " is valid. \n";
0061
0062
0063 std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> hits =
0064 collectHits(*BeforeDAFTraj, measurementCollector, &*measTk);
0065
0066
0067 CurrentTraj = fit(hits, theFitter, *BeforeDAFTraj);
0068
0069
0070 for (std::vector<double>::const_iterator ian = updator->getAnnealingProgram().begin();
0071 ian != updator->getAnnealingProgram().end();
0072 ian++) {
0073 if (CurrentTraj.isValid()) {
0074 LogDebug("DAFTrackProducerAlgorithm") << "Seed direction is " << CurrentTraj.seed().direction()
0075 << ".Traj direction is " << CurrentTraj.direction() << std::endl;
0076
0077
0078 std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> curiterationhits =
0079 updateHits(CurrentTraj, updator, &*measTk, *ian);
0080 if (curiterationhits.first.size() < 3) {
0081 LogDebug("DAFTrackProducerAlgorithm")
0082 << "Rejecting trajectory with " << curiterationhits.first.size() << " hits" << std::endl;
0083 CurrentTraj = Trajectory();
0084 break;
0085 }
0086
0087 CurrentTraj = fit(curiterationhits, theFitter, CurrentTraj);
0088
0089
0090 if (TrajAnnSaving_) {
0091 TrajAnnealing temp(CurrentTraj, *ian);
0092 trajann.push_back(temp);
0093 }
0094
0095 LogDebug("DAFTrackProducerAlgorithm") << "done annealing value " << (*ian);
0096
0097 } else
0098 break;
0099 }
0100
0101 int percOfHitsUnchangedAfterDAF =
0102 (1. * checkHits(*BeforeDAFTraj, CurrentTraj) / (1. * BeforeDAFTrajMeas.size())) * 100.;
0103 LogDebug("DAFTrackProducerAlgorithm")
0104 << "Ended annealing program with " << percOfHitsUnchangedAfterDAF << " unchanged." << std::endl;
0105
0106
0107 ndof = calculateNdof(CurrentTraj);
0108
0109
0110
0111 int goodHits = countingGoodHits(CurrentTraj);
0112
0113 if (goodHits >= minHits_) {
0114 bool ok = buildTrack(CurrentTraj, algoResults, ndof, bs, &BeforeDAFTrack);
0115
0116 if (ok)
0117 cont++;
0118
0119
0120 if ((100. - percOfHitsUnchangedAfterDAF) > 0.) {
0121 bool okBefore = buildTrack(*BeforeDAFTraj, algoResultsBeforeDAF, ndof, bs, &BeforeDAFTrack);
0122 bool okAfter = buildTrack(CurrentTraj, algoResultsAfterDAF, ndof, bs, &BeforeDAFTrack);
0123 if (okBefore && okAfter)
0124 nTracksChanged++;
0125 }
0126 } else {
0127 LogDebug("DAFTrackProducerAlgorithm") << "Rejecting trajectory with " << CurrentTraj.foundHits() << " hits";
0128 }
0129 }
0130 else {
0131 LogDebug("DAFTrackProducerAlgorithm") << "Rejecting empty trajectory" << std::endl;
0132 }
0133 }
0134
0135 LogDebug("DAFTrackProducerAlgorithm") << "Number of Tracks found: " << cont << "\n";
0136 LogDebug("DAFTrackProducerAlgorithm") << "Number of Tracks changed: " << nTracksChanged << "\n";
0137 }
0138
0139 std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> DAFTrackProducerAlgorithm::collectHits(
0140 const Trajectory vtraj,
0141 const MultiRecHitCollector* measurementCollector,
0142 const MeasurementTrackerEvent* measTk) const {
0143 LogDebug("DAFTrackProducerAlgorithm") << "Calling DAFTrackProducerAlgorithm::collectHits";
0144
0145
0146 TransientTrackingRecHit::RecHitContainer hits;
0147 std::vector<TrajectoryMeasurement> collectedmeas = measurementCollector->recHits(vtraj, measTk);
0148
0149
0150
0151 if (collectedmeas.empty())
0152 return std::make_pair(TransientTrackingRecHit::RecHitContainer(), TrajectoryStateOnSurface());
0153
0154 for (std::vector<TrajectoryMeasurement>::const_iterator iter = collectedmeas.begin(); iter != collectedmeas.end();
0155 iter++) {
0156 hits.push_back(iter->recHit());
0157 }
0158
0159
0160
0161
0162
0163
0164
0165 TrajectoryStateOnSurface initialStateFromTrack = collectedmeas.front().predictedState();
0166
0167 LogDebug("DAFTrackProducerAlgorithm")
0168 << "Pair (hits, TSOS) with TSOS predicted(collectedmeas.front().predictedState())";
0169 return std::make_pair(hits, initialStateFromTrack);
0170 }
0171
0172 std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface> DAFTrackProducerAlgorithm::updateHits(
0173 const Trajectory vtraj,
0174 const SiTrackerMultiRecHitUpdator* updator,
0175 const MeasurementTrackerEvent* theMTE,
0176 double annealing) const {
0177 LogDebug("DAFTrackProducerAlgorithm") << "Calling DAFTrackProducerAlgorithm::updateHits";
0178 TransientTrackingRecHit::RecHitContainer hits;
0179 std::vector<TrajectoryMeasurement> vmeas = vtraj.measurements();
0180 std::vector<TrajectoryMeasurement>::reverse_iterator imeas;
0181 unsigned int hitcounter = 1;
0182
0183
0184 for (imeas = vmeas.rbegin(); imeas != vmeas.rend(); imeas++, hitcounter++) {
0185 DetId id = imeas->recHit()->geographicalId();
0186 MeasurementDetWithData measDet = theMTE->idToDet(id);
0187
0188 TrajectoryStateCombiner combiner;
0189 TrajectoryStateOnSurface combtsos;
0190 if (hitcounter == vmeas.size())
0191 combtsos = imeas->predictedState();
0192 else if (hitcounter == 1)
0193 combtsos = imeas->backwardPredictedState();
0194 else
0195 combtsos = combiner(imeas->backwardPredictedState(), imeas->predictedState());
0196
0197 PrintHit(&*imeas->recHit(), combtsos);
0198 if (imeas->recHit()->isValid()) {
0199 TransientTrackingRecHit::RecHitPointer updated = updator->update(imeas->recHit(), combtsos, measDet, annealing);
0200 hits.push_back(updated);
0201 } else {
0202 hits.push_back(imeas->recHit());
0203 }
0204 }
0205
0206 TrajectoryStateOnSurface updatedStateFromTrack = vmeas.back().predictedState();
0207
0208
0209 LogDebug("DAFTrackProducerAlgorithm") << "Pair (hits, TSOS) with TSOS predicted (vmeas.back().predictedState())";
0210
0211 return std::make_pair(hits, updatedStateFromTrack);
0212 }
0213
0214 Trajectory DAFTrackProducerAlgorithm::fit(
0215 const std::pair<TransientTrackingRecHit::RecHitContainer, TrajectoryStateOnSurface>& hits,
0216 const TrajectoryFitter* theFitter,
0217 Trajectory vtraj) const {
0218
0219 Trajectory newVec = theFitter->fitOne(TrajectorySeed({}, {}, vtraj.seed().direction()), hits.first, hits.second);
0220
0221 if (newVec.isValid())
0222 return newVec;
0223 else {
0224 LogDebug("DAFTrackProducerAlgorithm") << "Fit no valid.";
0225 return Trajectory();
0226 }
0227 }
0228
0229 bool DAFTrackProducerAlgorithm::buildTrack(const Trajectory vtraj,
0230 AlgoProductCollection& algoResults,
0231 float ndof,
0232 const reco::BeamSpot& bs,
0233 const reco::TrackRef* BeforeDAFTrack) const {
0234 LogDebug("DAFTrackProducerAlgorithm") << " BUILDER " << std::endl;
0235 ;
0236 TrajectoryStateOnSurface innertsos;
0237
0238 if (vtraj.isValid()) {
0239 std::unique_ptr<Trajectory> theTraj(new Trajectory(vtraj));
0240
0241 if (vtraj.direction() == alongMomentum) {
0242
0243 innertsos = vtraj.firstMeasurement().updatedState();
0244 } else {
0245 innertsos = vtraj.lastMeasurement().updatedState();
0246 }
0247
0248 TSCBLBuilderNoMaterial tscblBuilder;
0249 TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(*(innertsos.freeState()), bs);
0250
0251 if (tscbl.isValid() == false)
0252 return false;
0253
0254 GlobalPoint v = tscbl.trackStateAtPCA().position();
0255 math::XYZPoint pos(v.x(), v.y(), v.z());
0256 GlobalVector p = tscbl.trackStateAtPCA().momentum();
0257 math::XYZVector mom(p.x(), p.y(), p.z());
0258
0259
0260
0261 auto algo = (*BeforeDAFTrack)->algo();
0262 std::unique_ptr<reco::Track> theTrack(new reco::Track(vtraj.chiSquared(),
0263 ndof,
0264 pos,
0265 mom,
0266 tscbl.trackStateAtPCA().charge(),
0267 tscbl.trackStateAtPCA().curvilinearError(),
0268 algo));
0269 theTrack->setQualityMask((*BeforeDAFTrack)->qualityMask());
0270
0271 AlgoProduct aProduct{theTraj.release(), theTrack.release(), vtraj.direction(), 0};
0272 algoResults.push_back(aProduct);
0273
0274 return true;
0275 } else {
0276 LogDebug("DAFTrackProducerAlgorithm") << " BUILDER NOT POSSIBLE: traj is not valid" << std::endl;
0277 ;
0278 return false;
0279 }
0280 }
0281
0282 int DAFTrackProducerAlgorithm::countingGoodHits(const Trajectory traj) const {
0283 int ngoodhits = 0;
0284 std::vector<TrajectoryMeasurement> vtm = traj.measurements();
0285
0286 for (std::vector<TrajectoryMeasurement>::const_iterator tm = vtm.begin(); tm != vtm.end(); tm++) {
0287
0288 if (tm->recHit()->isValid()) {
0289 SiTrackerMultiRecHit const& mHit = dynamic_cast<SiTrackerMultiRecHit const&>(*tm->recHit());
0290 std::vector<const TrackingRecHit*> components = mHit.recHits();
0291
0292 int iComp = 0;
0293
0294 for (std::vector<const TrackingRecHit*>::const_iterator iter = components.begin(); iter != components.end();
0295 iter++, iComp++) {
0296
0297 if (mHit.weight(iComp) > 1e-6) {
0298 ngoodhits++;
0299 iComp++;
0300 break;
0301 }
0302 }
0303 }
0304 }
0305
0306 LogDebug("DAFTrackProducerAlgorithm") << "Original number of valid hits " << traj.foundHits()
0307 << " -> hit with good weight (>1e-6) are " << ngoodhits;
0308 return ngoodhits;
0309 }
0310
0311
0312 void DAFTrackProducerAlgorithm::filter(const TrajectoryFitter* fitter,
0313 std::vector<Trajectory>& input,
0314 int minhits,
0315 std::vector<Trajectory>& output,
0316 const TransientTrackingRecHitBuilder* builder) const {
0317 if (input.empty())
0318 return;
0319
0320 int ngoodhits = 0;
0321 std::vector<TrajectoryMeasurement> vtm = input[0].measurements();
0322 TransientTrackingRecHit::RecHitContainer hits;
0323
0324
0325 for (std::vector<TrajectoryMeasurement>::reverse_iterator tm = vtm.rbegin(); tm != vtm.rend(); tm++) {
0326
0327 if (tm->recHit()->isValid()) {
0328 SiTrackerMultiRecHit const& mHit = dynamic_cast<SiTrackerMultiRecHit const&>(*tm->recHit());
0329 std::vector<const TrackingRecHit*> components = mHit.recHits();
0330 int iComp = 0;
0331 bool isGood = false;
0332 for (std::vector<const TrackingRecHit*>::const_iterator iter = components.begin(); iter != components.end();
0333 iter++, iComp++) {
0334
0335 if (mHit.weight(iComp) > 1e-6) {
0336 ngoodhits++;
0337 iComp++;
0338 isGood = true;
0339 break;
0340 }
0341 }
0342 if (isGood) {
0343 TkClonerImpl hc = static_cast<TkTransientTrackingRecHitBuilder const*>(builder)->cloner();
0344 auto tempHit = hc.makeShared(tm->recHit(), tm->updatedState());
0345 hits.push_back(tempHit);
0346 } else
0347 hits.push_back(std::make_shared<InvalidTrackingRecHit>(*tm->recHit()->det(), TrackingRecHit::missing));
0348 } else {
0349 TkClonerImpl hc = static_cast<TkTransientTrackingRecHitBuilder const*>(builder)->cloner();
0350 auto tempHit = hc.makeShared(tm->recHit(), tm->updatedState());
0351 hits.push_back(tempHit);
0352 }
0353 }
0354
0355 LogDebug("DAFTrackProducerAlgorithm") << "Original number of valid hits " << input[0].foundHits()
0356 << "; after filtering " << ngoodhits;
0357 if (ngoodhits > input[0].foundHits())
0358 edm::LogError("DAFTrackProducerAlgorithm")
0359 << "Something wrong: the number of good hits from DAFTrackProducerAlgorithm::filter " << ngoodhits
0360 << " is higher than the original one " << input[0].foundHits();
0361
0362 if (ngoodhits < minhits)
0363 return;
0364
0365 TrajectoryStateOnSurface curstartingTSOS = input.front().lastMeasurement().updatedState();
0366 LogDebug("DAFTrackProducerAlgorithm") << "starting tsos for final refitting " << curstartingTSOS;
0367
0368
0369 output = fitter->fit(TrajectorySeed({}, {}, input.front().seed().direction()),
0370 hits,
0371 TrajectoryStateWithArbitraryError()(curstartingTSOS));
0372
0373 LogDebug("DAFTrackProducerAlgorithm") << "After filtering " << output.size() << " trajectories";
0374 }
0375
0376 float DAFTrackProducerAlgorithm::calculateNdof(const Trajectory vtraj) const {
0377 if (!vtraj.isValid())
0378 return 0;
0379 float ndof = 0;
0380 const std::vector<TrajectoryMeasurement>& meas = vtraj.measurements();
0381 for (std::vector<TrajectoryMeasurement>::const_iterator iter = meas.begin(); iter != meas.end(); iter++) {
0382 if (iter->recHit()->isValid()) {
0383 SiTrackerMultiRecHit const& mHit = dynamic_cast<SiTrackerMultiRecHit const&>(*iter->recHit());
0384 std::vector<const TrackingRecHit*> components = mHit.recHits();
0385 int iComp = 0;
0386 for (std::vector<const TrackingRecHit*>::const_iterator iter2 = components.begin(); iter2 != components.end();
0387 iter2++, iComp++) {
0388 if ((*iter2)->isValid())
0389 ndof += ((*iter2)->dimension()) * mHit.weight(iComp);
0390 }
0391 }
0392 }
0393
0394 return ndof - 5;
0395 }
0396
0397 int DAFTrackProducerAlgorithm::checkHits(Trajectory iInitTraj, const Trajectory iFinalTraj) const {
0398 std::vector<TrajectoryMeasurement> initmeasurements = iInitTraj.measurements();
0399 std::vector<TrajectoryMeasurement> finalmeasurements = iFinalTraj.measurements();
0400 std::vector<TrajectoryMeasurement>::iterator jmeas;
0401 int nSame = 0;
0402 int ihit = 0;
0403
0404 if (initmeasurements.empty() || finalmeasurements.empty()) {
0405 LogDebug("DAFTrackProducerAlgorithm") << "Initial or Final Trajectory empty.";
0406 return 0;
0407 }
0408
0409 if (initmeasurements.size() != finalmeasurements.size()) {
0410 LogDebug("DAFTrackProducerAlgorithm")
0411 << "Initial Trajectory size(" << initmeasurements.size() << " hits) "
0412 << "is different to final traj size (" << finalmeasurements.size() << ")! No checkHits possible! ";
0413 return 0;
0414 }
0415
0416 for (jmeas = initmeasurements.begin(); jmeas != initmeasurements.end(); jmeas++) {
0417 const TrackingRecHit* initHit = jmeas->recHit()->hit();
0418 if (!initHit->isValid() && ihit == 0)
0419 continue;
0420
0421 if (initHit->isValid()) {
0422 TrajectoryMeasurement imeas = finalmeasurements.at(ihit);
0423 const TrackingRecHit* finalHit = imeas.recHit()->hit();
0424 const TrackingRecHit* MaxWeightHit = nullptr;
0425 float maxweight = 0;
0426
0427 const SiTrackerMultiRecHit* mrh = dynamic_cast<const SiTrackerMultiRecHit*>(finalHit);
0428 if (mrh) {
0429 std::vector<const TrackingRecHit*> components = mrh->recHits();
0430 std::vector<const TrackingRecHit*>::const_iterator icomp;
0431 int hitcounter = 0;
0432
0433 for (icomp = components.begin(); icomp != components.end(); icomp++) {
0434 if ((*icomp)->isValid()) {
0435 double weight = mrh->weight(hitcounter);
0436 if (weight > maxweight) {
0437 MaxWeightHit = *icomp;
0438 maxweight = weight;
0439 }
0440 }
0441 hitcounter++;
0442 }
0443 if (!MaxWeightHit)
0444 continue;
0445
0446 auto myref1 = reinterpret_cast<const BaseTrackerRecHit*>(initHit)->firstClusterRef();
0447 auto myref2 = reinterpret_cast<const BaseTrackerRecHit*>(MaxWeightHit)->firstClusterRef();
0448
0449 if (myref1 == myref2) {
0450 nSame++;
0451 } else {
0452 LogDebug("DAFTrackProducerAlgorithm") << "diverso hit!" << std::endl;
0453 TrajectoryStateOnSurface dummState;
0454 LogTrace("DAFTrackProducerAlgorithm") << " This hit was:\n ";
0455 PrintHit(initHit, dummState);
0456 LogTrace("DAFTrackProducerAlgorithm") << " instead now is:\n ";
0457 PrintHit(MaxWeightHit, dummState);
0458 }
0459 }
0460 } else {
0461 TrajectoryMeasurement imeas = finalmeasurements.at(ihit);
0462 const TrackingRecHit* finalHit = imeas.recHit()->hit();
0463 if (!finalHit->isValid()) {
0464 nSame++;
0465 }
0466 }
0467
0468 ihit++;
0469 }
0470
0471 return nSame;
0472 }
0473
0474 void DAFTrackProducerAlgorithm::PrintHit(const TrackingRecHit* const& hit, TrajectoryStateOnSurface& tsos) const {
0475 if (hit->isValid()) {
0476 LogTrace("DAFTrackProducerAlgorithm")
0477 << " Valid Hit with DetId " << hit->geographicalId().rawId() << " and dim:" << hit->dimension()
0478 << " local position " << hit->localPosition() << " global position " << hit->globalPosition() << " and r "
0479 << hit->globalPosition().perp();
0480 if (tsos.isValid())
0481 LogTrace("DAFTrackProducerAlgorithm") << " TSOS combtsos " << tsos.localPosition();
0482
0483 } else {
0484 LogTrace("DAFTrackProducerAlgorithm") << " Invalid Hit with DetId " << hit->geographicalId().rawId();
0485 }
0486 }