NuclearSeedGenerator/src/NuclearInteractionFinder.cc

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   // initialization
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   // Loop on all the RecHits.
0066   while (!NIfound) {
0067     if (it_meas == measurements.end())
0068       break;
0069
0070     // check only the maxHits outermost hits of the primary track
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     // don't check tracks which reach the end of the tracker if checkCompletedTrack==false
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     // Get correct parametrization of the helix of the primary track at the interaction point (to be used by improveCurrentSeed)
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   // This method uses the 3 last TM after the interaction point to calculate the helix parameters
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         // keep one dummy TM at the end, skip the others
0151         result.insert(result.end() - invalidHits, tmp.begin(), tmp.end());
0152       }
0153       invalidHits++;
0154     }
0155   }
0156
0157   // sort the final result, keep dummy measurements at the end
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   // This method returns the seeds calculated by the class SeedsFromNuclearInteraction
0168
0169   const TM& innerTM = tmPairs.first;
0170   const std::vector<TM>& outerTMs = tmPairs.second;
0171
0172   // Loop on all outer TM
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   // loop on all actual seeds
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     // find compatible TM in an outer layer
0209     std::vector<TM> thirdTMs =
0210         findMeasurementsFromTSOS(it_seed->updatedTSOS(), it_seed->outerHitDetId(), layerMeasurements);
0211
0212     // loop on those new TMs
0213     for (std::vector<TM>::const_iterator tm = thirdTMs.begin(); tm != thirdTMs.end(); tm++) {
0214       if (!tm->recHit()->isValid())
0215         continue;
0216
0217       // create new seeds collection using the circle equation
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   // we assume that the error on q/p is equal to 20% of q/p * rescale
0232   mr(0, 0) = (ltp.signedInverseMomentum() * 0.2 * rescale) * (ltp.signedInverseMomentum() * 0.2 * rescale);
0233
0234   // the error on dx/z and dy/dz is fixed to 10% (* rescale)
0235   mr(1, 1) = 1E-2 * rescale * rescale;
0236   mr(2, 2) = 1E-2 * rescale * rescale;
0237
0238   // the error on the local x and y positions are not modified.
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 }