/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_5_3_13_patch3/src/RecoPixelVertexing/PixelTriplets/src/CosmicHitTripletGeneratorFromLayerTriplet.cc

Go to the documentation of this file.

#include "RecoPixelVertexing/PixelTriplets/interface/CosmicHitTripletGeneratorFromLayerTriplet.h"
#include "RecoTracker/TkTrackingRegions/interface/TrackingRegion.h"
#include "TrackingTools/DetLayers/interface/DetLayer.h"
#include "RecoPixelVertexing/PixelTriplets/interface/OrderedHitTriplets.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
#include "TrackingTools/Records/interface/TransientRecHitRecord.h"

#include <cmath>

typedef TransientTrackingRecHit::ConstRecHitPointer TkHitPairsCachedHit;

CosmicHitTripletGeneratorFromLayerTriplet::CosmicHitTripletGeneratorFromLayerTriplet
(const LayerWithHits* inner, 
 const LayerWithHits* middle, 
 const LayerWithHits* outer,
 const edm::EventSetup& iSetup)
  : TTRHbuilder(0),trackerGeometry(0),
    //theLayerCache(*layerCache), 
    theOuterLayer(outer),theMiddleLayer(middle), theInnerLayer(inner)
{

  edm::ESHandle<TrackerGeometry> tracker;
  iSetup.get<TrackerDigiGeometryRecord>().get(tracker);
  trackerGeometry = tracker.product();
}
void CosmicHitTripletGeneratorFromLayerTriplet::hitTriplets(
  const TrackingRegion & region, OrderedHitTriplets & result,
  const edm::EventSetup& iSetup)
{
  
  if (theInnerLayer->recHits().empty()) return;
  if (theMiddleLayer->recHits().empty()) return;
  if (theOuterLayer->recHits().empty()) return;
  float radius1 =dynamic_cast<const BarrelDetLayer*>(theInnerLayer->layer())->specificSurface().radius();
  float radius2 =dynamic_cast<const BarrelDetLayer*>(theMiddleLayer->layer())->specificSurface().radius();
  float radius3 =dynamic_cast<const BarrelDetLayer*>(theOuterLayer->layer())->specificSurface().radius();
  bool seedfromoverlaps=((std::abs(radius1-radius2)<0.1)|| (std::abs(radius3-radius2)<0.1))? true : false;
  std::vector<const TrackingRecHit*>::const_iterator ohh;
  std::vector<const TrackingRecHit*>::const_iterator mhh;
  std::vector<const TrackingRecHit*>::const_iterator ihh;

  std::string builderName = "WithTrackAngle";
  edm::ESHandle<TransientTrackingRecHitBuilder> builder;
  iSetup.get<TransientRecHitRecord>().get(builderName, builder);

  if(!seedfromoverlaps){
    for(ohh=theOuterLayer->recHits().begin();ohh!=theOuterLayer->recHits().end();ohh++){
      TkHitPairsCachedHit oh= builder->build(*ohh);
      for(mhh=theMiddleLayer->recHits().begin();mhh!=theMiddleLayer->recHits().end();mhh++){
        TkHitPairsCachedHit mh=  builder->build(*mhh);
        float z_diff =mh->globalPosition().z()-oh->globalPosition().z();
        float midy=mh->globalPosition().y();
        float outy=oh->globalPosition().y();
        float midx=mh->globalPosition().x();
        float outx=oh->globalPosition().x();
        float dxdy=std::abs((outx-midx)/(outy-midy));
        if((std::abs(z_diff)<30) && (outy*midy>0) &&(dxdy<2))     
          {
            for(ihh=theInnerLayer->recHits().begin();ihh!=theInnerLayer->recHits().end();ihh++){
              TkHitPairsCachedHit ih= builder->build(*ihh);
              float z_diff =mh->globalPosition().z()-ih->globalPosition().z();
              float inny=ih->globalPosition().y();
              float innx=ih->globalPosition().x();
              float dxdy=std::abs((innx-midx)/(inny-midy));
              if ((std::abs(z_diff)<30) && (inny*midy>0) &&(dxdy<2)&&(!seedfromoverlaps))
                {
                  result.push_back( OrderedHitTriplet(ih,mh,oh) );
                }
            } 
          }
      }
    }
  } else {
    for(ohh=theOuterLayer->recHits().begin();ohh!=theOuterLayer->recHits().end();ohh++){
      TkHitPairsCachedHit oh= builder->build(*ohh);
      for(mhh=theMiddleLayer->recHits().begin();mhh!=theMiddleLayer->recHits().end();mhh++){
        TkHitPairsCachedHit mh= builder->build(*mhh);
        float z_diff =mh->globalPosition().z()-oh->globalPosition().z();
        float midy=mh->globalPosition().y();
        float outy=oh->globalPosition().y();
        float midx=mh->globalPosition().x();
        float outx=oh->globalPosition().x();
        float dxdy=std::abs((outx-midx)/(outy-midy));
        float DeltaR=oh->globalPosition().perp()-mh->globalPosition().perp();
        if((std::abs(z_diff)<18) && (std::abs(oh->globalPosition().phi()-mh->globalPosition().phi())<0.05) &&(DeltaR<0)&&(dxdy<2)){
          for(ihh=theInnerLayer->recHits().begin();ihh!=theInnerLayer->recHits().end();ihh++){
            TkHitPairsCachedHit ih= builder->build(*ihh);
            float z_diff =mh->globalPosition().z()-ih->globalPosition().z();
            float inny=ih->globalPosition().y();
            float innx=ih->globalPosition().x();
            float dxdy=std::abs((innx-midx)/(inny-midy));
            if ((std::abs(z_diff)<30) && (inny*midy>0) &&(dxdy<2))
              {
                result.push_back( OrderedHitTriplet(ih,mh,oh));
              }
          }
        }
      }
    }
  }
}