/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoTracker/TkHitPairs/src/HitPairGeneratorFromLayerPair.cc

Go to the documentation of this file.

#include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"

#include "RecoTracker/TkTrackingRegions/interface/HitRZCompatibility.h"
#include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
#include "RecoTracker/TkTrackingRegions/interface/TrackingRegionBase.h"
#include "RecoTracker/TkHitPairs/interface/OrderedHitPairs.h"
#include "RecoTracker/TkHitPairs/src/InnerDeltaPhi.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "RecoTracker/TkSeedingLayers/interface/SeedingLayer.h"

using namespace GeomDetEnumerators;
using namespace ctfseeding;
using namespace std;

typedef PixelRecoRange<float> Range;

namespace {
  template<class T> inline T sqr( T t) {return t*t;}
}


#include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"
#include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
#include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
#include "FWCore/Framework/interface/ESHandle.h"

HitPairGeneratorFromLayerPair::HitPairGeneratorFromLayerPair(
                                                             const Layer& inner, 
                                                             const Layer& outer, 
                                                             LayerCacheType* layerCache,
                                                             unsigned int nSize,
                                                             unsigned int max)
  : HitPairGenerator(nSize),
    theLayerCache(*layerCache), theOuterLayer(outer), theInnerLayer(inner)
{
  theMaxElement=max;
}


// devirtualizer
#include<tuple>
namespace {

  template<typename Algo>
  struct Kernel {
    using  Base = HitRZCompatibility;
    void set(Base const * a) {
      assert( a->algo()==Algo::me);
      checkRZ=reinterpret_cast<Algo const *>(a);
    }
    
    void operator()(int b, int e, const RecHitsSortedInPhi & innerHitsMap, bool * ok) const {
      constexpr float nSigmaRZ = std::sqrt(12.f);
      for (int i=b; i!=e; ++i) {
        Range allowed = checkRZ->range(innerHitsMap.u[i]);
        float vErr = nSigmaRZ * innerHitsMap.dv[i];
        Range hitRZ(innerHitsMap.v[i]-vErr, innerHitsMap.v[i]+vErr);
        Range crossRange = allowed.intersection(hitRZ);
        ok[i-b] = ! crossRange.empty() ;
      }
    }
    Algo const * checkRZ;
    
  };


  template<typename ... Args> using Kernels = std::tuple<Kernel<Args>...>;

}



void HitPairGeneratorFromLayerPair::hitPairs(
                                             const TrackingRegion & region, OrderedHitPairs & result,
                                             const edm::Event& iEvent, const edm::EventSetup& iSetup) {
  
  auto const & ds = doublets(region,iEvent,iSetup);
  for (std::size_t i=0; i!=ds.size(); ++i) {
    result.push_back( OrderedHitPair( ds.hit(i,HitDoublets::inner),ds.hit(i,HitDoublets::outer) ));
  }
  if (theMaxElement!=0 && result.size() >= theMaxElement){
     result.clear();
    edm::LogError("TooManyPairs")<<"number of pairs exceed maximum, no pairs produced";
  }
}


HitDoublets HitPairGeneratorFromLayerPair::doublets( const TrackingRegion& region, 
                                                    const edm::Event & iEvent, const edm::EventSetup& iSetup) {

  typedef OrderedHitPair::InnerRecHit InnerHit;
  typedef OrderedHitPair::OuterRecHit OuterHit;
  typedef RecHitsSortedInPhi::Hit Hit;

  const RecHitsSortedInPhi & innerHitsMap = theLayerCache(&theInnerLayer, region, iEvent, iSetup);
  if (innerHitsMap.empty()) return HitDoublets(innerHitsMap,innerHitsMap);
 
  const RecHitsSortedInPhi& outerHitsMap = theLayerCache(&theOuterLayer, region, iEvent, iSetup);
  if (outerHitsMap.empty()) return HitDoublets(innerHitsMap,outerHitsMap);

  HitDoublets result(innerHitsMap,outerHitsMap); result.reserve(std::max(innerHitsMap.size(),outerHitsMap.size()));

  InnerDeltaPhi deltaPhi(*theOuterLayer.detLayer(), *theInnerLayer.detLayer(), region, iSetup);

  // std::cout << "layers " << theInnerLayer.detLayer()->seqNum()  << " " << theOuterLayer.detLayer()->seqNum() << std::endl;

  // constexpr float nSigmaRZ = std::sqrt(12.f);
  constexpr float nSigmaPhi = 3.f;
  for (int io = 0; io!=int(outerHitsMap.theHits.size()); ++io) { 
    Hit const & ohit =  outerHitsMap.theHits[io].hit();
    PixelRecoRange<float> phiRange = deltaPhi(outerHitsMap.x[io], 
                                              outerHitsMap.y[io], 
                                              outerHitsMap.z[io], 
                                              nSigmaPhi*outerHitsMap.drphi[io]
                                              );    

    if (phiRange.empty()) continue;

    const HitRZCompatibility *checkRZ = region.checkRZ(theInnerLayer.detLayer(), ohit, iSetup,theOuterLayer.detLayer(), 
                                                       outerHitsMap.rv(io),outerHitsMap.z[io],
                                                       outerHitsMap.isBarrel ? outerHitsMap.du[io] :  outerHitsMap.dv[io],
                                                       outerHitsMap.isBarrel ? outerHitsMap.dv[io] :  outerHitsMap.du[io]
                                                       );
    if(!checkRZ) continue;

    Kernels<HitZCheck,HitRCheck,HitEtaCheck> kernels;

    auto innerRange = innerHitsMap.doubleRange(phiRange.min(), phiRange.max());
    LogDebug("HitPairGeneratorFromLayerPair")<<
      "preparing for combination of: "<< innerRange[1]-innerRange[0]+innerRange[3]-innerRange[2]
                                      <<" inner and: "<< outerHitsMap.theHits.size()<<" outter";
    for(int j=0; j<3; j+=2) {
      auto b = innerRange[j]; auto e=innerRange[j+1];
      bool ok[e-b];
      switch (checkRZ->algo()) {
        case (HitRZCompatibility::zAlgo) :
          std::get<0>(kernels).set(checkRZ);
          std::get<0>(kernels)(b,e,innerHitsMap, ok);
          break;
        case (HitRZCompatibility::rAlgo) :
          std::get<1>(kernels).set(checkRZ);
          std::get<1>(kernels)(b,e,innerHitsMap, ok);
          break;
        case (HitRZCompatibility::etaAlgo) :
          std::get<2>(kernels).set(checkRZ);
          std::get<2>(kernels)(b,e,innerHitsMap, ok);
          break;
      }
      for (int i=0; i!=e-b; ++i) {
        if (!ok[i]) continue;
        if (theMaxElement!=0 && result.size() >= theMaxElement){
          result.clear();
          edm::LogError("TooManyPairs")<<"number of pairs exceed maximum, no pairs produced";
          delete checkRZ;
          return result;
        }
        result.add(b+i,io);
      }
    }
    delete checkRZ;
  }
  LogDebug("HitPairGeneratorFromLayerPair")<<" total number of pairs provided back: "<<result.size();
  return result;
}