EgammaPhotonAlgos/src/OutInConversionTrackFinder.cc

0001 #include "RecoEgamma/EgammaPhotonAlgos/interface/OutInConversionTrackFinder.h"
0002 //
0003 #include "RecoTracker/CkfPattern/interface/CachingSeedCleanerBySharedInput.h"
0004 #include "RecoTracker/CkfPattern/interface/TransientInitialStateEstimator.h"
0005 //
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 //
0008
0009 #include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"
0010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0011 #include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
0012
0013 //
0014 #include "DataFormats/Common/interface/OwnVector.h"
0015 //
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 //
0018 #include "Utilities/General/interface/precomputed_value_sort.h"
0019
0020 OutInConversionTrackFinder::OutInConversionTrackFinder(const edm::ParameterSet& conf,
0021                                                        const BaseCkfTrajectoryBuilder* trajectoryBuilder,
0022                                                        edm::ConsumesCollector iC)
0023     : ConversionTrackFinder(conf, trajectoryBuilder, iC) {
0024   theTrajectoryCleaner_ = new TrajectoryCleanerBySharedHits(conf);
0025
0026   // get the seed cleaner
0027   std::string cleaner = conf.getParameter<std::string>("OutInRedundantSeedCleaner");
0028   if (cleaner == "CachingSeedCleanerBySharedInput") {
0029     theSeedCleaner_ = new CachingSeedCleanerBySharedInput();
0030   } else if (cleaner == "none") {
0031     theSeedCleaner_ = nullptr;
0032   } else {
0033     throw cms::Exception("OutInRedundantSeedCleaner not found, please use CachingSeedCleanerBySharedInput or none",
0034                          cleaner);
0035   }
0036 }
0037
0038 OutInConversionTrackFinder::~OutInConversionTrackFinder() {
0039   delete theTrajectoryCleaner_;
0040   if (theSeedCleaner_)
0041     delete theSeedCleaner_;
0042 }
0043
0044 std::vector<Trajectory> OutInConversionTrackFinder::tracks(const TrajectorySeedCollection& outInSeeds,
0045                                                            TrackCandidateCollection& output_p) const {
0046   //  std::cout  << "OutInConversionTrackFinder::tracks getting " <<  outInSeeds.size() << " Out-In seeds " << "\n";;
0047
0048   std::vector<Trajectory> tmpO;
0049   tmpO.erase(tmpO.begin(), tmpO.end());
0050
0051   std::vector<Trajectory> result;
0052   result.erase(result.begin(), result.end());
0053
0054   std::vector<Trajectory> rawResult;
0055   if (theSeedCleaner_)
0056     theSeedCleaner_->init(&rawResult);
0057
0058   ///// This loop is only for debugging
0059   /*
0060   for(TrajectorySeedCollection::const_iterator iSeed=outInSeeds.begin(); iSeed!=outInSeeds.end();iSeed++){
0061     DetId tmpId = DetId( iSeed->startingState().detId());
0062     const GeomDet* tmpDet  = theMeasurementTracker_->geomTracker()->idToDet( tmpId );
0063     GlobalVector gv = tmpDet->surface().toGlobal( iSeed->startingState().parameters().momentum() );
0064
0065     LogDebug("OutInConversionTrackFinder") << " OutInConversionTrackFinder::tracks hits in the seed " << iSeed->nHits() << "\n";
0066     LogDebug("OutInConversionTrackFinder")<< " OutInConversionTrackFinder::tracks seed starting state position  " << iSeed->startingState().parameters().position() << " momentum " <<  iSeed->startingState().parameters().momentum() << " charge " << iSeed->startingState().parameters().charge() << " R " << gv.perp() << " eta " << gv.eta() << " phi " << gv.phi() << "\n";
0067
0068     TrajectorySeed::range hitRange = iSeed->recHits();
0069     for (TrajectorySeed::const_iterator ihit = hitRange.first; ihit != hitRange.second; ihit++) {
0070
0071       if ( ihit->isValid() ) {
0072
0073     LogDebug("OutInConversionTrackFinder")  << " Valid hit global position " << theMeasurementTracker_->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()) << " R " << theMeasurementTracker_->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).perp() << " phi " << theMeasurementTracker_->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).phi() << " eta " << theMeasurementTracker_->geomTracker()->idToDet((ihit)->geographicalId())->surface().toGlobal((ihit)->localPosition()).eta() <<    "\n" ;
0074
0075       }
0076     }
0077   }
0078
0079   */
0080
0081   int goodSeed = 0;
0082   std::vector<Trajectory> theTmpTrajectories;
0083   for (TrajectorySeedCollection::const_iterator iSeed = outInSeeds.begin(); iSeed != outInSeeds.end(); iSeed++) {
0084     theTmpTrajectories.clear();
0085
0086     if (!theSeedCleaner_ || theSeedCleaner_->good(&(*iSeed))) {
0087       goodSeed++;
0088
0089       DetId tmpId = DetId(iSeed->startingState().detId());
0090       const GeomDet* tmpDet = theMeasurementTracker_->geomTracker()->idToDet(tmpId);
0091       GlobalVector gv = tmpDet->surface().toGlobal(iSeed->startingState().parameters().momentum());
0092
0093       //      std::cout << " OutInConversionTrackFinder::tracks hits in the seed " << iSeed->nHits() << "\n";
0094       LogDebug("OutInConversionTrackFinder") << " OutInConversionTrackFinder::tracks seed starting state position  "
0095                                              << iSeed->startingState().parameters().position() << " momentum "
0096                                              << iSeed->startingState().parameters().momentum() << " charge "
0097                                              << iSeed->startingState().parameters().charge() << " R " << gv.perp()
0098                                              << " eta " << gv.eta() << " phi " << gv.phi() << "\n";
0099
0100       theCkfTrajectoryBuilder_->trajectories(*iSeed, theTmpTrajectories);
0101
0102       LogDebug("OutInConversionTrackFinder")
0103           << "OutInConversionTrackFinder::track returned " << theTmpTrajectories.size() << " trajectories"
0104           << "\n";
0105
0106       theTrajectoryCleaner_->clean(theTmpTrajectories);
0107
0108       for (std::vector<Trajectory>::const_iterator it = theTmpTrajectories.begin(); it != theTmpTrajectories.end();
0109            it++) {
0110         if (it->isValid()) {
0111           rawResult.push_back(*it);
0112           if (theSeedCleaner_)
0113             theSeedCleaner_->add(&(*it));
0114         }
0115       }
0116     }
0117   }  // end loop over the seeds
0118   LogDebug("OutInConversionTrackFinder") << " OutInConversionTrackFinder::track Good seeds " << goodSeed << "\n";
0119   LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder::track rawResult size after cleaning "
0120                                          << rawResult.size() << "\n";
0121
0122   if (theSeedCleaner_)
0123     theSeedCleaner_->done();
0124
0125   std::vector<Trajectory> unsmoothedResult;
0126   theTrajectoryCleaner_->clean(rawResult);
0127
0128   for (std::vector<Trajectory>::const_iterator itraw = rawResult.begin(); itraw != rawResult.end(); itraw++) {
0129     if ((*itraw).isValid()) {
0130       tmpO.push_back(*itraw);
0131       LogDebug("OutInConversionTrackFinder") << " rawResult num hits " << (*itraw).foundHits() << "\n";
0132     }
0133   }
0134
0135   LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  tmpO size " << tmpO.size()
0136                                          << " before sorting "
0137                                          << "\n";
0138   //  for (std::vector<Trajectory>::const_iterator it =tmpO.begin(); it != tmpO.end(); it++) {
0139   // LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  tmpO num of hits " << (*it).foundHits() << " before ordering " << "\n";
0140   // }
0141
0142   precomputed_value_sort(tmpO.begin(), tmpO.end(), ExtractNumOfHits());
0143
0144   LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  tmpO after sorting "
0145                                          << "\n";
0146   //  for (std::vector<Trajectory>::const_iterator it =tmpO.begin(); it != tmpO.end(); it++) {
0147   // LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  tmpO  num of hits " << (*it).foundHits() << "\n";
0148   // }
0149
0150   for (int i = tmpO.size() - 1; i >= 0; i--) {
0151     unsmoothedResult.push_back(tmpO[i]);
0152   }
0153   LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  unsmoothedResult size  "
0154                                          << unsmoothedResult.size() << "\n";
0155
0156   // for (std::vector<Trajectory>::const_iterator it =  unsmoothedResult.begin(); it !=  unsmoothedResult.end(); it++) {
0157   // LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  unsmoothedResult  after reordering " <<(*it).foundHits() <<  "\n";
0158   //  }
0159
0160   // Check if the inner state is valid
0161   tmpO.clear();
0162   LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  tmpO size " << tmpO.size()
0163                                          << " after clearing "
0164                                          << "\n";
0165   for (std::vector<Trajectory>::const_iterator it = unsmoothedResult.begin(); it != unsmoothedResult.end(); it++) {
0166     if (!it->isValid())
0167       continue;
0168
0169     std::pair<TrajectoryStateOnSurface, const GeomDet*> initState = theInitialState_->innerState(*it);
0170     //  LogDebug("OutInConversionTrackFinder") << " Initial state parameters " << initState.first << "\n";
0171
0172     // temporary protection againt invalid initial states
0173     if (!initState.first.isValid() || initState.second == nullptr) {
0174       LogDebug("OutInConversionTrackFinder") << "invalid innerState, will not make TrackCandidate"
0175                                              << "\n";
0176       ;
0177       continue;
0178     }
0179     tmpO.push_back(*it);
0180   }
0181
0182   LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  tmpO size " << tmpO.size() << " after filling "
0183                                          << "\n";
0184   if (!tmpO.empty()) {
0185     std::vector<Trajectory>::iterator it = tmpO.begin();
0186
0187     // only send out the two best tracks
0188     result.push_back(*it);
0189     if (tmpO.size() > 1)
0190       result.push_back(*(++it));
0191   }
0192
0193   //  for (std::vector<Trajectory>::const_iterator it =result.begin(); it != result.end(); it++) {
0194   // LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  Result  num of hits " << (*it).foundHits() << "\n";
0195   //}
0196
0197   // Convert to TrackCandidates and fill in the output_p
0198   Traj2TrackHits t2t(theCkfTrajectoryBuilder_->hitBuilder(), true);
0199   for (std::vector<Trajectory>::const_iterator it = result.begin(); it != result.end(); it++) {
0200     edm::OwnVector<TrackingRecHit> recHits;
0201     if (it->direction() == alongMomentum)
0202       LogDebug("OutInConversionTrackFinder") << "OutInConv along momentum... " << std::endl;
0203     t2t(*it, recHits, useSplitHits_);
0204
0205     std::pair<TrajectoryStateOnSurface, const GeomDet*> initState = theInitialState_->innerState(*it);
0206
0207     assert(initState.second == recHits.front().det());
0208
0209     // temporary protection againt invalid initial states
0210     if (!initState.first.isValid() || initState.second == nullptr) {
0211       LogDebug("OutInConversionTrackFinder") << "invalid innerState, will not make TrackCandidate" << std::endl;
0212       continue;
0213     }
0214
0215     PTrajectoryStateOnDet state;
0216     if (useSplitHits_ && (initState.second != recHits.front().det()) && recHits.front().det()) {
0217       TrajectoryStateOnSurface propagated =
0218           thePropagator_->propagate(initState.first, recHits.front().det()->surface());
0219       if (!propagated.isValid())
0220         continue;
0221       state = trajectoryStateTransform::persistentState(propagated, recHits.front().rawId());
0222     } else
0223       state = trajectoryStateTransform::persistentState(initState.first, initState.second->geographicalId().rawId());
0224
0225     LogDebug("OutInConversionTrackFinder") << "OutInConversionTrackFinder  Number of hits for the track candidate "
0226                                            << recHits.size() << " TSOS charge " << initState.first.charge() << "\n";
0227
0228     output_p.push_back(TrackCandidate(recHits, it->seed(), state));
0229   }
0230
0231   //  std::cout << "  Returning " << result.size() << "Out In Trajectories  " << "\n";
0232
0233   return result;
0234 }