Phase2EndcapRing.cc

Go to the documentation of this file.

#include "Phase2EndcapRing.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "TrackingTools/DetLayers/interface/DetLayerException.h"
#include "TrackingTools/DetLayers/interface/MeasurementEstimator.h"
#include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
#include "TrackingTools/DetLayers/interface/rangesIntersect.h"
#include "TrackingTools/DetLayers/interface/ForwardRingDiskBuilderFromDet.h"

#include "LayerCrossingSide.h"
#include "DetGroupMerger.h"
#include "CompatibleDetToGroupAdder.h"

#include "TkDetUtil.h"
#include "DataFormats/GeometryVector/interface/VectorUtil.h"
#include <boost/function.hpp>

using namespace std;

typedef GeometricSearchDet::DetWithState DetWithState;

class DetGroupElementZLess {
public:
  bool operator()(DetGroup a,DetGroup b)
  {
    return (fabs(a.front().det()->position().z()) < fabs(b.front().det()->position().z()));
  }
};

Phase2EndcapRing::Phase2EndcapRing(vector<const GeomDet*>& innerDets,
                   vector<const GeomDet*>& outerDets,
                   const vector<const GeomDet*>& innerDetBrothers,
                   const vector<const GeomDet*>& outerDetBrothers):
  GeometricSearchDet(true),
  theFrontDets(innerDets.begin(),innerDets.end()), 
  theBackDets(outerDets.begin(),outerDets.end()),
  theFrontDetBrothers(innerDetBrothers.begin(),innerDetBrothers.end()), 
  theBackDetBrothers(outerDetBrothers.begin(),outerDetBrothers.end())
{
  theDets.assign(theFrontDets.begin(),theFrontDets.end());
  theDets.insert(theDets.end(),theBackDets.begin(),theBackDets.end());
  theDets.insert(theDets.end(),theFrontDetBrothers.begin(),theFrontDetBrothers.end());
  theDets.insert(theDets.end(),theBackDetBrothers.begin(),theBackDetBrothers.end());


  // the dets should be already phi-ordered. TO BE CHECKED
  //sort( theFrontDets.begin(), theFrontDets.end(), DetLessPhi() );
  //sort( theBackDets.begin(), theBackDets.end(), DetLessPhi() );

  theDisk = ForwardRingDiskBuilderFromDet()( theDets );

  theFrontDisk = ForwardRingDiskBuilderFromDet()( theFrontDets );
  theBackDisk  = ForwardRingDiskBuilderFromDet()( theBackDets );

  theFrontBinFinder = BinFinderType( theFrontDets.front()->surface().position().phi(),
                     theFrontDets.size());
  theBackBinFinder  = BinFinderType( theBackDets.front()->surface().position().phi(),
                     theBackDets.size());  


#ifdef EDM_ML_DEBUG
  LogDebug("TkDetLayers") << "DEBUG INFO for Phase2EndcapRing" ;
  for(vector<const GeomDet*>::const_iterator it=theFrontDets.begin(); 
      it!=theFrontDets.end(); it++){
    LogDebug("TkDetLayers") << "frontDet detId,phi,z,r: "
                            << (*it)->geographicalId().rawId()  << " , "
                << (*it)->surface().position().phi()  << " , "
                << (*it)->surface().position().z()    << " , "
                << (*it)->surface().position().perp() ;
  }

  if(!theFrontDetBrothers.empty()){
    for(vector<const GeomDet*>::const_iterator it=theFrontDetBrothers.begin(); 
        it!=theFrontDetBrothers.end(); it++){
      LogDebug("TkDetLayers") << "frontDet brothers detId,phi,z,r: "
                            << (*it)->geographicalId().rawId()  << " , "
                << (*it)->surface().position().phi()  << " , "
                << (*it)->surface().position().z()    << " , "
                << (*it)->surface().position().perp() ;
    }
  }

  for(vector<const GeomDet*>::const_iterator it=theBackDets.begin(); 
      it!=theBackDets.end(); it++){
    LogDebug("TkDetLayers") << "backDet detId,phi,z,r: "
                            << (*it)->geographicalId().rawId()  << " , "
                << (*it)->surface().position().phi() << " , "
                << (*it)->surface().position().z()   << " , "
                << (*it)->surface().position().perp() ;
  }

  if(!theBackDetBrothers.empty()){
    for(vector<const GeomDet*>::const_iterator it=theBackDetBrothers.begin(); 
        it!=theBackDetBrothers.end(); it++){
      LogDebug("TkDetLayers") << "backDet brothers detId,phi,z,r: "
                            << (*it)->geographicalId().rawId()  << " , "
                << (*it)->surface().position().phi() << " , "
                << (*it)->surface().position().z()   << " , "
                << (*it)->surface().position().perp() ;
    }
  }
#endif

}

Phase2EndcapRing::~Phase2EndcapRing(){

} 

const vector<const GeometricSearchDet*>& 
Phase2EndcapRing::components() const 
{
  throw DetLayerException("Phase2EndcapRing doesn't have GeometricSearchDet components");
}

  
pair<bool, TrajectoryStateOnSurface>
Phase2EndcapRing::compatible( const TrajectoryStateOnSurface&, const Propagator&, 
          const MeasurementEstimator&) const{
  edm::LogError("TkDetLayers") << "temporary dummy implementation of Phase2EndcapRing::compatible()!!" ;
  return pair<bool,TrajectoryStateOnSurface>();
}



void 
Phase2EndcapRing::groupedCompatibleDetsV( const TrajectoryStateOnSurface& tsos,
                 const Propagator& prop,
                 const MeasurementEstimator& est,
                 std::vector<DetGroup>& result) const
{
  SubLayerCrossings  crossings; 
  crossings = computeCrossings( tsos, prop.propagationDirection());
  if(! crossings.isValid()) return;


  std::vector<DetGroup> closestResult;
  std::vector<DetGroup> closestBrotherResult;
  addClosest( tsos, prop, est, crossings.closest(), closestResult,closestBrotherResult); 
  if (closestResult.empty())     return;
  
  DetGroupElement closestGel( closestResult.front().front());  
  int crossingSide = LayerCrossingSide().endcapSide( closestGel.trajectoryState(), prop);
  float phiWindow =  tkDetUtil::computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est); 
  searchNeighbors( tsos, prop, est, crossings.closest(), phiWindow,
           closestResult, closestBrotherResult, false); 

  vector<DetGroup> closestCompleteResult;
  DetGroupMerger::orderAndMergeTwoLevels(std::move(closestResult),std::move(closestBrotherResult),closestCompleteResult,
                     0, crossingSide);

  vector<DetGroup> nextResult;
  vector<DetGroup> nextBrotherResult;
  searchNeighbors( tsos, prop, est, crossings.other(), phiWindow,
           nextResult, nextBrotherResult, true); 

  vector<DetGroup> nextCompleteResult;
  DetGroupMerger::orderAndMergeTwoLevels(std::move(nextResult),std::move(nextBrotherResult),nextCompleteResult,
                     0, crossingSide);

  DetGroupMerger::orderAndMergeTwoLevels( std::move(closestCompleteResult), std::move(nextCompleteResult), result,
                      crossings.closestIndex(), crossingSide);

  //due to propagator problems, when we add single pt sub modules, we should order them in z (endcap)
  if(!theFrontDetBrothers.empty() && !theBackDetBrothers.empty())
    sort(result.begin(),result.end(),DetGroupElementZLess());

#ifdef EDM_ML_DEBUG
  LogTrace("TkDetLayers") <<"Number of groups : " << result.size() << std::endl;
  for (auto&  grp : result) {
    if ( grp.empty() )  continue;
      LogTrace("TkDetLayers") <<"New group in Phase2EndcapRing made by : " << std::endl;
      for (auto const & det : grp) {
          LogTrace("TkDetLayers") <<" geom det at r: " << det.det()->position().perp() <<" id:" << det.det()->geographicalId().rawId()
                    <<" tsos at:" << det.trajectoryState().globalPosition() << std::endl;
      }
    }
#endif

}


SubLayerCrossings 
Phase2EndcapRing::computeCrossings(const TrajectoryStateOnSurface& startingState,
              PropagationDirection propDir) const
{
  auto rho = startingState.transverseCurvature();
  
  HelixPlaneCrossing::PositionType startPos( startingState.globalPosition() );
  HelixPlaneCrossing::DirectionType startDir( startingState.globalMomentum() );
  HelixForwardPlaneCrossing crossing(startPos,startDir,rho,propDir);

  pair<bool,double> frontPath = crossing.pathLength( *theFrontDisk);
  if (!frontPath.first) return SubLayerCrossings();

  pair<bool,double> backPath = crossing.pathLength( *theBackDisk);
  if (!backPath.first) return SubLayerCrossings();

  GlobalPoint gFrontPoint(crossing.position(frontPath.second));
  GlobalPoint gBackPoint( crossing.position(backPath.second));

  int frontIndex = theFrontBinFinder.binIndex(gFrontPoint.barePhi()); 
  SubLayerCrossing frontSLC( 0, frontIndex, gFrontPoint);

  int backIndex = theBackBinFinder.binIndex(gBackPoint.barePhi());
  SubLayerCrossing backSLC( 1, backIndex, gBackPoint);

  
  // 0ss: frontDisk has index=0, backDisk has index=1
  float frontDist = std::abs(Geom::deltaPhi( gFrontPoint.barePhi(), 
                         theFrontDets[frontIndex]->surface().phi()));
  float backDist = std::abs(Geom::deltaPhi( gBackPoint.barePhi(), 
                        theBackDets[backIndex]->surface().phi()));


  if (frontDist < backDist) {
    return SubLayerCrossings( frontSLC, backSLC, 0);
  }
  else {
    return SubLayerCrossings( backSLC, frontSLC, 1);
  } 
}

bool Phase2EndcapRing::addClosest( const TrajectoryStateOnSurface& tsos,
                 const Propagator& prop,
                 const MeasurementEstimator& est,
                 const SubLayerCrossing& crossing,
                 vector<DetGroup>& result,
                 vector<DetGroup>& brotherresult) const
{
  const vector<const GeomDet*>& sub( subLayer( crossing.subLayerIndex()));
  const GeomDet* det(sub[crossing.closestDetIndex()]);
  bool firstgroup = CompatibleDetToGroupAdder::add( *det, tsos, prop, est, result); 
  if(theFrontDetBrothers.empty() && theBackDetBrothers.empty())   return firstgroup;
  // it assumes that the closestDetIndex is ok also for the brother detectors: the crossing is NOT recomputed
  const vector<const GeomDet*>& subBrothers( subLayerBrothers( crossing.subLayerIndex()));
  const GeomDet* detBrother(subBrothers[crossing.closestDetIndex()]);
  bool brothergroup = CompatibleDetToGroupAdder::add( *detBrother, tsos, prop, est, brotherresult); 
  return firstgroup || brothergroup;
}



void Phase2EndcapRing::searchNeighbors( const TrajectoryStateOnSurface& tsos,
                     const Propagator& prop,
                     const MeasurementEstimator& est,
                     const SubLayerCrossing& crossing,
                     float window, 
                     vector<DetGroup>& result,
                     vector<DetGroup>& brotherresult,
                     bool checkClosest) const
{
  const GlobalPoint& gCrossingPos = crossing.position();

  const vector<const GeomDet*>& sLayer( subLayer( crossing.subLayerIndex()));
  // It assumes that what is ok for the front modules in the pt modules is ok also for the back module
  const vector<const GeomDet*>& sBrotherLayer( subLayerBrothers( crossing.subLayerIndex()));
 
  int closestIndex = crossing.closestDetIndex();
  int negStartIndex = closestIndex-1;
  int posStartIndex = closestIndex+1;

  if (checkClosest) { // must decide if the closest is on the neg or pos side
    if ( Geom::phiLess( gCrossingPos.barePhi(), sLayer[closestIndex]->surface().phi())) {
      posStartIndex = closestIndex;
    }
    else {
      negStartIndex = closestIndex;
    }
  }

  const BinFinderType& binFinder = (crossing.subLayerIndex()==0 ? theFrontBinFinder : theBackBinFinder);

  typedef CompatibleDetToGroupAdder Adder;
  int half = sLayer.size()/2;  // to check if dets are called twice....
  for (int idet=negStartIndex; idet >= negStartIndex - half; idet--) {
    const GeomDet & neighborDet = *sLayer[binFinder.binIndex(idet)];
    if (!tkDetUtil::overlapInPhi( gCrossingPos, neighborDet, window)) break;
    if (!Adder::add( neighborDet, tsos, prop, est, result)) break;
    if(theFrontDetBrothers.empty() && theBackDetBrothers.empty()) break;
    // If the two above checks are passed also the brother module will be added with no further checks
    const GeomDet & neighborBrotherDet = *sBrotherLayer[binFinder.binIndex(idet)];
    Adder::add( neighborBrotherDet, tsos, prop, est, brotherresult);
    // maybe also add shallow crossing angle test here???
  }
  for (int idet=posStartIndex; idet < posStartIndex + half; idet++) {
    const GeomDet & neighborDet = *sLayer[binFinder.binIndex(idet)];
    if (!tkDetUtil::overlapInPhi( gCrossingPos, neighborDet, window)) break;
    if (!Adder::add( neighborDet, tsos, prop, est, result)) break;
    if(theFrontDetBrothers.empty() && theBackDetBrothers.empty()) break;
    // If the two above checks are passed also the brother module will be added with no further checks
    const GeomDet & neighborBrotherDet = *sBrotherLayer[binFinder.binIndex(idet)];
    Adder::add( neighborBrotherDet, tsos, prop, est, brotherresult);
    // maybe also add shallow crossing angle test here???
  }
}