#include <PixelTripletLargeTipGenerator.h>

Inheritance diagram for PixelTripletLargeTipGenerator:

Public Member Functions
virtual void	hitTriplets (const TrackingRegion &region, OrderedHitTriplets &trs, const edm::Event &ev, const edm::EventSetup &es, SeedingLayerSetsHits::SeedingLayerSet pairLayers, const std::vector< SeedingLayerSetsHits::SeedingLayer > &thirdLayers) override

	PixelTripletLargeTipGenerator (const edm::ParameterSet &cfg, edm::ConsumesCollector &iC)

virtual	~PixelTripletLargeTipGenerator ()

Public Member Functions inherited from HitTripletGeneratorFromPairAndLayers
	HitTripletGeneratorFromPairAndLayers (unsigned int maxElement=0)

	HitTripletGeneratorFromPairAndLayers (const edm::ParameterSet &pset)

void	init (std::unique_ptr< HitPairGeneratorFromLayerPair > &&pairs, LayerCacheType *layerCache)

const HitPairGeneratorFromLayerPair &	pairGenerator () const

virtual	~HitTripletGeneratorFromPairAndLayers ()

Private Types
typedef CombinedHitTripletGenerator::LayerCacheType	LayerCacheType

Private Attributes
const float	dphi

const float	extraHitRPhitolerance

const float	extraHitRZtolerance

const bool	useBend

const bool	useFixedPreFiltering

const bool	useMScat

Additional Inherited Members
Public Types inherited from HitTripletGeneratorFromPairAndLayers
typedef LayerHitMapCache	LayerCacheType

Protected Attributes inherited from HitTripletGeneratorFromPairAndLayers
LayerCacheType *	theLayerCache

const unsigned int	theMaxElement

std::unique_ptr < HitPairGeneratorFromLayerPair >	thePairGenerator

Detailed Description

A HitTripletGenerator from HitPairGenerator and vector of Layers. The HitPairGenerator provides a set of hit pairs. For each pair the search for compatible hit(s) is done among provided Layers

Definition at line 19 of file PixelTripletLargeTipGenerator.h.

Member Typedef Documentation

typedef CombinedHitTripletGenerator::LayerCacheType PixelTripletLargeTipGenerator::LayerCacheType

private

Definition at line 21 of file PixelTripletLargeTipGenerator.h.

Constructor & Destructor Documentation

PixelTripletLargeTipGenerator::PixelTripletLargeTipGenerator	(	const edm::ParameterSet &	cfg,
		edm::ConsumesCollector &	iC
	)

Definition at line 48 of file PixelTripletLargeTipGenerator.cc.

   : HitTripletGeneratorFromPairAndLayers(cfg),
     useFixedPreFiltering(cfg.getParameter<bool>("useFixedPreFiltering")),
     extraHitRZtolerance(cfg.getParameter<double>("extraHitRZtolerance")),
     extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")),
     useMScat(cfg.getParameter<bool>("useMultScattering")),
     useBend(cfg.getParameter<bool>("useBending")),
     dphi(useFixedPreFiltering ? cfg.getParameter<double>("phiPreFiltering") : 0) {}

PixelTripletLargeTipGenerator::~PixelTripletLargeTipGenerator ( )

virtual

Definition at line 57 of file PixelTripletLargeTipGenerator.cc.

57 {}

Member Function Documentation

void PixelTripletLargeTipGenerator::hitTriplets	(	const TrackingRegion &	region,
		OrderedHitTriplets &	trs,
		const edm::Event &	ev,
		const edm::EventSetup &	es,
		SeedingLayerSetsHits::SeedingLayerSet	pairLayers,
		const std::vector< SeedingLayerSetsHits::SeedingLayer > &	thirdLayers
	)

overridevirtual

Implements HitTripletGeneratorFromPairAndLayers.

Definition at line 73 of file PixelTripletLargeTipGenerator.cc.

References angle(), checkPhiInRange(), PixelRecoUtilities::curvature(), declareDynArray, dphi, relativeConstraints::empty, ev, extraHitRPhitolerance, extraHitRZtolerance, f, plotBeamSpotDB::first, fnSigmaRZ, Geom::ftwoPi(), edm::EventSetup::get(), i, HitDoublets::inner, DetLayer::isBarrel(), geometryCSVtoXML::line, M_PI, bookConverter::max, ThirdHitPredictionFromCircle::HelixRZ::maxCurvature(), min(), normalizedPhi(), nSigmaPhi, nSigmaRZ, PixelRecoLineRZ::origin(), TrackingRegion::originRBound(), HitDoublets::outer, p3, PV3DBase< T, PVType, FrameType >::perp(), GeometricSearchDet::position(), edm::ESHandle< class >::product(), proxim(), TrackingRegion::ptMin(), Basic2DVector< T >::r(), CosmicsPD_Skims::radius, HLT_FULL_cff::region, OrderedHitTriplets::size(), findQualityFiles::size, std::swap(), HitTripletGeneratorFromPairAndLayers::theMaxElement, HitTripletGeneratorFromPairAndLayers::thePairGenerator, patCandidatesForDimuonsSequences_cff::tracker, useBend, useFixedPreFiltering, useMScat, findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

 { 
   edm::ESHandle<TrackerGeometry> tracker;
   es.get<TrackerDigiGeometryRecord>().get(tracker);
 
   //Retrieve tracker topology from geometry
   edm::ESHandle<TrackerTopology> tTopoHand;
   es.get<TrackerTopologyRcd>().get(tTopoHand);
   const TrackerTopology *tTopo=tTopoHand.product();
 
   auto const & doublets = thePairGenerator->doublets(region,ev,es, pairLayers);
   
   if (doublets.empty()) return;
    
   auto outSeq =  doublets.detLayer(HitDoublets::outer)->seqNum();
 
 
   int size = thirdLayers.size();
 
 
   using NodeInfo = KDTreeNodeInfo<unsigned int>;
   std::vector<NodeInfo > layerTree; // re-used throughout
   std::vector<unsigned int> foundNodes; // re-used throughout
   foundNodes.reserve(100);
 
   declareDynArray(KDTreeLinkerAlgo<unsigned int>, size, hitTree);
   declareDynArray(LayerRZPredictions, size, mapPred);
 
   float rzError[size]; //save maximum errors
 
   const float maxDelphi = region.ptMin() < 0.3f ? float(M_PI)/4.f : float(M_PI)/8.f; // FIXME move to config??
   const float maxphi = M_PI+maxDelphi, minphi = -maxphi; // increase to cater for any range
   const float safePhi = M_PI-maxDelphi; // sideband
 
 
   const RecHitsSortedInPhi * thirdHitMap[size];
 
   for(int il = 0; il < size; il++) {
     thirdHitMap[il] = &(*theLayerCache)(thirdLayers[il], region, ev, es);
     auto const & hits = *thirdHitMap[il];
  
     const DetLayer *layer = thirdLayers[il].detLayer();
     LayerRZPredictions &predRZ = mapPred[il];
     predRZ.line.initLayer(layer);
     predRZ.helix1.initLayer(layer);
     predRZ.helix2.initLayer(layer);
     predRZ.line.initTolerance(extraHitRZtolerance);
     predRZ.helix1.initTolerance(extraHitRZtolerance);
     predRZ.helix2.initTolerance(extraHitRZtolerance);
     predRZ.rzPositionFixup = MatchedHitRZCorrectionFromBending(layer,tTopo);
     predRZ.correction.init(es, region.ptMin(), *doublets.detLayer(HitDoublets::inner), *doublets.detLayer(HitDoublets::outer), *thirdLayers[il].detLayer(), useMScat, false);
 
 
     layerTree.clear();
     float minv=999999.0; float maxv = -999999.0; // Initialise to extreme values in case no hits
     float maxErr=0.0f;
     for (unsigned int i=0; i!=hits.size(); ++i) {
       auto angle = hits.phi(i);
       auto v =  hits.gv(i);
       //use (phi,r) for endcaps rather than (phi,z)
       minv = std::min(minv,v);  maxv = std::max(maxv,v);
       float myerr = hits.dv[i];
       maxErr = std::max(maxErr,myerr);
       layerTree.emplace_back(i, angle, v); // save it
       // populate side-bands
       if (angle>safePhi) layerTree.emplace_back(i, angle-Geom::ftwoPi(), v);
       else if (angle<-safePhi) layerTree.emplace_back(i, angle+Geom::ftwoPi(), v);
     }
     KDTreeBox phiZ(minphi, maxphi, minv-0.01f, maxv+0.01f);  // declare our bounds
     //add fudge factors in case only one hit and also for floating-point inaccuracy
     hitTree[il].build(layerTree, phiZ); // make KDtree
     rzError[il] = maxErr; //save error
   }
 
   double curv = PixelRecoUtilities::curvature(1. / region.ptMin(), es);
   
   for (std::size_t ip =0;  ip!=doublets.size(); ip++) {
     auto xi = doublets.x(ip,HitDoublets::inner);
     auto yi = doublets.y(ip,HitDoublets::inner);
     auto zi = doublets.z(ip,HitDoublets::inner);
     // auto rvi = doublets.rv(ip,HitDoublets::inner);
     auto xo = doublets.x(ip,HitDoublets::outer);
     auto yo = doublets.y(ip,HitDoublets::outer);
     auto zo = doublets.z(ip,HitDoublets::outer);
     // auto rvo = doublets.rv(ip,HitDoublets::outer);
     GlobalPoint gp1(xi,yi,zi);
     GlobalPoint gp2(xo,yo,zo);
 
     auto toPos = std::signbit(zo-zi); 
 
     PixelRecoLineRZ line(gp1, gp2);
     PixelRecoPointRZ point2(gp2.perp(), zo);
     ThirdHitPredictionFromCircle predictionRPhi(gp1, gp2, extraHitRPhitolerance);
 
     Range generalCurvature = predictionRPhi.curvature(region.originRBound());
     if (!intersect(generalCurvature, Range(-curv, curv))) continue;
 
     for(int il = 0; il < size; il++) {
       if (hitTree[il].empty()) continue; // Don't bother if no hits
       const DetLayer *layer = thirdLayers[il].detLayer();
       bool barrelLayer = layer->isBarrel();
 
       if ( (!barrelLayer) & (toPos != std::signbit(layer->position().z())) ) continue;
 
       
       Range curvature = generalCurvature;
 
       LayerRZPredictions &predRZ = mapPred[il];
       predRZ.line.initPropagator(&line);
 
       auto & correction = predRZ.correction;
       correction.init(line, point2,  outSeq);
 
 
       Range rzRange;
       if (useBend) {
         // For the barrel region:
         // swiping the helix passing through the two points across from
         // negative to positive bending, can give us a sort of U-shaped
         // projection onto the phi-z (barrel) or r-z plane (forward)
         // so we checking minimum/maximum of all three possible extrema
         // 
         // For the endcap region:
         // Checking minimum/maximum radius of the helix projection
         // onto an endcap plane, here we have to guard against
         // looping tracks, when phi(delta z) gets out of control.
         // HelixRZ::rAtZ should not follow looping tracks, but clamp
         // to the minimum reachable r with the next-best lower |curvature|.
         // So same procedure as for the barrel region can be applied.
         //
         // In order to avoid looking for potential looping tracks at all
         // we also clamp the allowed curvature range for this layer,
         // and potentially fail the layer entirely
     
         if (!barrelLayer) {
           Range z3s = predRZ.line.detRange();
           double z3 = z3s.first < 0 ? std::max(z3s.first, z3s.second)
         : std::min(z3s.first, z3s.second);
           double maxCurvature = HelixRZ::maxCurvature(&predictionRPhi,
                                                       gp1.z(), gp2.z(), z3);
           if (!intersect(curvature, Range(-maxCurvature, maxCurvature)))
             continue;
         }
     
         HelixRZ helix1(&predictionRPhi, gp1.z(), gp2.z(), curvature.first);
         HelixRZ helix2(&predictionRPhi, gp1.z(), gp2.z(), curvature.second);
     
         predRZ.helix1.initPropagator(&helix1);
         predRZ.helix2.initPropagator(&helix2);
     
         Range rzRanges[2] = { predRZ.helix1(), predRZ.helix2() };
         predRZ.helix1.initPropagator(nullptr);
         predRZ.helix2.initPropagator(nullptr);
 
         rzRange.first = std::min(rzRanges[0].first, rzRanges[1].first);
         rzRange.second = std::max(rzRanges[0].second, rzRanges[1].second);
     
         // if the allowed curvatures include a straight line,
         // this can give us another extremum for allowed r/z
         if (curvature.first * curvature.second < 0.0) {
           Range rzLineRange = predRZ.line();
           rzRange.first = std::min(rzRange.first, rzLineRange.first);
           rzRange.second = std::max(rzRange.second, rzLineRange.second);
         }
       } else {
         rzRange = predRZ.line();
       }
 
       if (rzRange.first >= rzRange.second)
         continue;
 
       correction.correctRZRange(rzRange);
 
       Range phiRange;
       if (useFixedPreFiltering) {
     float phi0 = doublets.phi(ip,HitDoublets::outer);
         phiRange = Range(phi0 - dphi, phi0 + dphi);
       } else {
         Range radius;
     
         if (barrelLayer) {
           radius = predRZ.line.detRange();
           if (!intersect(rzRange, predRZ.line.detSize()))
             continue;
         } else {
           radius = rzRange;
           if (!intersect(radius, predRZ.line.detSize()))
             continue;
         }
          
         //gc: predictionRPhi uses the cosine rule to find the phi of the 3rd point at radius, assuming the pairCurvature range [-c,+c]
         if ( (curvature.first<0.0f) & (curvature.second<0.0f) ) {
           radius.swap();
         } else if ( (curvature.first>=0.0f) & (curvature.second>=0.0f) ) {;}
         else {
           radius.first=radius.second;
         }
         auto phi12 = predictionRPhi.phi(curvature.first,radius.first);
         auto phi22 = predictionRPhi.phi(curvature.second,radius.second);
         phi12 = normalizedPhi(phi12);
         phi22 = proxim(phi22,phi12);
         phiRange = Range(phi12,phi22); phiRange.sort();
         auto rmean = radius.mean();
         phiRange.first *= rmean;
         phiRange.second *= rmean;
     correction.correctRPhiRange(phiRange);
         phiRange.first /= rmean;
         phiRange.second /= rmean;
 
       }
       
       foundNodes.clear(); // Now recover hits in bounding box...
       float prmin=phiRange.min(), prmax=phiRange.max(); //get contiguous range
 
       if (prmax-prmin>maxDelphi) {
         auto prm = phiRange.mean();
         prmin = prm - 0.5f*maxDelphi;
         prmax = prm + 0.5f*maxDelphi;
       }
 
       if (barrelLayer) {
     Range regMax = predRZ.line.detRange();
     Range regMin = predRZ.line(regMax.min());
     regMax = predRZ.line(regMax.max());
     correction.correctRZRange(regMin);
     correction.correctRZRange(regMax);
     if (regMax.min() < regMin.min()) { std::swap(regMax, regMin);}
     KDTreeBox phiZ(prmin, prmax,
                regMin.min()-fnSigmaRZ*rzError[il],
                regMax.max()+fnSigmaRZ*rzError[il]);
     hitTree[il].search(phiZ, foundNodes);
       }
       else {
     KDTreeBox phiZ(prmin, prmax,
                rzRange.min()-fnSigmaRZ*rzError[il],
                rzRange.max()+fnSigmaRZ*rzError[il]);
     hitTree[il].search(phiZ, foundNodes);
       }
       
       MatchedHitRZCorrectionFromBending l2rzFixup(doublets.hit(ip,HitDoublets::outer)->det()->geographicalId(), tTopo);
       MatchedHitRZCorrectionFromBending l3rzFixup = predRZ.rzPositionFixup;
 
       thirdHitMap[il] = &(*theLayerCache)(thirdLayers[il], region, ev, es);
       auto const & hits = *thirdHitMap[il];
       for (auto KDdata : foundNodes) {
     GlobalPoint p3 = hits.gp(KDdata); 
     double p3_r = p3.perp();
     double p3_z = p3.z();
     float p3_phi =  hits.phi(KDdata); 
 
     Range rangeRPhi = predictionRPhi(curvature, p3_r);
     correction.correctRPhiRange(rangeRPhi);
 
     float ir = 1.f/p3_r;
     float phiErr = nSigmaPhi *  hits.drphi[KDdata]*ir;
     if (!checkPhiInRange(p3_phi, rangeRPhi.first*ir-phiErr, rangeRPhi.second*ir+phiErr))
       continue;
     
     Basic2DVector<double> thc(p3.x(), p3.y());
     
     auto curv_ = predictionRPhi.curvature(thc);
     double p2_r = point2.r(); double p2_z = point2.z(); // they will be modified!
     
     l2rzFixup(predictionRPhi, curv_, *doublets.hit(ip,HitDoublets::outer), p2_r, p2_z, tTopo);
     l3rzFixup(predictionRPhi, curv_, *hits.theHits[KDdata].hit(), p3_r, p3_z, tTopo);
     
     Range rangeRZ;
     if (useBend) {
       HelixRZ updatedHelix(&predictionRPhi, gp1.z(), p2_z, curv_);
       rangeRZ = predRZ.helix1(barrelLayer ? p3_r : p3_z, updatedHelix);
     } else {
       float tIP = predictionRPhi.transverseIP(thc);
       PixelRecoPointRZ updatedPoint2(p2_r, p2_z);
       PixelRecoLineRZ updatedLine(line.origin(), point2, tIP);
       rangeRZ = predRZ.line(barrelLayer ? p3_r : p3_z, line);
     }
     correction.correctRZRange(rangeRZ);
     
     double err = nSigmaRZ * hits.dv[KDdata];
     
     rangeRZ.first -= err, rangeRZ.second += err;
     
     if (!rangeRZ.inside(barrelLayer ? p3_z : p3_r)) continue;
 
     if (theMaxElement!=0 && result.size() >= theMaxElement) {
       result.clear();
       edm::LogError("TooManyTriplets")<<" number of triples exceed maximum. no triplets produced.";
       return;
     }
     result.emplace_back( doublets.hit(ip,HitDoublets::inner), doublets.hit(ip,HitDoublets::outer), hits.theHits[KDdata].hit()); 
       }
     }
   }
   // std::cout << "found triplets " << result.size() << std::endl;
 }