File indexing completed on 2024-04-06 12:28:27
0001 #include "RecoTracker/NuclearSeedGenerator/interface/NuclearInteractionFinder.h"
0002 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
0003 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0004
0005 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0006
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0011 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0012
0013 NuclearInteractionFinder::NuclearInteractionFinder(const Config& iConfig,
0014 const TrackerGeometry* iTrackerGeom,
0015 const Propagator* iPropagator,
0016 const MeasurementEstimator* iEstimator,
0017 const MeasurementTracker* iMeasurementTracker,
0018 const GeometricSearchTracker* iGeomSearchTracker,
0019 const NavigationSchool* iNavigationSchool)
0020 : maxHits(iConfig.maxHits),
0021 rescaleErrorFactor(iConfig.rescaleErrorFactor),
0022 checkCompletedTrack(iConfig.checkCompletedTrack) {
0023 thePropagator = iPropagator;
0024 theEstimator = iEstimator;
0025 theMeasurementTracker = iMeasurementTracker;
0026 theGeomSearchTracker = iGeomSearchTracker;
0027 theNavigationSchool = iNavigationSchool;
0028
0029 LogDebug("NuclearSeedGenerator") << "New NuclearInteractionFinder instance with parameters : \n"
0030 << "maxHits : " << maxHits << "\n"
0031 << "rescaleErrorFactor : " << rescaleErrorFactor << "\n"
0032 << "checkCompletedTrack : " << checkCompletedTrack << "\n";
0033
0034 nuclTester = std::make_unique<NuclearTester>(maxHits, theEstimator, iTrackerGeom);
0035
0036 currentSeed = std::make_unique<SeedFromNuclearInteraction>(thePropagator, iTrackerGeom, iConfig.ptMin);
0037
0038 thePrimaryHelix = std::make_unique<TangentHelix>();
0039 }
0040
0041
0042 bool NuclearInteractionFinder::run(const Trajectory& traj, const MeasurementTrackerEvent& event) {
0043 if (traj.empty() || !traj.isValid())
0044 return false;
0045
0046 LogDebug("NuclearSeedGenerator") << "Analyzis of a new trajectory with a number of valid hits = " << traj.foundHits();
0047
0048 std::vector<TrajectoryMeasurement> measurements = traj.measurements();
0049
0050
0051 nuclTester->reset(measurements.size());
0052 allSeeds.clear();
0053
0054 LayerMeasurements layerMeasurements(*theMeasurementTracker, event);
0055
0056 if (traj.direction() == alongMomentum) {
0057 std::reverse(measurements.begin(), measurements.end());
0058 }
0059
0060 std::vector<TrajectoryMeasurement>::const_iterator it_meas = measurements.begin();
0061
0062 std::vector<double> ncompatibleHits;
0063 bool NIfound = false;
0064
0065
0066 while (!NIfound) {
0067 if (it_meas == measurements.end())
0068 break;
0069
0070
0071 if (nuclTester->nHitsChecked() > maxHits)
0072 break;
0073
0074 nuclTester->push_back(*it_meas, findCompatibleMeasurements(*it_meas, rescaleErrorFactor, layerMeasurements));
0075
0076 LogDebug("NuclearSeedGenerator") << "Number of compatible meas:" << (nuclTester->back()).size() << "\n"
0077 << "Mean distance between hits :" << nuclTester->meanHitDistance() << "\n"
0078 << "Forward estimate :" << nuclTester->fwdEstimate() << "\n";
0079
0080
0081 if (checkCompletedTrack == false && (nuclTester->compatibleHits()).front() == 0)
0082 break;
0083
0084 if (nuclTester->isNuclearInteraction())
0085 NIfound = true;
0086
0087 ++it_meas;
0088 }
0089
0090 if (NIfound) {
0091 LogDebug("NuclearSeedGenerator") << "NUCLEAR INTERACTION FOUND at index : " << nuclTester->nuclearIndex() << "\n";
0092
0093
0094 definePrimaryHelix(measurements.begin() + nuclTester->nuclearIndex() - 1);
0095
0096 this->fillSeeds(nuclTester->goodTMPair());
0097
0098 return true;
0099 }
0100
0101 return false;
0102 }
0103
0104 void NuclearInteractionFinder::definePrimaryHelix(std::vector<TrajectoryMeasurement>::const_iterator it_meas) {
0105
0106
0107 GlobalPoint pt[3];
0108 for (int i = 0; i < 3; i++) {
0109 pt[i] = (it_meas->updatedState()).globalParameters().position();
0110 it_meas++;
0111 }
0112 thePrimaryHelix = std::make_unique<TangentHelix>(pt[0], pt[1], pt[2]);
0113 }
0114
0115 std::vector<TrajectoryMeasurement> NuclearInteractionFinder::findCompatibleMeasurements(
0116 const TM& lastMeas, double rescale, const LayerMeasurements& layerMeasurements) const {
0117 const TSOS& currentState = lastMeas.updatedState();
0118 LogDebug("NuclearSeedGenerator") << "currentState :" << currentState << "\n";
0119
0120 TSOS newState = rescaleError(rescale, currentState);
0121 return findMeasurementsFromTSOS(newState, lastMeas.recHit()->geographicalId(), layerMeasurements);
0122 }
0123
0124 std::vector<TrajectoryMeasurement> NuclearInteractionFinder::findMeasurementsFromTSOS(
0125 const TSOS& currentState, DetId detid, const LayerMeasurements& layerMeasurements) const {
0126 using namespace std;
0127 int invalidHits = 0;
0128 vector<TM> result;
0129 const DetLayer* lastLayer = theGeomSearchTracker->detLayer(detid);
0130 vector<const DetLayer*> nl;
0131
0132 if (lastLayer) {
0133 nl = theNavigationSchool->nextLayers(*lastLayer, *currentState.freeState(), alongMomentum);
0134 } else {
0135 edm::LogError("NuclearInteractionFinder") << "In findCompatibleMeasurements : lastLayer not accessible";
0136 return result;
0137 }
0138
0139 if (nl.empty()) {
0140 LogDebug("NuclearSeedGenerator") << "In findCompatibleMeasurements : no compatible layer found";
0141 return result;
0142 }
0143
0144 for (vector<const DetLayer*>::iterator il = nl.begin(); il != nl.end(); il++) {
0145 vector<TM> tmp = layerMeasurements.measurements((**il), currentState, *thePropagator, *theEstimator);
0146 if (!tmp.empty()) {
0147 if (result.empty())
0148 result = tmp;
0149 else {
0150
0151 result.insert(result.end() - invalidHits, tmp.begin(), tmp.end());
0152 }
0153 invalidHits++;
0154 }
0155 }
0156
0157
0158 if (result.size() > 1) {
0159 sort(result.begin(), result.end() - invalidHits, TrajMeasLessEstim());
0160 }
0161 return result;
0162 }
0163
0164
0165 void NuclearInteractionFinder::fillSeeds(
0166 const std::pair<TrajectoryMeasurement, std::vector<TrajectoryMeasurement> >& tmPairs) {
0167
0168
0169 const TM& innerTM = tmPairs.first;
0170 const std::vector<TM>& outerTMs = tmPairs.second;
0171
0172
0173 for (std::vector<TM>::const_iterator outtm = outerTMs.begin(); outtm != outerTMs.end(); outtm++) {
0174 if ((innerTM.recHit())->isValid() && (outtm->recHit())->isValid()) {
0175 currentSeed->setMeasurements(innerTM.updatedState(), innerTM.recHit(), outtm->recHit());
0176 allSeeds.push_back(*currentSeed);
0177 } else
0178 LogDebug("NuclearSeedGenerator") << "The initial hits for seeding are invalid"
0179 << "\n";
0180 }
0181 return;
0182 }
0183
0184 std::unique_ptr<TrajectorySeedCollection> NuclearInteractionFinder::getPersistentSeeds() {
0185 auto output = std::make_unique<TrajectorySeedCollection>();
0186 for (std::vector<SeedFromNuclearInteraction>::const_iterator it_seed = allSeeds.begin(); it_seed != allSeeds.end();
0187 it_seed++) {
0188 if (it_seed->isValid()) {
0189 output->push_back(it_seed->TrajSeed());
0190 } else
0191 LogDebug("NuclearSeedGenerator") << "The seed is invalid"
0192 << "\n";
0193 }
0194 return output;
0195 }
0196
0197 void NuclearInteractionFinder::improveSeeds(const MeasurementTrackerEvent& event) {
0198 std::vector<SeedFromNuclearInteraction> newSeedCollection;
0199
0200 LayerMeasurements layerMeasurements(*theMeasurementTracker, event);
0201
0202
0203 for (std::vector<SeedFromNuclearInteraction>::const_iterator it_seed = allSeeds.begin(); it_seed != allSeeds.end();
0204 it_seed++) {
0205 if (!it_seed->isValid())
0206 continue;
0207
0208
0209 std::vector<TM> thirdTMs =
0210 findMeasurementsFromTSOS(it_seed->updatedTSOS(), it_seed->outerHitDetId(), layerMeasurements);
0211
0212
0213 for (std::vector<TM>::const_iterator tm = thirdTMs.begin(); tm != thirdTMs.end(); tm++) {
0214 if (!tm->recHit()->isValid())
0215 continue;
0216
0217
0218 currentSeed->setMeasurements(*thePrimaryHelix, it_seed->initialTSOS(), it_seed->outerHit(), tm->recHit());
0219 newSeedCollection.push_back(*currentSeed);
0220 }
0221 }
0222 allSeeds.clear();
0223 allSeeds = newSeedCollection;
0224 }
0225
0226 TrajectoryStateOnSurface NuclearInteractionFinder::rescaleError(float rescale, const TSOS& state) const {
0227 AlgebraicSymMatrix55 m(state.localError().matrix());
0228 AlgebraicSymMatrix55 mr;
0229 LocalTrajectoryParameters ltp = state.localParameters();
0230
0231
0232 mr(0, 0) = (ltp.signedInverseMomentum() * 0.2 * rescale) * (ltp.signedInverseMomentum() * 0.2 * rescale);
0233
0234
0235 mr(1, 1) = 1E-2 * rescale * rescale;
0236 mr(2, 2) = 1E-2 * rescale * rescale;
0237
0238
0239 mr(3, 3) = m(3, 3);
0240 mr(4, 4) = m(4, 4);
0241
0242 return TSOS(ltp, mr, state.surface(), &(state.globalParameters().magneticField()), state.surfaceSide());
0243 }