GroupedCkfTrajectoryBuilder.cc

Go to the documentation of this file.

#include "RecoTracker/CkfPattern/interface/GroupedCkfTrajectoryBuilder.h"
#include "TrajectorySegmentBuilder.h"


#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"

#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
#include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
#include "TrackingTools/TrackFitters/interface/KFTrajectoryFitter.h"
#include "RecoTracker/CkfPattern/interface/GroupedTrajCandLess.h"
#include "TrackingTools/TrajectoryFiltering/interface/RegionalTrajectoryFilter.h"
#include "TrackingTools/PatternTools/interface/TempTrajectory.h"
#include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
#include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
#include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
#include "TrackingTools/DetLayers/interface/DetGroup.h"
#include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"

#include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
#include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
#include "TrackingTools/TrajectoryState/interface/BasicSingleTrajectoryState.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"

// only included for RecHit comparison operator:
#include "TrackingTools/TrajectoryCleaning/interface/TrajectoryCleanerBySharedHits.h"

// for looper reconstruction
#include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
#include "TrackingTools/GeomPropagators/interface/HelixBarrelPlaneCrossingByCircle.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"


#include <algorithm> 

using namespace std;

//#define DBG2_GCTB

//#define STANDARD_INTERMEDIARYCLEAN

#ifdef STANDARD_INTERMEDIARYCLEAN
#include "RecoTracker/CkfPattern/interface/IntermediateTrajectoryCleaner.h"
#endif

/* ====== B.M. to be ported layer ===========
#ifdef DBG_GCTB
#include "RecoTracker/CkfPattern/src/ShowCand.h"
#endif
// #define DBG2_GCTB
#ifdef DBG2_GCTB
#include "RecoTracker/CkfPattern/src/SimIdPrinter.h"
#include "Tracker/TkDebugTools/interface/LayerFinderByDet.h"
#include "Tracker/TkLayout/interface/TkLayerName.h"
#endif
=================================== */


GroupedCkfTrajectoryBuilder::
GroupedCkfTrajectoryBuilder(const edm::ParameterSet&              conf,
                const TrajectoryStateUpdator*         updator,
                const Propagator*                     propagatorAlong,
                const Propagator*                     propagatorOpposite,
                const Chi2MeasurementEstimatorBase*   estimator,
                const TransientTrackingRecHitBuilder* recHitBuilder,
                const MeasurementTracker*             measurementTracker,
                const TrajectoryFilter*               filter,
                const TrajectoryFilter*               inOutFilter):


  BaseCkfTrajectoryBuilder(conf,
               updator, propagatorAlong,propagatorOpposite,
               estimator, recHitBuilder, measurementTracker, filter, inOutFilter)
{
  // fill data members from parameters (eventually data members could be dropped)
  //
  theMaxCand                  = conf.getParameter<int>("maxCand");

  theLostHitPenalty           = conf.getParameter<double>("lostHitPenalty");
  theFoundHitBonus            = conf.getParameter<double>("foundHitBonus");
  theIntermediateCleaning     = conf.getParameter<bool>("intermediateCleaning");
  theAlwaysUseInvalid         = conf.getParameter<bool>("alwaysUseInvalidHits");
  theLockHits                 = conf.getParameter<bool>("lockHits");
  theBestHitOnly              = conf.getParameter<bool>("bestHitOnly");
  theMinNrOf2dHitsForRebuild  = 2;
  theRequireSeedHitsInRebuild = conf.getParameter<bool>("requireSeedHitsInRebuild");
  theMinNrOfHitsForRebuild    = max(0,conf.getParameter<int>("minNrOfHitsForRebuild"));
  maxPtForLooperReconstruction     = conf.existsAs<double>("maxPtForLooperReconstruction") ? 
    conf.getParameter<double>("maxPtForLooperReconstruction") : 0;
  maxDPhiForLooperReconstruction     = conf.existsAs<double>("maxDPhiForLooperReconstruction") ? 
    conf.getParameter<double>("maxDPhiForLooperReconstruction") : 2.0;


  /* ======= B.M. to be ported layer ===========
  bool setOK = thePropagator->setMaxDirectionChange(1.6);
  if (!setOK) 
    cout  << "GroupedCkfTrajectoryBuilder WARNING: "
      << "propagator does not support setMaxDirectionChange" 
      << endl;
  //   addStopCondition(theMinPtStopCondition);

  theConfigurableCondition = createAlgo<TrajectoryFilter>(componentConfig("StopCondition"));
  ===================================== */

}

/*
  void GroupedCkfTrajectoryBuilder::setEvent(const edm::Event& event) const
  {
  theMeasurementTracker->update(event);
}
*/

GroupedCkfTrajectoryBuilder::TrajectoryContainer 
GroupedCkfTrajectoryBuilder::trajectories (const TrajectorySeed& seed) const 
{
  TrajectoryContainer ret; 
  ret.reserve(10);
  buildTrajectories(seed, ret, 0);
  return ret; 
}

GroupedCkfTrajectoryBuilder::TrajectoryContainer 
GroupedCkfTrajectoryBuilder::trajectories (const TrajectorySeed& seed, 
                       const TrackingRegion& region) const
{
  TrajectoryContainer ret; 
  ret.reserve(10);
  RegionalTrajectoryFilter regionalCondition(region);
  buildTrajectories(seed, ret, &regionalCondition);
  return ret; 
}

void 
GroupedCkfTrajectoryBuilder::trajectories (const TrajectorySeed& seed, GroupedCkfTrajectoryBuilder::TrajectoryContainer &ret) const 
{
  buildTrajectories(seed,ret,0);
}

void
GroupedCkfTrajectoryBuilder::trajectories (const TrajectorySeed& seed, 
                                            GroupedCkfTrajectoryBuilder::TrajectoryContainer &ret,
                        const TrackingRegion& region) const
{
  RegionalTrajectoryFilter regionalCondition(region);
  buildTrajectories(seed,ret,&regionalCondition);
}

void  
GroupedCkfTrajectoryBuilder::rebuildSeedingRegion(const TrajectorySeed& seed,
                          TrajectoryContainer& result) const {    
  TempTrajectory startingTraj = createStartingTrajectory( seed);
  TempTrajectoryContainer work;

  TrajectoryContainer final;

  work.reserve(result.size());
  for (TrajectoryContainer::iterator traj=result.begin();
       traj!=result.end(); ++traj) {
    if(traj->isValid()) work.push_back(TempTrajectory(*traj));
  }

  rebuildSeedingRegion(startingTraj,work);
  final.reserve(work.size());

  for (TempTrajectoryContainer::iterator traj=work.begin();
       traj!=work.end(); ++traj) {
    final.push_back(traj->toTrajectory());
  }
  
  result.swap(final);
  
}

void
GroupedCkfTrajectoryBuilder::buildTrajectories (const TrajectorySeed& seed,
                                                GroupedCkfTrajectoryBuilder::TrajectoryContainer &result,
                        const TrajectoryFilter* regionalCondition) const
{
  //
  // Build trajectory outwards from seed
  //

  analyseSeed( seed);

  TempTrajectory startingTraj = createStartingTrajectory( seed);

  work_.clear();
  const bool inOut = true;
  groupedLimitedCandidates( startingTraj, regionalCondition, theForwardPropagator, inOut, work_);
  if ( work_.empty() )  return ;



  /*  rebuilding is de-coupled from standard building
  //
  // try to additional hits in the seeding region
  //
  if ( theMinNrOfHitsForRebuild>0 ) {
    // reverse direction
    //thePropagator->setPropagationDirection(oppositeDirection(seed.direction()));
    // rebuild part of the trajectory
    rebuildSeedingRegion(startingTraj,work);
  }

  */

  result.reserve(work_.size());
  for (TempTrajectoryContainer::const_iterator it = work_.begin(), ed = work_.end(); it != ed; ++it) {
      result.push_back( it->toTrajectory() );
  }

  work_.clear(); 
  if (work_.capacity() > work_MaxSize_) {  TempTrajectoryContainer().swap(work_); work_.reserve(work_MaxSize_/2); }

  analyseResult(result);

  LogDebug("CkfPattern")<< "GroupedCkfTrajectoryBuilder: returning result of size " << result.size();

}


void 
GroupedCkfTrajectoryBuilder::groupedLimitedCandidates (TempTrajectory& startingTraj, 
                               const TrajectoryFilter* regionalCondition,
                               const Propagator* propagator, 
                                                       bool inOut,
                               TempTrajectoryContainer& result) const
{
  unsigned int nIter=1;
  TempTrajectoryContainer candidates;
  TempTrajectoryContainer newCand;
  candidates.push_back( startingTraj);

  while ( !candidates.empty()) {

    newCand.clear();
    for (TempTrajectoryContainer::iterator traj=candidates.begin();
     traj!=candidates.end(); traj++) {
      if ( !advanceOneLayer(*traj, regionalCondition, propagator, inOut, newCand, result) ) {
    LogDebug("CkfPattern")<< "GCTB: terminating after advanceOneLayer==false";
    continue;
      }

      LogDebug("CkfPattern")<<"newCand(1): after advanced one layer:\n"<<PrintoutHelper::dumpCandidates(newCand);

      if ((int)newCand.size() > theMaxCand) {
    //ShowCand()(newCand);

    sort( newCand.begin(), newCand.end(), GroupedTrajCandLess(theLostHitPenalty,theFoundHitBonus));
    newCand.erase( newCand.begin()+theMaxCand, newCand.end());
      }
      LogDebug("CkfPattern")<<"newCand(2): after removing extra candidates.\n"<<PrintoutHelper::dumpCandidates(newCand);
    }

    LogDebug("CkfPattern") << "newCand.size() at end = " << newCand.size();
/*
    if (theIntermediateCleaning) {
      candidates.clear();
      candidates = groupedIntermediaryClean(newCand);
    } else {
      candidates.swap(newCand);
    }
*/
    if (theIntermediateCleaning) {
#ifdef STANDARD_INTERMEDIARYCLEAN
    IntermediateTrajectoryCleaner::clean(newCand);  
#else 
    groupedIntermediaryClean(newCand);
#endif  

    }   
    candidates.swap(newCand);

    LogDebug("CkfPattern") <<"candidates(3): "<<result.size()<<" candidates after "<<nIter++<<" groupedCKF iteration: \n"
                   <<PrintoutHelper::dumpCandidates(result)
               <<"\n "<<candidates.size()<<" running candidates are: \n"
               <<PrintoutHelper::dumpCandidates(candidates);
  }
}

std::string whatIsTheNextStep(TempTrajectory& traj , std::pair<TrajectoryStateOnSurface,std::vector<const DetLayer*> >& stateAndLayers){
  std::stringstream buffer;
  vector<const DetLayer*> & nl = stateAndLayers.second;
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
  //B.M. TkLayerName layerName;
  //B.M. buffer << "Started from " << layerName(traj.lastLayer()) 
  const BarrelDetLayer* sbdl = dynamic_cast<const BarrelDetLayer*>(traj.lastLayer());
  const ForwardDetLayer* sfdl = dynamic_cast<const ForwardDetLayer*>(traj.lastLayer());
  if (sbdl) {
    buffer << "Started from " << traj.lastLayer() << " r=" << sbdl->specificSurface().radius() 
       << " phi=" << sbdl->specificSurface().phi() << endl;
  } else if (sfdl) {
    buffer << "Started from " << traj.lastLayer() << " z " << sfdl->specificSurface().position().z()
       << " phi " << sfdl->specificSurface().phi() << endl;
  }
  buffer << "Trying to go to";
  for ( vector<const DetLayer*>::iterator il=nl.begin();
    il!=nl.end(); il++){ 
    //B.M. buffer << " " << layerName(*il)  << " " << *il << endl;
    const BarrelDetLayer* bdl = dynamic_cast<const BarrelDetLayer*>(*il);
    const ForwardDetLayer* fdl = dynamic_cast<const ForwardDetLayer*>(*il);
    
    if (bdl) buffer << " r " << bdl->specificSurface().radius() << endl;
    if (fdl) buffer << " z " << fdl->specificSurface().position().z() << endl;
    //buffer << " " << *il << endl;   
  }
  return buffer.str();
}

std::string whatIsTheStateToUse(TrajectoryStateOnSurface &initial, TrajectoryStateOnSurface & stateToUse, const DetLayer * l){
  std::stringstream buffer;
  buffer << "GCTB: starting from " 
         << " r / phi / z = " << stateToUse.globalPosition().perp()
     << " / " << stateToUse.globalPosition().phi()
     << " / " << stateToUse.globalPosition().z() 
         << " , pt / phi / pz /charge = " 
     << stateToUse.globalMomentum().perp() << " / "  
     << stateToUse.globalMomentum().phi() << " / " 
     << stateToUse.globalMomentum().z() << " / " 
     << stateToUse.charge()
     << " for layer at "<< l << endl;
  buffer << "     errors:";
  for ( int i=0; i<5; i++ )  buffer << " " << sqrt(stateToUse.curvilinearError().matrix()(i,i));
  buffer << endl;
  
  //buffer << "GCTB: starting from r / phi / z = " << initial.globalPosition().perp()
  //<< " / " << initial.globalPosition().phi()
  //<< " / " << initial.globalPosition().z() << " , pt / pz = " 
  //<< initial.globalMomentum().perp() << " / " 
  //<< initial.globalMomentum().z() << " for layer at "
  //<< l << endl;
  //buffer << "     errors:";
  //for ( int i=0; i<5; i++ )  buffer << " " << sqrt(initial.curvilinearError().matrix()(i,i));
  //buffer << endl;
  return buffer.str();
}


bool 
GroupedCkfTrajectoryBuilder::advanceOneLayer (TempTrajectory& traj, 
                          const TrajectoryFilter* regionalCondition, 
                          const Propagator* propagator,
                                              bool inOut,
                          TempTrajectoryContainer& newCand, 
                          TempTrajectoryContainer& result) const
{
  std::pair<TSOS,std::vector<const DetLayer*> > stateAndLayers = findStateAndLayers(traj);
  if(maxPtForLooperReconstruction>0){
    if(
       //stateAndLayers.second.size()==0 &&
       traj.lastLayer()->location()==0 && 
       stateAndLayers.first.globalMomentum().perp()<maxPtForLooperReconstruction &&
       stateAndLayers.first.globalMomentum().perp()>0.3)
      stateAndLayers.second.push_back(traj.lastLayer());
  }

  vector<const DetLayer*>::iterator layerBegin = stateAndLayers.second.begin();
  vector<const DetLayer*>::iterator layerEnd   = stateAndLayers.second.end();

  //   if (nl.empty()) {
  //     addToResult(traj,result,inOut);
  //     return false;
  //   }
  
  LogDebug("CkfPattern")<<whatIsTheNextStep(traj, stateAndLayers);
  
  bool foundSegments(false);
  bool foundNewCandidates(false);
  for ( vector<const DetLayer*>::iterator il=layerBegin; 
    il!=layerEnd; il++) {

    TSOS stateToUse = stateAndLayers.first;

    double dPhiCacheForLoopersReconstruction(0);
    if ((*il)==traj.lastLayer()){
      if(maxPtForLooperReconstruction>0){
    // ------ For loopers reconstruction
    //cout<<" self propagating in advanceOneLayer (for loopers) \n";
    const BarrelDetLayer* sbdl = dynamic_cast<const BarrelDetLayer*>(traj.lastLayer());
    if(sbdl){
      HelixBarrelCylinderCrossing cylinderCrossing(stateToUse.globalPosition(),
                               stateToUse.globalMomentum(),
                               stateToUse.transverseCurvature(),
                               propagator->propagationDirection(),
                               sbdl->specificSurface());
      if(!cylinderCrossing.hasSolution()) continue;
      GlobalPoint starting = stateToUse.globalPosition();
      GlobalPoint target1 = cylinderCrossing.position1();
      GlobalPoint target2 = cylinderCrossing.position2();
      
      GlobalPoint farther = fabs(starting.phi()-target1.phi()) > fabs(starting.phi()-target2.phi()) ?
        target1 : target2;
    
      const Bounds& bounds( sbdl->specificSurface().bounds());
      float length = bounds.length() / 2.f;
    
      /*
        cout << "starting: " << starting << endl;
        cout << "target1: " << target1 << endl;
        cout << "target2: " << target2 << endl;
        cout << "dphi: " << (target1.phi()-target2.phi()) << endl;
        cout << "length: " << length << endl;
      */
    
      /*
      float deltaZ = bounds.thickness()/2.f/fabs(tan(stateToUse.globalDirection().theta()) ) ;
      if(stateToUse.hasError())
        deltaZ += 3*sqrt(stateToUse.cartesianError().position().czz());
      if( fabs(farther.z()) > length + deltaZ ) continue;
      */
      if(fabs(farther.z())*0.95>length) continue;

      Geom::Phi<double> tmpDphi = target1.phi()-target2.phi();
      if(fabs(tmpDphi)>maxDPhiForLooperReconstruction) continue;
      GlobalPoint target((target1.x()+target2.x())/2,
                 (target1.y()+target2.y())/2,
                 (target1.z()+target2.z())/2);
      
      //cout << "target: " << target << endl;
      

      
      TransverseImpactPointExtrapolator extrapolator;
      stateToUse = extrapolator.extrapolate(stateToUse, target, *propagator);
      if (!stateToUse.isValid()) continue; //SK: consider trying the original? probably not

      //dPhiCacheForLoopersReconstruction = fabs(target1.phi()-target2.phi())*2.;
      dPhiCacheForLoopersReconstruction = fabs(tmpDphi);
      traj.incrementLoops();
    }else{
      continue;
    }
      }else{
    // ------ For cosmics reconstruction
    LogDebug("CkfPattern")<<" self propagating in advanceOneLayer.\n from: \n"<<stateToUse;
    //self navigation case
    // go to a middle point first
    TransverseImpactPointExtrapolator middle;
    GlobalPoint center(0,0,0);
    stateToUse = middle.extrapolate(stateToUse, center, *theForwardPropagator);
    
    if (!stateToUse.isValid()) continue;
    LogDebug("CkfPattern")<<"to: "<<stateToUse;
      }
    }
    
    unsigned int maxCandidates = theMaxCand > 21 ? theMaxCand*2 : 42; //limit the number of returned segments
    TrajectorySegmentBuilder layerBuilder(theMeasurementTracker,
                      theLayerMeasurements,
                      **il,*propagator,
                      *theUpdator,*theEstimator,
                      theLockHits,theBestHitOnly, maxCandidates);

    LogDebug("CkfPattern")<<whatIsTheStateToUse(stateAndLayers.first,stateToUse,*il);
    
    TempTrajectoryContainer segments=
      layerBuilder.segments(stateToUse);

    LogDebug("CkfPattern")<< "GCTB: number of segments = " << segments.size();

    if ( !segments.empty() )  foundSegments = true;

    for ( TempTrajectoryContainer::const_iterator is=segments.begin();
      is!=segments.end(); is++ ) {
      //
      // assume "invalid hit only" segment is last in list
      //
      const TempTrajectory::DataContainer & measurements = is->measurements();
      if ( !theAlwaysUseInvalid && is!=segments.begin() && measurements.size()==1 && 
       (measurements.front().recHit()->getType() == TrackingRecHit::missing) )  break;
      //
      // create new candidate
      //
      TempTrajectory newTraj(traj);
      traj.setDPhiCacheForLoopersReconstruction(dPhiCacheForLoopersReconstruction);
      
      //----  avoid to add the same hits more than once in the trajectory ----
      bool toBeRejected(false);
      for(const TempTrajectory::DataContainer::const_iterator revIt = measurements.rbegin(); 
      revIt!=measurements.rend(); --revIt){
    int tmpCounter(0);
    for(const TempTrajectory::DataContainer::const_iterator newTrajMeasIt = newTraj.measurements().rbegin(); 
        newTrajMeasIt != newTraj.measurements().rend(); --newTrajMeasIt){
      //if(tmpCounter==2) break;
      if(revIt->recHit()->geographicalId()==newTrajMeasIt->recHit()->geographicalId()){
        toBeRejected=true;
        break;
      }
      tmpCounter++;
    }
      }

      if(!toBeRejected){
    //newTraj.push(*is);
    //std::cout << "DEBUG: newTraj after push found,lost: " 
    //    << newTraj.foundHits() << " , " 
    //    << newTraj.lostHits() << " , "
    //    << newTraj.measurements().size() << std::endl;
      }else{    
    /*cout << "WARNING: neglect candidate because it contains the same hit twice \n";
    cout << "-- discarded track's pt,eta,#found: " 
         << newTraj.lastMeasurement().updatedState().globalMomentum().perp() << " , "
         << newTraj.lastMeasurement().updatedState().globalMomentum().eta() << " , "
         << newTraj.foundHits() << "\n";
    */
    continue; //Are we sure about this????
      }
      // ------------------------


      newTraj.push(*is);
      //GIO// for ( vector<TM>::const_iterator im=measurements.begin();
      //GIO//        im!=measurements.end(); im++ )  newTraj.push(*im);
      //if ( toBeContinued(newTraj,regionalCondition) ) { TOBE FIXED
      if ( toBeContinued(newTraj, inOut) ) {
    // Have added one more hit to track candidate
    
    LogDebug("CkfPattern")<<"GCTB: adding updated trajectory to candidates: inOut="<<inOut<<" hits="<<newTraj.foundHits();

    newCand.push_back(newTraj);
    foundNewCandidates = true;
      }
      else {
    // Have finished building this track. Check if it passes cuts.

    LogDebug("CkfPattern")<< "GCTB: adding completed trajectory to results if passes cuts: inOut="<<inOut<<" hits="<<newTraj.foundHits();

    addToResult(newTraj, result, inOut);
      }
    }
  }

  if ( !foundSegments ){
    LogDebug("CkfPattern")<< "GCTB: adding input trajectory to result";
    addToResult(traj, result, inOut);
  }
  return foundNewCandidates;
}

//TempTrajectoryContainer
void
GroupedCkfTrajectoryBuilder::groupedIntermediaryClean (TempTrajectoryContainer& theTrajectories) const 
{
  //if (theTrajectories.empty()) return TrajectoryContainer();
  //TrajectoryContainer result;
  if (theTrajectories.empty()) return;  
  //RecHitEqualByChannels recHitEqualByChannels(false, false);
  int firstLayerSize, secondLayerSize;
  vector<const DetLayer*> firstLayers, secondLayers;

  for (TempTrajectoryContainer::iterator firstTraj=theTrajectories.begin();
       firstTraj!=(theTrajectories.end()-1); firstTraj++) {

    if ( (!firstTraj->isValid()) ||
         (!firstTraj->lastMeasurement().recHit()->isValid()) ) continue;
    const TempTrajectory::DataContainer & firstMeasurements = firstTraj->measurements();
    layers(firstMeasurements, firstLayers);
    firstLayerSize = firstLayers.size();
    if ( firstLayerSize<4 )  continue;

    for (TempTrajectoryContainer::iterator secondTraj=(firstTraj+1);
       secondTraj!=theTrajectories.end(); secondTraj++) {

      if ( (!secondTraj->isValid()) ||
           (!secondTraj->lastMeasurement().recHit()->isValid()) ) continue;
      const TempTrajectory::DataContainer & secondMeasurements = secondTraj->measurements();
      layers(secondMeasurements, secondLayers);
      secondLayerSize = secondLayers.size();
      //
      // only candidates using the same last 3 layers are compared
      //
      if ( firstLayerSize!=secondLayerSize )  continue;
      if ( firstLayers[0]!=secondLayers[0] ||
       firstLayers[1]!=secondLayers[1] ||
       firstLayers[2]!=secondLayers[2] )  continue;

      TempTrajectory::DataContainer::const_iterator im1 = firstMeasurements.rbegin();
      TempTrajectory::DataContainer::const_iterator im2 = secondMeasurements.rbegin();
      //
      // check for identical hits in the last layer
      //
      bool unequal(false);
      const DetLayer* layerPtr = firstLayers[0];
      while ( im1!=firstMeasurements.rend()&&im2!=secondMeasurements.rend() ) {
    if ( im1->layer()!=layerPtr || im2->layer()!=layerPtr )  break;
    if ( !(im1->recHit()->isValid()) || !(im2->recHit()->isValid()) ||
         !im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::some) ) {
      unequal = true;
      break;
    }
    --im1;
    --im2;
      }
      if ( im1==firstMeasurements.rend() || im2==secondMeasurements.rend() ||
       im1->layer()==layerPtr || im2->layer()==layerPtr || unequal )  continue;
      //
      // check for invalid hits in the layer -2
      // compare only candidates with invalid / valid combination
      //
      layerPtr = firstLayers[1];
      bool firstValid(true);
      while ( im1!=firstMeasurements.rend()&&im1->layer()==layerPtr ) {
    if ( !im1->recHit()->isValid() )  firstValid = false;
    --im1;
      }
      bool secondValid(true);
      while ( im2!=secondMeasurements.rend()&&im2->layer()==layerPtr ) {
    if ( !im2->recHit()->isValid() )  secondValid = false;
    --im2;
      }
      if ( !tkxor(firstValid,secondValid) )  continue;
      //
      // ask for identical hits in layer -3
      //
      unequal = false;
      layerPtr = firstLayers[2];
      while ( im1!=firstMeasurements.rend()&&im2!=secondMeasurements.rend() ) {
    if ( im1->layer()!=layerPtr || im2->layer()!=layerPtr )  break;
    if ( !(im1->recHit()->isValid()) || !(im2->recHit()->isValid()) ||
         !im1->recHit()->hit()->sharesInput(im2->recHit()->hit(), TrackingRecHit::some) ) {
      unequal = true;
      break;
    }
    --im1;
    --im2;
      }
      if ( im1==firstMeasurements.rend() || im2==secondMeasurements.rend() ||
       im1->layer()==layerPtr || im2->layer()==layerPtr || unequal )  continue;

      if ( !firstValid ) {
    firstTraj->invalidate();
    break;
      }
      else {
    secondTraj->invalidate();
    break;
      }
    }
  }
/*
  for (TempTrajectoryContainer::const_iterator it = theTrajectories.begin();
       it != theTrajectories.end(); it++) {
    if(it->isValid()) result.push_back( *it);
  }

  return result;
*/
  theTrajectories.erase(std::remove_if( theTrajectories.begin(),theTrajectories.end(),
                                        std::not1(std::mem_fun_ref(&TempTrajectory::isValid))),
 //                                     boost::bind(&TempTrajectory::isValid,_1)), 
                        theTrajectories.end());
}

void
GroupedCkfTrajectoryBuilder::layers (const TempTrajectory::DataContainer& measurements,
                                     vector<const DetLayer*> &result) const 
{
  result.clear();

  if ( measurements.empty() )  return ;

  result.push_back(measurements.back().layer());
  TempTrajectory::DataContainer::const_iterator ifirst = measurements.rbegin();
  --ifirst;  
  for ( TempTrajectory::DataContainer::const_iterator im=ifirst;
    im!=measurements.rend(); --im ) {
    if ( im->layer()!=result.back() )  result.push_back(im->layer());
  }

  for (vector<const DetLayer*>::const_iterator iter = result.begin(); iter != result.end(); iter++){
    if (!*iter) edm::LogWarning("CkfPattern")<< "Warning: null det layer!! ";
  }
}

void
GroupedCkfTrajectoryBuilder::rebuildSeedingRegion(TempTrajectory& startingTraj,
                          TempTrajectoryContainer& result) const
{
  //
  // Rebuilding of trajectories. Candidates are taken from result,
  // which will be replaced with the solutions after rebuild
  // (assume vector::swap is more efficient than building new container)
  //
  LogDebug("CkfPattern")<< "Starting to rebuild " << result.size() << " tracks";
  //
  // Fitter (need to create it here since the propagation direction
  // might change between different starting trajectories)
  //
  KFTrajectoryFitter fitter(&(*theBackwardPropagator),&updator(),&estimator());
  //
  TempTrajectoryContainer reFitted;
  TrajectorySeed::range rseedHits = startingTraj.seed().recHits();
  std::vector<const TrackingRecHit*> seedHits;
  //seedHits.insert(seedHits.end(), rseedHits.first, rseedHits.second);
  //for (TrajectorySeed::recHitContainer::const_iterator iter = rseedHits.first; iter != rseedHits.second; iter++){
  //    seedHits.push_back(&*iter);
  //}

  //unsigned int nSeed(seedHits.size());
  unsigned int nSeed(rseedHits.second-rseedHits.first);
  //seedHits.reserve(nSeed);
  TempTrajectoryContainer rebuiltTrajectories;
  for ( TempTrajectoryContainer::iterator it=result.begin();
    it!=result.end(); it++ ) {
    //
    // skip candidates which are not exceeding the seed size 
    // (e.g. because no Tracker layers outside seeding region) 
    //

    if ( it->measurements().size()<=startingTraj.measurements().size() ) {
      rebuiltTrajectories.push_back(*it);
      LogDebug("CkfPattern")<< "RebuildSeedingRegion skipped as in-out trajectory does not exceed seed size.";
      continue;
    }
    //
    // Refit - keep existing trajectory in case fit is not possible
    // or fails
    //

    backwardFit(*it,nSeed,fitter,reFitted,seedHits);
    if ( reFitted.size()!=1 ) {
      rebuiltTrajectories.push_back(*it);
      LogDebug("CkfPattern")<< "RebuildSeedingRegion skipped as backward fit failed";
      //                << "after reFitted.size() " << reFitted.size();
      continue;
    }
    //LogDebug("CkfPattern")<<"after reFitted.size() " << reFitted.size();
    //
    // Rebuild seeding part. In case it fails: keep initial trajectory
    // (better to drop it??)
    //
    int nRebuilt =
      rebuildSeedingRegion (seedHits,reFitted.front(),rebuiltTrajectories);

    if ( nRebuilt==0 ) it->invalidate();  // won't use original in-out track

    if ( nRebuilt<=0 ) rebuiltTrajectories.push_back(*it);

  }
  //
  // Replace input trajectories with new ones
  //
  result.swap(rebuiltTrajectories);
  result.erase(std::remove_if( result.begin(),result.end(),
                   std::not1(std::mem_fun_ref(&TempTrajectory::isValid))),
           result.end());
}

int
GroupedCkfTrajectoryBuilder::rebuildSeedingRegion(const std::vector<const TrackingRecHit*>& seedHits, 
                          TempTrajectory& candidate,
                          TempTrajectoryContainer& result) const 
{
  //
  // Starting from track found by in-out tracking phase, extrapolate it inwards through
  // the seeding region if possible in towards smaller Tracker radii, searching for additional
  // hits.
  // The resulting trajectories are returned in result,
  // the count is the return value.
  //
  TempTrajectoryContainer rebuiltTrajectories;
#ifdef DBG2_GCTB
/*  const LayerFinderByDet layerFinder;
  if ( !seedHits.empty() && seedHits.front().isValid() ) {
    DetLayer* seedLayer = layerFinder(seedHits.front().det());
    cout << "Seed hit at " << seedHits.front().globalPosition()
     << " " << seedLayer << endl;
    cout << "Started from " 
     << candidate.lastMeasurement().updatedState().globalPosition().perp() << " "
     << candidate.lastMeasurement().updatedState().globalPosition().z() << endl;
    pair<bool,TrajectoryStateOnSurface> layerComp(false,TrajectoryStateOnSurface());
    if ( seedLayer ) layerComp =
      seedLayer->compatible(candidate.lastMeasurement().updatedState(),
                  propagator(),estimator());
    pair<bool,TrajectoryStateOnSurface> detComp =
      seedHits.front().det().compatible(candidate.lastMeasurement().updatedState(),
                    propagator(),estimator());
    cout << "  layer compatibility = " << layerComp.first;
    cout << "  det compatibility = " << detComp.first;
    if ( detComp.first ) {
      cout << "  estimate = " 
       << estimator().estimate(detComp.second,seedHits.front()).second ;
    }
    cout << endl;
  }*/
  cout << "Before backward building: #measurements = " 
       << candidate.measurements().size() ; //<< endl;;
#endif
  //
  // Use standard building with standard cuts. Maybe better to use different
  // cuts from "forward" building (e.g. no check on nr. of invalid hits)?
  //
  const bool inOut = false;
  groupedLimitedCandidates(candidate, (const TrajectoryFilter*)0, theBackwardPropagator, inOut, rebuiltTrajectories);

  LogDebug("CkfPattern")<<" After backward building: "<<PrintoutHelper::dumpCandidates(rebuiltTrajectories);

  //
  // Check & count resulting candidates
  //
  int nrOfTrajectories(0);
  bool orig_ok = false;
  //const RecHitEqualByChannels recHitEqual(false,false);
  //vector<TM> oldMeasurements(candidate.measurements());
  for ( TempTrajectoryContainer::iterator it=rebuiltTrajectories.begin();
    it!=rebuiltTrajectories.end(); it++ ) {

    TempTrajectory::DataContainer newMeasurements(it->measurements());
    //
    // Verify presence of seeding hits?
    //
    if ( theRequireSeedHitsInRebuild ) {
      orig_ok = true;
      // no hits found (and possibly some invalid hits discarded): drop track
      if ( newMeasurements.size()<=candidate.measurements().size() ){  
    LogDebug("CkfPattern") << "newMeasurements.size()<=candidate.measurements().size()";
    continue;
      } 
      // verify presence of hits
      //GIO//if ( !verifyHits(newMeasurements.begin()+candidate.measurements().size(),
      //GIO//              newMeasurements.end(),seedHits) ){
      if ( !verifyHits(newMeasurements.rbegin(), 
                       newMeasurements.size() - candidate.measurements().size(),
               seedHits) ){
    LogDebug("CkfPattern")<< "seed hits not found in rebuild";
    continue; 
      }
    }
    //
    // construct final trajectory in the right order
    //
    TempTrajectory reversedTrajectory(it->seed(),it->seed().direction());
    reversedTrajectory.setNLoops(it->nLoops());
    for (TempTrajectory::DataContainer::const_iterator im=newMeasurements.rbegin(), ed = newMeasurements.rend();
      im != ed; --im ) {
      reversedTrajectory.push(*im);
    }
    // save & count result
    result.push_back(reversedTrajectory);
    nrOfTrajectories++;

    LogDebug("CkgPattern")<<"New traj direction = " << reversedTrajectory.direction()<<"\n"
              <<PrintoutHelper::dumpMeasurements(reversedTrajectory.measurements());
  }
  // If nrOfTrajectories = 0 and orig_ok = true, this means that a track was actually found on the
  // out-in step (meeting those requirements) but did not have the seed hits in it.
  // In this case when we return we will go ahead and use the original in-out track.

  // If nrOfTrajectories = 0 and orig_ok = false, this means that the out-in step failed to
  // find any track.  Two cases are a technical failure in fitting the original seed hits or
  // because the track did not meet the out-in criteria (which may be stronger than the out-in
  // criteria).  In this case we will NOT allow the original in-out track to be used.

  if ( (nrOfTrajectories == 0) && orig_ok ) {
    nrOfTrajectories = -1;
  }
  return nrOfTrajectories;
}

void
GroupedCkfTrajectoryBuilder::backwardFit (TempTrajectory& candidate, unsigned int nSeed,
                      const TrajectoryFitter& fitter,
                      TempTrajectoryContainer& fittedTracks,
                      std::vector<const TrackingRecHit*>& remainingHits) const
{
  //
  // clear array of non-fitted hits
  //
  remainingHits.clear();
  fittedTracks.clear();
  //
  // skip candidates which are not exceeding the seed size
  // (e.g. Because no Tracker layers exist outside seeding region)
  //
  if ( candidate.measurements().size()<=nSeed ) {
    fittedTracks.clear();
    return;
  }

  LogDebug("CkfPattern")<<"nSeed " << nSeed << endl
            << "Old traj direction = " << candidate.direction() << endl
            <<PrintoutHelper::dumpMeasurements(candidate.measurements());

  //
  // backward fit trajectory.
  // (Will try to fit only hits outside the seeding region. However,
  // if there are not enough of these, it will also use the seeding hits).
  //
  TempTrajectory::DataContainer oldMeasurements(candidate.measurements());
//   int nOld(oldMeasurements.size());
//   const unsigned int nHitAllMin(5);
//   const unsigned int nHit2dMin(2);
  unsigned int nHit(0);    // number of valid hits after seeding region
  //unsigned int nHit2d(0);  // number of valid hits after seeding region with 2D info
  // use all hits except the first n (from seed), but require minimum
  // specified in configuration.
  //  Swapped over next two lines.
  unsigned int nHitMin = max(candidate.foundHits()-nSeed,theMinNrOfHitsForRebuild);
  //  unsigned int nHitMin = oldMeasurements.size()-nSeed;
  // we want to rebuild only if the number of VALID measurements excluding the seed measurements is higher than the cut
  if (nHitMin<theMinNrOfHitsForRebuild){
    fittedTracks.clear();
        return;
  }

  LogDebug("CkfPattern")/* << "nHitMin " << nHitMin*/ <<"Sizes: " << oldMeasurements.size() << " / ";
  //
  // create input trajectory for backward fit
  //
  Trajectory fwdTraj(candidate.seed(),oppositeDirection(candidate.direction()));
  fwdTraj.setNLoops(candidate.nLoops());
  //const TrajectorySeed seed = TrajectorySeed(PTrajectoryStateOnDet(), TrajectorySeed::recHitContainer(), oppositeDirection(candidate.direction()));
  //Trajectory fwdTraj(seed, oppositeDirection(candidate.direction()));
  std::vector<const DetLayer*> bwdDetLayer; 
  for ( TempTrajectory::DataContainer::const_iterator im=oldMeasurements.rbegin();
    im!=oldMeasurements.rend(); --im) {
    const TrackingRecHit* hit = im->recHit()->hit();
    //
    // add hits until required number is reached
    //
    if ( nHit<nHitMin ){//|| nHit2d<theMinNrOf2dHitsForRebuild ) {
      fwdTraj.push(*im);
      bwdDetLayer.push_back(im->layer());
      //
      // count valid / 2D hits
      //
      if ( hit->isValid() ) {
    nHit++;
    //if ( hit.isMatched() ||
    //     hit.det().detUnits().front()->type().module()==pixel )
        //nHit2d++;
      }
    }
    //if (nHit==nHitMin) lastBwdDetLayer=im->layer();   
    //
    // keep remaining (valid) hits for verification
    //
    else if ( hit->isValid() ) {
      //std::cout << "Adding a remaining hit" << std::endl;
      remainingHits.push_back(hit);
    }
  }
  //
  // Fit only if required number of valid hits can be used
  //
  if ( nHit<nHitMin ){  //|| nHit2d<theMinNrOf2dHitsForRebuild ) {
    fittedTracks.clear();
    return;
  }
  //
  // Do the backward fit (important: start from scaled, not random cov. matrix!)
  //
  TrajectoryStateOnSurface firstTsos(fwdTraj.firstMeasurement().updatedState());
  //cout << "firstTsos "<< firstTsos << endl;
  firstTsos.rescaleError(10.);
  //TrajectoryContainer bwdFitted(fitter.fit(fwdTraj.seed(),fwdTraj.recHits(),firstTsos));
  TrajectoryContainer bwdFitted(fitter.fit(
        TrajectorySeed(PTrajectoryStateOnDet(), TrajectorySeed::recHitContainer(), oppositeDirection(candidate.direction())),
        fwdTraj.recHits(),firstTsos));
  if (bwdFitted.size()){
    LogDebug("CkfPattern")<<"Obtained " << bwdFitted.size() << " bwdFitted trajectories with measurement size " << bwdFitted.front().measurements().size();
    TempTrajectory fitted(fwdTraj.seed(), fwdTraj.direction());
    fitted.setNLoops(fwdTraj.nLoops());
        vector<TM> tmsbf = bwdFitted.front().measurements();
    int iDetLayer=0;
    //this is ugly but the TM in the fitted track do not contain the DetLayer.
    //So we have to cache the detLayer pointers and replug them in.
    //For the backward building it would be enaugh to cache the last DetLayer, 
    //but for the intermediary cleaning we need all
    for ( vector<TM>::const_iterator im=tmsbf.begin();im!=tmsbf.end(); im++ ) {
        fitted.push(TM( (*im).forwardPredictedState(),
                (*im).backwardPredictedState(),
                (*im).updatedState(),
                (*im).recHit(),
                (*im).estimate(),
                bwdDetLayer[iDetLayer]));

        LogDebug("CkfPattern")<<PrintoutHelper::dumpMeasurement(*im);
        iDetLayer++;
    }
/*
    TM lastMeas = bwdFitted.front().lastMeasurement();
    fitted.pop();
    fitted.push(TM(lastMeas.forwardPredictedState(), 
                   lastMeas.backwardPredictedState(), 
                   lastMeas.updatedState(),
                   lastMeas.recHit(),
                   lastMeas.estimate(),
                               lastBwdDetLayer));*/
    fittedTracks.push_back(fitted);
  }
  //
  // save result
  //
  //fittedTracks.swap(bwdFitted);
  //cout << "Obtained " << fittedTracks.size() << " fittedTracks trajectories with measurement size " << fittedTracks.front().measurements().size() << endl;
}

bool
GroupedCkfTrajectoryBuilder::verifyHits (TempTrajectory::DataContainer::const_iterator rbegin,
                                         size_t maxDepth,
                     const std::vector<const TrackingRecHit*>& hits) const
{
  //
  // verify presence of the seeding hits
  //
  LogDebug("CkfPattern")<<"Checking for " << hits.size() << " hits in "
            << maxDepth << " measurements" << endl;

  TempTrajectory::DataContainer::const_iterator rend = rbegin; 
  while (maxDepth > 0) { --maxDepth; --rend; }
  for ( vector<const TrackingRecHit*>::const_iterator ir=hits.begin();
    ir!=hits.end(); ir++ ) {
    // assume that all seeding hits are valid!
    bool foundHit(false);
    for ( TempTrajectory::DataContainer::const_iterator im=rbegin; im!=rend; --im ) {
      if ( im->recHit()->isValid() && (*ir)->sharesInput(im->recHit()->hit(), TrackingRecHit::some) ) {
    foundHit = true;
    break;
      }
    }
    if ( !foundHit )  return false;
  }
  return true;
}