TIDLayer.cc

Go to the documentation of this file.

#include "TIDLayer.h"

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

#include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"

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

#include<array>
#include "DetGroupMerger.h"


//#include "CommonDet/DetLayout/src/DetLessR.h"


using namespace std;

typedef GeometricSearchDet::DetWithState DetWithState;

namespace {

// 2 0 1
/*
class TIDringLess  {
  // z-position ordering of TID rings indices
public:
  bool operator()( const pair<vector<DetGroup> const *,int> & a,const pair<vector<DetGroup>const *,int> & b) const {
    if(a.second==2) {return true;}
    else if(a.second==0) {
      if(b.second==2) return false;
      if(b.second==1) return true;
    }
    return false;    
  };
};
*/

  // groups in correct order: one may be empty
  inline
  void mergeOutward(std::array<vector<DetGroup>,3> & groups, 
            std::vector<DetGroup> & result ) {
    typedef DetGroupMerger Merger;
    Merger::orderAndMergeTwoLevels(std::move(groups[0]),
                   std::move(groups[1]),result,1,1);
    if(!groups[2].empty()) {
      std::vector<DetGroup> tmp;
      tmp.swap(result);
      Merger::orderAndMergeTwoLevels(std::move(tmp),std::move(groups[2]),result,1,1);      
    }
    
  }
  
  inline
  void mergeInward(std::array<vector<DetGroup>,3> & groups,
           std::vector<DetGroup> & result ) {
    typedef DetGroupMerger Merger;
    Merger::orderAndMergeTwoLevels(std::move(groups[2]),
                   std::move(groups[1]),result,1,1);
    if(!groups[0].empty()) {
      std::vector<DetGroup> tmp;
      tmp.swap(result);
      Merger::orderAndMergeTwoLevels(std::move(tmp),std::move(groups[0]),result,1,1);      
    }
    
  }
  
  
  void
  orderAndMergeLevels(const TrajectoryStateOnSurface& tsos,
              const Propagator& prop,
              std::array<vector<DetGroup>,3> & groups,
              std::vector<DetGroup> & result ) {
    
    float zpos = tsos.globalPosition().z();
    if(tsos.globalMomentum().z()*zpos>0){ // momentum points outwards
      if(prop.propagationDirection() == alongMomentum)
    mergeOutward(groups,result);
      else
    mergeInward(groups,result);
    }
    else{ //  momentum points inwards
      if(prop.propagationDirection() == oppositeToMomentum)
    mergeOutward(groups,result);
      else
      mergeInward(groups,result);    
    }  
    
  }
  
}

//hopefully is never called!
const std::vector<const GeometricSearchDet*>& TIDLayer::components() const{
  if( not theComponents) {
    auto temp = std::make_unique<std::vector<const GeometricSearchDet*>>();
    temp->reserve(3);
    for ( auto c: theComps) temp->push_back(c);
    std::vector<const GeometricSearchDet*>* expected = nullptr;
    if(theComponents.compare_exchange_strong(expected,temp.get())) {
      //this thread set the value
      temp.release();
    }
  }

  return *theComponents;
 }


void
TIDLayer::fillRingPars(int i) {
  const BoundDisk& ringDisk = static_cast<const BoundDisk&>(theComps[i]->surface());
  float ringMinZ = std::abs( ringDisk.position().z()) - ringDisk.bounds().thickness()/2.;
  float ringMaxZ = std::abs( ringDisk.position().z()) + ringDisk.bounds().thickness()/2.; 
  ringPars[i].thetaRingMin =  ringDisk.innerRadius()/ ringMaxZ;
  ringPars[i].thetaRingMax =  ringDisk.outerRadius()/ ringMinZ;
  ringPars[i].theRingR=( ringDisk.innerRadius() +
             ringDisk.outerRadius())/2.;

}


TIDLayer::TIDLayer(vector<const TIDRing*>& rings) :
  RingedForwardLayer(true),
  theComponents{nullptr}
{
  //They should be already R-ordered. TO BE CHECKED!!
  //sort( theRings.begin(), theRings.end(), DetLessR());

  if ( rings.size() != 3) throw DetLayerException("Number of rings in TID layer is not equal to 3 !!");
  setSurface( computeDisk( rings ) );

  for(int i=0; i!=3; ++i) {
    theComps[i]=rings[i];
    fillRingPars(i);
    theBasicComps.insert(theBasicComps.end(),   
             (*rings[i]).basicComponents().begin(),
             (*rings[i]).basicComponents().end());
  }

 
  LogDebug("TkDetLayers") << "==== DEBUG TIDLayer =====" ; 
  LogDebug("TkDetLayers") << "r,zed pos  , thickness, innerR, outerR: " 
              << this->position().perp() << " , "
              << this->position().z() << " , "
              << this->specificSurface().bounds().thickness() << " , "
              << this->specificSurface().innerRadius() << " , "
              << this->specificSurface().outerRadius() ;
}


BoundDisk* 
TIDLayer::computeDisk( const vector<const TIDRing*>& rings) const
{
  float theRmin = rings.front()->specificSurface().innerRadius();
  float theRmax = rings.front()->specificSurface().outerRadius();
  float theZmin = rings.front()->position().z() -
    rings.front()->surface().bounds().thickness()/2;
  float theZmax = rings.front()->position().z() +
    rings.front()->surface().bounds().thickness()/2;
  
  for (vector<const TIDRing*>::const_iterator i = rings.begin(); i != rings.end(); i++) {
    float rmin = (**i).specificSurface().innerRadius();
    float rmax = (**i).specificSurface().outerRadius();
    float zmin = (**i).position().z() - (**i).surface().bounds().thickness()/2.;
    float zmax = (**i).position().z() + (**i).surface().bounds().thickness()/2.;
    theRmin = min( theRmin, rmin);
    theRmax = max( theRmax, rmax);
    theZmin = min( theZmin, zmin);
    theZmax = max( theZmax, zmax);
  }
  
  float zPos = (theZmax+theZmin)/2.;
  PositionType pos(0.,0.,zPos);
  RotationType rot;

  return new BoundDisk( pos, rot, new SimpleDiskBounds(theRmin, theRmax,    
                      theZmin-zPos, theZmax-zPos));

}


TIDLayer::~TIDLayer(){
  for (auto c : theComps) delete c;

  delete theComponents.load();
} 

  

void
TIDLayer::groupedCompatibleDetsV( const TrajectoryStateOnSurface& startingState,
                 const Propagator& prop,
                 const MeasurementEstimator& est,
                 std::vector<DetGroup> & result) const
{
  std::array<int,3> const & ringIndices = ringIndicesByCrossingProximity(startingState,prop);
  if ( ringIndices[0]==-1 ) {
    edm::LogError("TkDetLayers") << "TkRingedForwardLayer::groupedCompatibleDets : error in CrossingProximity";
    return;
  }

  std::array<vector<DetGroup>,3> groupsAtRingLevel;
  //order is ring3,ring1,ring2 i.e. 2 0 1
  //                                0 1 2  
#ifdef __INTEL_COMPILER
  const int ringOrder[3]{1,2,0};
#else
  constexpr int ringOrder[3]{1,2,0};
#endif
  auto index = [&ringIndices,& ringOrder](int i) { return ringOrder[ringIndices[i]];};

  auto & closestResult =  groupsAtRingLevel[index(0)];
  theComps[ringIndices[0]]->groupedCompatibleDetsV( startingState, prop, est, closestResult);       
  if ( closestResult.empty() ){
    theComps[ringIndices[1]]->groupedCompatibleDetsV( startingState, prop, est, result); 
    return;
  }

  DetGroupElement closestGel( closestResult.front().front());  
  float rWindow = computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est); 

  if(!overlapInR(closestGel.trajectoryState(),ringIndices[1],rWindow)) {
    result.swap(closestResult);
    return;
  }

 
  auto & nextResult =  groupsAtRingLevel[index(1)];
  theComps[ringIndices[1]]->groupedCompatibleDetsV( startingState, prop, est, nextResult);
  if(nextResult.empty()) {
    result.swap(closestResult);
    return;
  }
 
  if(!overlapInR(closestGel.trajectoryState(),ringIndices[2],rWindow) ) {
    //then merge 2 levels & return 
    orderAndMergeLevels(closestGel.trajectoryState(),prop,groupsAtRingLevel,result);      
    return;
  }

  auto & nextNextResult =  groupsAtRingLevel[index(2)];
  theComps[ringIndices[2]]->groupedCompatibleDetsV( startingState, prop, est, nextNextResult);   
  if(nextNextResult.empty()) {
    // then merge 2 levels and return 
    orderAndMergeLevels(closestGel.trajectoryState(),prop,groupsAtRingLevel,result); // 
    return;
  }

  // merge 3 level and return merged   
  orderAndMergeLevels(closestGel.trajectoryState(),prop,groupsAtRingLevel, result);  
}


std::array<int,3> 
TIDLayer::ringIndicesByCrossingProximity(const TrajectoryStateOnSurface& startingState,
                     const Propagator& prop ) const
{
  typedef HelixForwardPlaneCrossing Crossing; 
  typedef MeasurementEstimator::Local2DVector Local2DVector;

  HelixPlaneCrossing::PositionType startPos( startingState.globalPosition());
  HelixPlaneCrossing::DirectionType startDir( startingState.globalMomentum());
  PropagationDirection propDir( prop.propagationDirection());
  float rho( startingState.transverseCurvature());

  // calculate the crossings with the ring surfaces
  // rings are assumed to be sorted in R !
  
  Crossing myXing(  startPos, startDir, rho, propDir );

  GlobalPoint   ringCrossings[3];
  // vector<GlobalVector>  ringXDirections;

  for (int i = 0; i < 3 ; i++ ) {
    const BoundDisk & theRing  = static_cast<const BoundDisk &>(theComps[i]->surface());
    pair<bool,double> pathlen = myXing.pathLength( theRing);
    if ( pathlen.first ) { 
      ringCrossings[i] = GlobalPoint( myXing.position(pathlen.second ));
      // ringXDirections.push_back( GlobalVector( myXing.direction(pathlen.second )));
    } else {
      // TO FIX.... perhaps there is something smarter to do
      //throw DetLayerException("trajectory doesn't cross TID rings");
      ringCrossings[i] = GlobalPoint( 0.,0.,0.);
      //  ringXDirections.push_back( GlobalVector( 0.,0.,0.));
    }
  }

  int closestIndex = findClosest(ringCrossings);
  int nextIndex    = findNextIndex(ringCrossings,closestIndex);
  if ( closestIndex<0 || nextIndex<0 )  return std::array<int,3>{{-1,-1,-1}};
  int nextNextIndex = -1;
  for(int i=0; i<3 ; i++){
    if(i!= closestIndex && i!=nextIndex) {
      nextNextIndex = i;
      break;
    }
  }
  
  std::array<int,3> indices{{closestIndex,nextIndex,nextNextIndex}};
  return indices;
}




float 
TIDLayer::computeWindowSize( const GeomDet* det, 
                 const TrajectoryStateOnSurface& tsos, 
                 const MeasurementEstimator& est) const
{
  const Plane& startPlane = det->surface();  
  MeasurementEstimator::Local2DVector maxDistance = 
    est.maximalLocalDisplacement( tsos, startPlane);
  return maxDistance.y();
}


int
TIDLayer::findClosest(const GlobalPoint ringCrossing[3] ) const
{
  int theBin = 0;
  float initialR =  ringPars[0].theRingR;
  float rDiff = std::abs( ringCrossing[0].perp() - initialR);
  for (int i = 1; i < 3 ; i++){
    float ringR =  ringPars[i].theRingR;
    float testDiff = std::abs( ringCrossing[i].perp() - ringR);
    if ( testDiff<rDiff ) {
      rDiff = testDiff;
      theBin = i;
    }
  }
  return theBin;
}

int
TIDLayer::findNextIndex(const GlobalPoint ringCrossing[3], int closest ) const
{

  int firstIndexToCheck = (closest != 0)? 0 : 1; 
  float initialR =  ringPars[firstIndexToCheck].theRingR;        
  float rDiff = std::abs( ringCrossing[firstIndexToCheck].perp() - initialR);
  int theBin = firstIndexToCheck;
  for (int i = firstIndexToCheck+1; i < 3 ; i++){
    if ( i != closest) {
      float ringR =  ringPars[i].theRingR;
      float testDiff = std::abs( ringCrossing[i].perp() - ringR);
      if ( testDiff<rDiff ) {
    rDiff = testDiff;
    theBin = i;
      }
    }
  }
  return theBin;
}



bool
TIDLayer::overlapInR( const TrajectoryStateOnSurface& tsos, int index, double ymax ) const 
{
  // assume "fixed theta window", i.e. margin in local y = r is changing linearly with z
  float tsRadius = tsos.globalPosition().perp();
  float thetamin = ( max(0.,tsRadius-ymax))/(std::abs(tsos.globalPosition().z())+10.f); // add 10 cm contingency 
  float thetamax = ( tsRadius + ymax)/(std::abs(tsos.globalPosition().z())-10.f);
  
  // do the theta regions overlap ?

  return !( thetamin > ringPars[index].thetaRingMax || ringPars[index].thetaRingMin > thetamax);
}