61 totSeed=totTraj=totRebuilt=totInvCand=0;
63 void traj(
long long t) {
66 void seed() {++totSeed;}
67 void rebuilt(
long long t) {
72 std::cout <<
"GroupedCkfTrajectoryBuilder stat\nSeed/Traj/Rebuilt " 73 << totSeed <<
'/'<< totTraj <<
'/'<< totRebuilt
76 StatCount() { zero();}
77 ~StatCount() {
print();}
83 void traj(
long long){}
85 void rebuilt(
long long) {}
88 [[cms::thread_safe]] StatCount statCount;
102 #ifdef STANDARD_INTERMEDIARYCLEAN 129 conf.getParameter<
bool>(
"useSameTrajFilter") ?
148 conf.
getParameter<
double>(
"maxPtForLooperReconstruction") : 0;
151 conf.
getParameter<
double>(
"maxDPhiForLooperReconstruction") : 2.0;
232 boost::shared_ptr<const TrajectorySeed> sharedSeed;
235 else sharedSeed = result.front().sharedSeed();
238 work.reserve(result.size());
239 for (
auto && traj : result)
240 if(traj.isValid()) work.emplace_back(
std::move(traj));
249 for (
auto const & it : work)
if (it.isValid()) {
250 final.push_back( it.toTrajectory() );
final.back().setSharedSeed(sharedSeed);
255 statCount.rebuilt(result.size());
262 unsigned int& nCandPerSeed,
266 throw cms::Exception(
"LogicError") <<
"Asking to create trajectories to an un-initialized GroupedCkfTrajectoryBuilder.\nYou have to call clone(const MeasurementTrackerEvent *data) and then call trajectories on it instead.\n";
279 const bool inOut =
true;
281 if ( work_.empty() )
return startingTraj;
285 cleaner.
clean(work_);
287 boost::shared_ptr<const TrajectorySeed> pseed(
new TrajectorySeed(seed));
288 for (
auto const & it : work_)
if (it.isValid()) {
289 result.push_back( it.toTrajectory() ); result.back().setSharedSeed(pseed);
296 LogDebug(
"CkfPattern")<<
"GroupedCkfTrajectoryBuilder: returning result of size " << result.size();
297 statCount.traj(result.size());
301 for (
auto const & traj : result) {
303 for (
auto const & tm : traj.measurements()) {
304 auto const &
hit = tm.recHitR();
305 if (!
hit.isValid()) ++chit[0];
306 if (
hit.det()==
nullptr) ++chit[1];
308 if (
hit.dimension()!=2 ) {
313 if (clus.isPixel()) ++chit[3];
314 else if (thit.isMatched()) {
316 }
else if (thit.isProjected()) {
341 unsigned int nIter=1;
342 unsigned int nCands=0;
343 unsigned int prevNewCandSize=0;
346 candidates.push_back( startingTraj);
348 while ( !candidates.empty()) {
351 for (TempTrajectoryContainer::iterator traj=candidates.begin();
352 traj!=candidates.end(); traj++) {
353 if ( !
advanceOneLayer(seed, *traj, regionalCondition, propagator, inOut, newCand, result) ) {
354 LogDebug(
"CkfPattern")<<
"GCTB: terminating after advanceOneLayer==false";
362 nCands += newCand.size() - prevNewCandSize;
363 prevNewCandSize = newCand.size();
369 newCand.erase( newCand.begin()+
theMaxCand, newCand.end());
374 LogDebug(
"CkfPattern") <<
"newCand.size() at end = " << newCand.size();
384 #ifdef STANDARD_INTERMEDIARYCLEAN 391 candidates.swap(newCand);
393 LogDebug(
"CkfPattern") <<
"candidates(3): "<<result.size()<<
" candidates after "<<nIter++<<
" groupedCKF iteration: \n" 395 <<
"\n "<<candidates.size()<<
" running candidates are: \n" 405 vector<const DetLayer*> & nl = stateAndLayers.second;
419 buffer <<
"Trying to go to";
420 for ( vector<const DetLayer*>::iterator il=nl.begin();
427 if (fdl) buffer <<
" z " << fdl->
specificSurface().position().z() << endl;
435 buffer <<
"GCTB: starting from " 439 <<
" , pt / phi / pz /charge = " 444 <<
" for layer at "<< l << endl;
445 buffer <<
" errors:";
471 std::pair<TSOS,std::vector<const DetLayer*> > && stateAndLayers =
findStateAndLayers(seed,traj);
478 float pt2 = stateAndLayers.first.globalMomentum().perp2();
482 stateAndLayers.second.push_back(traj.
lastLayer());
486 auto layerBegin = stateAndLayers.second.begin();
487 auto layerEnd = stateAndLayers.second.end();
495 LogDebug(
"CkfPattern")<<whatIsTheNextStep(traj, stateAndLayers);
498 bool foundSegments(
false);
499 bool foundNewCandidates(
false);
500 for (
auto il=layerBegin; il!=layerEnd; il++) {
502 TSOS stateToUse = stateAndLayers.first;
504 double dPhiCacheForLoopersReconstruction(0);
517 if(!cylinderCrossing.hasSolution())
continue;
519 GlobalPoint target1 = cylinderCrossing.position1();
520 GlobalPoint target2 = cylinderCrossing.position2();
526 float length = 0.5f*bounds.
length();
542 if(fabs(farther.
z())*0.95
f>length)
continue;
552 stateToUse = extrapolator.
extrapolate(stateToUse, target, *propagator);
553 if (!stateToUse.
isValid())
continue;
556 dPhiCacheForLoopersReconstruction =
std::abs(tmpDphi);
563 LogDebug(
"CkfPattern")<<
" self propagating in advanceOneLayer.\n from: \n"<<stateToUse;
570 if (!stateToUse.
isValid())
continue;
571 LogDebug(
"CkfPattern")<<
"to: "<<stateToUse;
583 LogDebug(
"CkfPattern")<<whatIsTheStateToUse(stateAndLayers.first,stateToUse,*il);
586 auto && segments= layerBuilder.segments(stateToUse);
588 LogDebug(
"CkfPattern")<<
"GCTB: number of segments = " << segments.size();
590 if ( !segments.empty() ) foundSegments =
true;
592 for (
auto is=segments.begin(); is!=segments.end(); is++ ) {
596 auto const & measurements = is->measurements();
602 bool toBeRejected(
false);
603 for(
auto revIt = measurements.rbegin(); revIt!=measurements.rend(); --revIt){
608 if(revIt->recHitR().geographicalId()==newTrajMeasIt->recHitR().geographicalId()
609 && (revIt->recHitR().geographicalId() !=
DetId(0)) ){
620 cout <<
"WARNING: neglect candidate because it contains the same hit twice \n";
621 cout <<
"-- discarded track's pt,eta,#found/lost: " 651 LogDebug(
"CkfPattern")<<
"GCTB: adding updated trajectory to candidates: inOut="<<inOut<<
" hits="<<newTraj.
foundHits();
655 foundNewCandidates =
true;
660 LogDebug(
"CkfPattern")<<
"GCTB: adding completed trajectory to results if passes cuts: inOut="<<inOut<<
" hits="<<newTraj.
foundHits();
666 if ( !foundSegments ){
667 LogDebug(
"CkfPattern")<<
"GCTB: adding input trajectory to result";
668 if (!stateAndLayers.second.empty())
672 return foundNewCandidates;
678 struct LayersInTraj {
681 std::array<DetLayer const *,N>
layers;
688 auto currl = layers[tot] = measurements.
back().
layer();
692 im!=measurements.
rend(); --im ) {
693 if ( im->layer()!=currl ) { ++tot; currl = im->layer();
if (tot<N) layers[tot] = currl;}
712 if (theTrajectories.empty())
return;
714 LayersInTraj
layers[theTrajectories.size()];
716 for (
auto &
t : theTrajectories) {
717 if (
t.isValid() &&
t.lastMeasurement().recHitR().isValid() )
718 layers[ntraj++].
fill(
t);
723 for (
int ifirst=0; ifirst!=ntraj-1; ++ifirst) {
724 auto firstTraj = layers[ifirst].traj;
725 if (!firstTraj->isValid())
continue;
728 int firstLayerSize = layers[ifirst].tot;
729 if ( firstLayerSize<4 )
continue;
730 auto const & firstLayers = layers[ifirst].layers;
732 for (
int isecond= ifirst+1; isecond!=ntraj; ++isecond) {
733 auto secondTraj = layers[isecond].traj;
734 if (!secondTraj->isValid())
continue;
738 int secondLayerSize = layers[isecond].tot;
742 if ( firstLayerSize!=secondLayerSize )
continue;
743 auto const & secondLayers = layers[isecond].layers;
744 if ( firstLayers[0]!=secondLayers[0] ||
745 firstLayers[1]!=secondLayers[1] ||
746 firstLayers[2]!=secondLayers[2] )
continue;
754 const DetLayer* layerPtr = firstLayers[0];
755 while ( im1!=firstMeasurements.
rend()&&im2!=secondMeasurements.
rend() ) {
756 if ( im1->layer()!=layerPtr || im2->layer()!=layerPtr )
break;
757 if ( !(im1->recHit()->isValid()) || !(im2->recHit()->isValid()) ||
766 if ( im1==firstMeasurements.
rend() || im2==secondMeasurements.
rend() ||
767 im1->layer()==layerPtr || im2->layer()==layerPtr || unequal )
continue;
772 layerPtr = firstLayers[1];
774 while ( im1!=firstMeasurements.
rend()&&im1->layer()==layerPtr ) {
775 if ( !im1->recHit()->isValid() ) firstValid =
false;
778 bool secondValid(
true);
779 while ( im2!=secondMeasurements.
rend()&&im2->layer()==layerPtr ) {
780 if ( !im2->recHit()->isValid() ) secondValid =
false;
783 if ( !
tkxor(firstValid,secondValid) )
continue;
788 layerPtr = firstLayers[2];
789 while ( im1!=firstMeasurements.
rend()&&im2!=secondMeasurements.
rend() ) {
790 if ( im1->layer()!=layerPtr || im2->layer()!=layerPtr )
break;
791 if ( !(im1->recHit()->isValid()) || !(im2->recHit()->isValid()) ||
800 if ( im1==firstMeasurements.
rend() || im2==secondMeasurements.
rend() ||
801 im1->layer()==layerPtr || im2->layer()==layerPtr || unequal )
continue;
804 firstTraj->invalidate();
808 secondTraj->invalidate();
821 theTrajectories.erase(std::remove_if( theTrajectories.begin(),theTrajectories.end(),
824 theTrajectories.end());
838 LogDebug(
"CkfPattern")<<
"Starting to rebuild " << result.size() <<
" tracks";
847 std::vector<const TrackingRecHit*> seedHits;
854 unsigned int nSeed(rseedHits.second-rseedHits.first);
858 for ( TempTrajectoryContainer::iterator it=result.begin();
859 it!=result.end(); it++ ) {
865 auto && reFitted =
backwardFit(*it,nSeed,fitter,seedHits);
866 if unlikely( !reFitted.isValid() ) {
867 rebuiltTrajectories.push_back(
std::move(*it));
868 LogDebug(
"CkfPattern")<<
"RebuildSeedingRegion skipped as backward fit failed";
882 for (
size_t i = rebuiltTrajectories.size() - 1;
i < rebuiltTrajectories.size() - nRebuilt - 1; --
i) {
883 rebuiltTrajectories[
i].setStopReason(it->stopReason());
888 if ( nRebuilt<0 ) rebuiltTrajectories.push_back(
std::move(*it));
893 result.swap(rebuiltTrajectories);
894 result.erase(std::remove_if( result.begin(),result.end(),
901 const std::vector<const TrackingRecHit*>& seedHits,
937 cout <<
"Before backward building: #measurements = " 944 const bool inOut =
false;
952 int nrOfTrajectories(0);
953 bool orig_ok =
false;
956 for ( TempTrajectoryContainer::iterator it=rebuiltTrajectories.begin();
957 it!=rebuiltTrajectories.end(); it++ ) {
967 LogDebug(
"CkfPattern") <<
"newMeasurements.size()<=candidate.measurements().size()";
976 LogDebug(
"CkfPattern")<<
"seed hits not found in rebuild";
987 reversedTrajectory.
setNLoops(it->nLoops());
990 reversedTrajectory.
push(*im);
993 LogDebug(
"CkgPattern")<<
"New traj direction = " << reversedTrajectory.
direction()<<
"\n" 1006 if ( (nrOfTrajectories == 0) && orig_ok ) {
1007 nrOfTrajectories = -1;
1009 return nrOfTrajectories;
1015 std::vector<const TrackingRecHit*>& remainingHits)
const 1020 remainingHits.clear();
1022 LogDebug(
"CkfPattern")<<
"nSeed " << nSeed << endl
1023 <<
"Old traj direction = " << candidate.
direction() << endl
1033 unsigned int nHit(0);
1059 if ( nHit<nHitMin ){
1061 bwdDetLayer[nl++]=tm.layer();
1078 remainingHits.push_back(hit);
1084 if unlikely( nHit<nHitMin )
return TempTrajectory();
1098 LogDebug(
"CkfPattern")<<
"Obtained bwdFitted trajectory with measurement size " << bwdFitted.
measurements().size();
1099 TempTrajectory fitted(fwdTraj.
direction(),nSeed);
1100 fitted.setNLoops(fwdTraj.
nLoops());
1107 for ( vector<TM>::const_iterator im=tmsbf.begin();im!=tmsbf.end(); im++ ) {
1108 fitted.emplace( (*im).forwardPredictedState(),
1109 (*im).backwardPredictedState(),
1110 (*im).updatedState(),
1113 bwdDetLayer[iDetLayer]
1137 const std::vector<const TrackingRecHit*>&
hits)
const 1142 LogDebug(
"CkfPattern")<<
"Checking for " << hits.size() <<
" hits in " 1143 << maxDepth <<
" measurements" << endl;
1146 while (maxDepth > 0) { --
maxDepth; --rend; }
1147 for (
auto ir=hits.begin(); ir!=hits.end(); ir++ ) {
1149 bool foundHit(
false);
1150 for (
auto im=rbegin; im!=rend; --im ) {
1151 if ( im->recHit()->isValid() && (*ir)->sharesInput(im->recHit()->hit(),
TrackingRecHit::some) ) {
1156 if ( !foundHit )
return false;
PropagationDirection direction() const
void setEvent_(const edm::Event &iEvent, const edm::EventSetup &iSetup) override
void rescaleError(double factor)
T getParameter(std::string const &) const
const_iterator rend() const
unsigned int groupedLimitedCandidates(const TrajectorySeed &seed, TempTrajectory const &startingTraj, const TrajectoryFilter *regionalCondition, const Propagator *propagator, bool inOut, TempTrajectoryContainer &result) const
virtual float length() const =0
static std::string dumpCandidates(collection &candidates)
void rebuildSeedingRegion(const TrajectorySeed &, TrajectoryContainer &result) const override
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
bool existsAs(std::string const ¶meterName, bool trackiness=true) const
checks if a parameter exists as a given type
void clean(TempTrajectoryContainer &) const override
void join(TempTrajectory &segment)
TrackCharge charge() const
void rebuildTrajectories(TempTrajectory const &startingTraj, const TrajectorySeed &, TrajectoryContainer &result) const override
bool empty() const
True if trajectory has no measurements.
TempTrajectory buildTrajectories(const TrajectorySeed &seed, TrajectoryContainer &ret, unsigned int &nCandPerSeed, const TrajectoryFilter *) const override
common part of both public trajectory building methods
const Propagator * forwardPropagator(const TrajectorySeed &seed) const
std::string print(const Track &, edm::Verbosity=edm::Concise)
Track print utility.
const CurvilinearTrajectoryError & curvilinearError() const
const DataContainer & measurements() const
virtual Location location() const =0
Which part of the detector (barrel, endcap)
Geom::Phi< T > phi() const
void setNLoops(signed char value)
void setNLoops(signed char value)
const TrajectoryStateUpdator * theUpdator
virtual void analyseSeed(const TrajectorySeed &seed) const
GlobalPoint globalPosition() const
ConstRecHitContainer recHits() const
virtual const BoundCylinder & specificSurface() const final
Extension of the interface.
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is rejected(acceptEvent retursns false).There are 3 relevant cases of types of criteria
float maxDPhiForLooperReconstruction
static std::string dumpMeasurement(const TrajectoryMeasurement &tm)
static std::string dumpMeasurements(const std::vector< TrajectoryMeasurement > &v)
unsigned int theMinNrOf2dHitsForRebuild
signed char nLoops() const
virtual Trajectory fitOne(const Trajectory &traj, fitType type=standard) const =0
const TrajectoryMeasurement & lastMeasurement() const
float maxPt2ForLooperReconstruction
PropagationDirection const & direction() const
PropagationDirection direction() const
void setDPhiCacheForLoopersReconstruction(float dphi)
DataContainer const & measurements() const
const TransientTrackingRecHitBuilder * hitBuilder() const
virtual void analyseResult(const TrajectoryContainer &result) const
virtual PropagationDirection propagationDirection() const final
static PropagationDirection oppositeDirection(PropagationDirection dir)
change of propagation direction
GroupedCkfTrajectoryBuilder(const edm::ParameterSet &conf, edm::ConsumesCollector &iC)
constructor from ParameterSet
bool verifyHits(TempTrajectory::DataContainer::const_iterator rbegin, size_t maxDepth, const std::vector< const TrackingRecHit * > &hits) const
Verifies presence of a RecHits in a range of TrajectoryMeasurements.
void groupedIntermediaryClean(TempTrajectoryContainer &theTrajectories) const
intermediate cleaning in the case of grouped measurements
const Chi2MeasurementEstimatorBase & estimator() const
Abs< T >::type abs(const T &t)
std::pair< const_iterator, const_iterator > range
TempTrajectory backwardFit(TempTrajectory &candidate, unsigned int nSeed, const TrajectoryFitter &fitter, std::vector< const TrackingRecHit * > &remainingHits) const
const DetLayer * layer() const
bool advanceOneLayer(const TrajectorySeed &seed, TempTrajectory &traj, const TrajectoryFilter *regionalCondition, const Propagator *propagator, bool inOut, TempTrajectoryContainer &newCand, TempTrajectoryContainer &result) const
void addToResult(boost::shared_ptr< const TrajectorySeed > const &seed, TempTrajectory &traj, TrajectoryContainer &result, bool inOut=false) const
TrajectoryContainer trajectories(const TrajectorySeed &) const override
set Event for the internal MeasurementTracker data member
const_iterator rbegin() const
StateAndLayers findStateAndLayers(const TrajectorySeed &seed, const TempTrajectory &traj) const
const MeasurementTrackerEvent * theMeasurementTracker
bool theKeepOriginalIfRebuildFails
TrajectoryMeasurement const & firstMeasurement() const
virtual const BoundDisk & specificSurface() const final
unsigned int theMinNrOfHitsForRebuild
std::vector< TempTrajectory > TempTrajectoryContainer
const TrajectoryStateUpdator & updator() const
void moveToResult(TempTrajectory &&traj, TempTrajectoryContainer &result, bool inOut=false) const
std::vector< std::vector< double > > tmp
unsigned int nHits() const
const AlgebraicSymMatrix55 & matrix() const
GlobalVector globalMomentum() const
bool toBeContinued(TempTrajectory &traj, bool inOut=false) const
void setStopReason(StopReason s)
virtual OmniClusterRef const & firstClusterRef() const =0
bool tkxor(bool a, bool b) const
bool theIntermediateCleaning
bool isUndef(TrackingRecHit const &hit)
const Propagator * backwardPropagator(const TrajectorySeed &seed) const
TrajectoryStateOnSurface const & updatedState() const
bool theRequireSeedHitsInRebuild
signed char nLoops() const
const DetLayer * lastLayer() const
Redundant method, returns the layer of lastMeasurement() .
const BasicVectorType & basicVector() const
void push(const TrajectoryMeasurement &tm)
std::vector< Trajectory > TrajectoryContainer
const Chi2MeasurementEstimatorBase * theEstimator
TempTrajectory createStartingTrajectory(const TrajectorySeed &seed) const
void push(const TrajectoryMeasurement &tm)
double transverseCurvature() const