#include <MultiHitGeneratorFromChi2.h>

Inheritance diagram for MultiHitGeneratorFromChi2:

Public Member Functions
virtual void	hitSets (const TrackingRegion &region, OrderedMultiHits &trs, const edm::Event &ev, const edm::EventSetup &es)

void	init (const HitPairGenerator &pairs, LayerCacheType *layerCache) override

void	initES (const edm::EventSetup &es) override

	MultiHitGeneratorFromChi2 (const edm::ParameterSet &cfg)

const HitPairGenerator &	pairGenerator () const

void	setSeedingLayers (SeedingLayerSetsHits::SeedingLayerSet pairLayers, std::vector< SeedingLayerSetsHits::SeedingLayer > thirdLayers) override

virtual	~MultiHitGeneratorFromChi2 ()

Public Member Functions inherited from MultiHitGeneratorFromPairAndLayers
virtual	~MultiHitGeneratorFromPairAndLayers ()

Public Member Functions inherited from MultiHitGenerator
virtual void	clear ()

virtual void	hitSets (const TrackingRegion &reg, OrderedMultiHits &prs, const edm::EventSetup &es)

	MultiHitGenerator (unsigned int size=400)

	MultiHitGenerator (MultiHitGenerator const &other)

virtual const OrderedMultiHits &	run (const TrackingRegion &region, const edm::Event &ev, const edm::EventSetup &es) final

virtual	~MultiHitGenerator ()

Public Member Functions inherited from OrderedHitsGenerator
	OrderedHitsGenerator ()

virtual	~OrderedHitsGenerator ()

Private Types
using	HitOwnPtr = mayown_ptr< BaseTrackerRecHit >

typedef CombinedMultiHitGenerator::LayerCacheType	LayerCacheType

Private Member Functions
bool	checkPhiInRange (float phi, float phi1, float phi2) const

std::pair< float, float >	mergePhiRanges (const std::pair< float, float > &r1, const std::pair< float, float > &r2) const

void	refit2Hits (HitOwnPtr &hit0, HitOwnPtr &hit1, TrajectoryStateOnSurface &tsos0, TrajectoryStateOnSurface &tsos1, const TrackingRegion &region, float nomField, bool isDebug)

Private Attributes
const MagneticField *	bfield

TkTransientTrackingRecHitBuilder const *	builder

std::string	builderName_

std::vector< double >	chi2_cuts

bool	chi2VsPtCut

TkClonerImpl	cloner

bool	debug

std::vector< int >	detIdsToDebug

float	dphi

float	extraHitRPhitolerance

float	extraHitRZtolerance

float	extraPhiKDBox

float	extraRKDBox

float	extraZKDBox

const ClusterShapeHitFilter *	filter

std::string	filterName_

double	fnSigmaRZ

double	maxChi2

std::string	mfName_

float	nomField

double	nSigmaPhi

double	nSigmaRZ

std::vector< double >	pt_interv

bool	refitHits

LayerCacheType *	theLayerCache

std::vector < SeedingLayerSetsHits::SeedingLayer >	theLayers

HitPairGenerator *	thePairGenerator

bool	useFixedPreFiltering

bool	useSimpleMF_

Additional Inherited Members
Public Types inherited from MultiHitGeneratorFromPairAndLayers
typedef LayerHitMapCache	LayerCacheType

Public Attributes inherited from OrderedHitsGenerator
unsigned int	theMaxElement

Protected Types inherited from MultiHitGenerator
using	cacheHitPointer = std::unique_ptr< BaseTrackerRecHit >

using	cacheHits = std::vector< cacheHitPointer >

Protected Attributes inherited from MultiHitGenerator
cacheHits	cache

edm::RunningAverage	localRA

Detailed Description

A MultiHitGenerator 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 27 of file MultiHitGeneratorFromChi2.h.

Member Typedef Documentation

using MultiHitGeneratorFromChi2::HitOwnPtr = mayown_ptr<BaseTrackerRecHit>

private

Definition at line 49 of file MultiHitGeneratorFromChi2.h.

typedef CombinedMultiHitGenerator::LayerCacheType MultiHitGeneratorFromChi2::LayerCacheType

private

Definition at line 29 of file MultiHitGeneratorFromChi2.h.

Constructor & Destructor Documentation

MultiHitGeneratorFromChi2::MultiHitGeneratorFromChi2 ( const edm::ParameterSet & cfg )

Definition at line 50 of file MultiHitGeneratorFromChi2.cc.

References bfield, chi2_cuts, chi2VsPtCut, debug, detIdsToDebug, dphi, edm::ParameterSet::exists(), filter, edm::ParameterSet::getParameter(), mfName_, nomField, pt_interv, AlCaHLTBitMon_QueryRunRegistry::string, OrderedHitsGenerator::theMaxElement, useFixedPreFiltering, and useSimpleMF_.

   : thePairGenerator(0),
     theLayerCache(0),
     useFixedPreFiltering(cfg.getParameter<bool>("useFixedPreFiltering")),
     extraHitRZtolerance(cfg.getParameter<double>("extraHitRZtolerance")),//extra window in ThirdHitRZPrediction range 
     extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")),//extra window in ThirdHitPredictionFromCircle range (divide by R to get phi) 
     extraZKDBox(cfg.getParameter<double>("extraZKDBox")),//extra windown in Z when building the KDTree box (used in barrel)
     extraRKDBox(cfg.getParameter<double>("extraRKDBox")),//extra windown in R when building the KDTree box (used in endcap)
     extraPhiKDBox(cfg.getParameter<double>("extraPhiKDBox")),//extra windown in Phi when building the KDTree box
     fnSigmaRZ(cfg.getParameter<double>("fnSigmaRZ")),//this multiplies the max hit error on the layer when building the KDTree box
     chi2VsPtCut(cfg.getParameter<bool>("chi2VsPtCut")),
     maxChi2(cfg.getParameter<double>("maxChi2")),
     refitHits(cfg.getParameter<bool>("refitHits")),
     debug(cfg.getParameter<bool>("debug")),
     filterName_(cfg.getParameter<std::string>("ClusterShapeHitFilterName")),
     builderName_(cfg.existsAs<std::string>("TTRHBuilder") ? cfg.getParameter<std::string>("TTRHBuilder") : std::string("WithTrackAngle")),
     useSimpleMF_(false),
     mfName_("")
 {    
   theMaxElement=cfg.getParameter<unsigned int>("maxElement");
   if (useFixedPreFiltering)
     dphi = cfg.getParameter<double>("phiPreFiltering");
   if (chi2VsPtCut) {
     pt_interv = cfg.getParameter<std::vector<double> >("pt_interv");
     chi2_cuts = cfg.getParameter<std::vector<double> >("chi2_cuts");    
   }  
   if (debug) {
     detIdsToDebug = cfg.getParameter<std::vector<int> >("detIdsToDebug");
     //if (detIdsToDebug.size()<3) //fixme
   } else {
     detIdsToDebug.push_back(0);
     detIdsToDebug.push_back(0);
     detIdsToDebug.push_back(0);
   }
   // 2014/02/11 mia:
   // we should get rid of the boolean parameter useSimpleMF,
   // and use only a string magneticField [instead of SimpleMagneticField]
   // or better an edm::ESInputTag (at the moment HLT does not handle ESInputTag)
   if (cfg.exists("SimpleMagneticField")) {
     useSimpleMF_ = true;
     mfName_ = cfg.getParameter<std::string>("SimpleMagneticField");
   }
   filter = 0;
   bfield = 0;
   nomField = -1.;
 }

virtual MultiHitGeneratorFromChi2::~MultiHitGeneratorFromChi2 ( )

inlinevirtual

Definition at line 34 of file MultiHitGeneratorFromChi2.h.

34 { delete thePairGenerator; }

MultiHitGeneratorFromChi2::thePairGenerator

HitPairGenerator * thePairGenerator

Definition: MultiHitGeneratorFromChi2.h:70

Member Function Documentation

bool MultiHitGeneratorFromChi2::checkPhiInRange	(	float	phi,
		float	phi1,
		float	phi2
	)		const

private

Definition at line 648 of file MultiHitGeneratorFromChi2.cc.

References M_PI.

 { while (phi > phi2) phi -= 2. * M_PI;
   while (phi < phi1) phi += 2. * M_PI;
   return phi <= phi2;
 }  

void MultiHitGeneratorFromChi2::hitSets	(	const TrackingRegion &	region,
		OrderedMultiHits &	trs,
		const edm::Event &	ev,
		const edm::EventSetup &	es
	)

virtual

Implements MultiHitGenerator.

Definition at line 143 of file MultiHitGeneratorFromChi2.cc.

References RecHitsSortedInPhi::all(), angle(), assert(), GeomDetEnumerators::barrel, bfield, KDTreeLinkerAlgo< DATA >::build(), MultiHitGenerator::cache, RZLine::chi2(), chi2_cuts, chi2VsPtCut, cloner, gather_cfg::cout, ThirdHitPredictionFromCircle::curvature(), PixelRecoUtilities::curvature(), debug, detIdsToDebug, dphi, relativeConstraints::empty, mayown_ptr< T, N >::empty(), ev, extraHitRPhitolerance, extraHitRZtolerance, extraPhiKDBox, extraRKDBox, extraZKDBox, filter, RZLine::fit(), fnSigmaRZ, ClusterShapeHitFilter::isCompatible(), edm::isNotFinite(), mayown_ptr< T, N >::isOwn(), geometryCSVtoXML::line, TrajectoryStateOnSurface::localMomentum(), DetLayer::location(), M_PI, bookConverter::max, maxChi2, SiStripMatchedRecHit2D::monoCluster(), SiStripMatchedRecHit2D::monoId(), mergeVDriftHistosByStation::name, nomField, ProjectedSiStripRecHit2D::originalHit(), TrackingRegion::originRBound(), PV3DBase< T, PVType, FrameType >::perp(), ThirdHitPredictionFromCircle::phi(), Geom::pi(), PixelSubdetector::PixelBarrel, EnergyCorrector::pt, pt_interv, TrackingRegion::ptMin(), CosmicsPD_Skims::radius, refit2Hits(), refitHits, HLT_50ns_5e33_v3_cff::region, mayown_ptr< T, N >::release(), mayown_ptr< T, N >::reset(), ProxyBase11< T >::reset(), rho, KDTreeLinkerAlgo< DATA >::search(), OrderedHitPairs::size(), OrderedMultiHits::size(), findQualityFiles::size, SiStripMatchedRecHit2D::stereoCluster(), SiStripMatchedRecHit2D::stereoId(), theLayers, OrderedHitsGenerator::theMaxElement, thePairGenerator, StripSubdetector::TIB, SiStripDetId::TIB, SiStripDetId::TID, StripSubdetector::TOB, Geom::twoPi(), unlikely, useFixedPreFiltering, and PV3DBase< T, PVType, FrameType >::z().

 { 
 
   unsigned int debug_Id0 = detIdsToDebug[0];
   unsigned int debug_Id1 = detIdsToDebug[1];
   unsigned int debug_Id2 = detIdsToDebug[2];
 
   if (debug) cout << "pair: " << ((HitPairGeneratorFromLayerPair*) thePairGenerator)->innerLayer().name() << "+" <<  ((HitPairGeneratorFromLayerPair*) thePairGenerator)->outerLayer().name() << " 3rd lay size: " << theLayers.size() << endl;
 
   //gc: first get the pairs
   OrderedHitPairs pairs;
   pairs.reserve(30000);
   thePairGenerator->hitPairs(region,pairs,ev,es);
   if (debug) cout << endl;
   if (pairs.empty()) {
     //cout << "empy pairs" << endl;
     return;
   }
   
   //gc: these are all the layers compatible with the layer pairs (as defined in the config file)
   int size = theLayers.size();
 
 
   //gc: initialize a KDTree per each 3rd layer
   std::vector<KDTreeNodeInfo<RecHitsSortedInPhi::HitIter> > layerTree; // re-used throughout
   std::vector<RecHitsSortedInPhi::HitIter> foundNodes; // re-used thoughout
   foundNodes.reserve(100);
   #ifdef __clang__
   std::vector<KDTreeLinkerAlgo<RecHitsSortedInPhi::HitIter>> hitTree(size);
   #else
   KDTreeLinkerAlgo<RecHitsSortedInPhi::HitIter> hitTree[size];
   #endif
   float rzError[size]; //save maximum errors
   double maxphi = Geom::twoPi(), minphi = -maxphi; //increase to cater for any range
 
   map<std::string, LayerRZPredictions> mapPred;//need to use the name as map key since we may have more than one SeedingLayer per DetLayer (e.g. TID and MTID)
   const RecHitsSortedInPhi * thirdHitMap[size];//gc: this comes from theLayerCache
 
   //gc: loop over each layer
   for(int il = 0; il < size; il++) {
     thirdHitMap[il] = &(*theLayerCache)(theLayers[il], region, ev, es);
     if (debug) cout << "considering third layer: " << theLayers[il].name() << " with hits: " << thirdHitMap[il]->all().second-thirdHitMap[il]->all().first << endl;
     const DetLayer *layer = theLayers[il].detLayer();
     LayerRZPredictions &predRZ = mapPred[theLayers[il].name()];
     predRZ.line.initLayer(layer);
     predRZ.line.initTolerance(extraHitRZtolerance);
 
     //gc: now we take all hits in the layer and fill the KDTree    
     RecHitsSortedInPhi::Range hitRange = thirdHitMap[il]->all(); // Get iterators
     layerTree.clear();
     double minz=999999.0, maxz= -999999.0; // Initialise to extreme values in case no hits
     float maxErr=0.0f;
     bool barrelLayer = (theLayers[il].detLayer()->location() == GeomDetEnumerators::barrel);
     if (hitRange.first != hitRange.second)
       { minz = barrelLayer? hitRange.first->hit()->globalPosition().z() : hitRange.first->hit()->globalPosition().perp();
     maxz = minz; //In case there's only one hit on the layer
     for (RecHitsSortedInPhi::HitIter hi=hitRange.first; hi != hitRange.second; ++hi)
       { double angle = hi->phi();
         double myz = barrelLayer? hi->hit()->globalPosition().z() : hi->hit()->globalPosition().perp();
 
         if (debug && hi->hit()->rawId()==debug_Id2) {
           cout << "filling KDTree with hit in id=" << debug_Id2 
            << " with pos: " << hi->hit()->globalPosition() 
            << " phi=" << hi->hit()->globalPosition().phi() 
            << " z=" << hi->hit()->globalPosition().z() 
            << " r=" << hi->hit()->globalPosition().perp() 
            << endl;
         }
         //use (phi,r) for endcaps rather than (phi,z)
         if (myz < minz) { minz = myz;} else { if (myz > maxz) {maxz = myz;}}
         float myerr = barrelLayer? hi->hit()->errorGlobalZ(): hi->hit()->errorGlobalR();
         if (myerr > maxErr) { maxErr = myerr;}
         layerTree.push_back(KDTreeNodeInfo<RecHitsSortedInPhi::HitIter>(hi, angle, myz)); // save it
         if (angle < 0)  // wrap all points in phi
           { layerTree.push_back(KDTreeNodeInfo<RecHitsSortedInPhi::HitIter>(hi, angle+Geom::twoPi(), myz));}
         else
           { layerTree.push_back(KDTreeNodeInfo<RecHitsSortedInPhi::HitIter>(hi, angle-Geom::twoPi(), myz));}
       }
       }
     KDTreeBox phiZ(minphi, maxphi, minz-0.01, maxz+0.01);  // 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
   }
   //gc: now we have initialized the KDTrees and we are out of the layer loop
   
   //gc: this sets the minPt of the triplet
   double curv = PixelRecoUtilities::curvature(1. / region.ptMin(), es);
 
   if (debug) std::cout << "pair size=" << pairs.size() << std::endl;
 
   //gc: now we loop over all pairs
   for (OrderedHitPairs::const_iterator ip = pairs.begin(); ip != pairs.end(); ++ip) {
 
     int foundTripletsFromPair = 0;
     bool usePair = false;
     cacheHitPointer bestH2;
     float minChi2 = std::numeric_limits<float>::max();
 
     SeedingHitSet::ConstRecHitPointer oriHit0 = ip->inner();
     SeedingHitSet::ConstRecHitPointer oriHit1 = ip->outer();
 
     HitOwnPtr hit0(*oriHit0);
     HitOwnPtr hit1(*oriHit1);
     GlobalPoint gp0 = hit0->globalPosition();
     GlobalPoint gp1 = hit1->globalPosition();
 
     bool debugPair = debug && ip->inner()->rawId()==debug_Id0 && ip->outer()->rawId()==debug_Id1;
 
     if (debugPair) {
       cout << endl << endl
        << "found new pair with ids "<<debug_Id0<<" "<<debug_Id1<<" with pos: " << gp0 << " " << gp1 
            << endl;
     }
 
     if (refitHits) {
 
       TrajectoryStateOnSurface tsos0, tsos1;
       assert(!hit0.isOwn()); assert(!hit1.isOwn());
       refit2Hits(hit0,hit1,tsos0,tsos1,region,nomField,debugPair);
       assert(hit0.isOwn()); assert(hit1.isOwn());
 
       //fixme add pixels
       bool passFilterHit0 = true;
       if (//hit0->geographicalId().subdetId() > 2
       hit0->geographicalId().subdetId()==SiStripDetId::TIB 
       || hit0->geographicalId().subdetId()==SiStripDetId::TID
       //|| hit0->geographicalId().subdetId()==SiStripDetId::TOB
       //|| hit0->geographicalId().subdetId()==SiStripDetId::TEC
       ) {   
     const std::type_info &tid = typeid(*hit0->hit());
     if (tid == typeid(SiStripMatchedRecHit2D)) {
       const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D *>(hit0->hit());
       if (filter->isCompatible(DetId(matchedHit->monoId()), matchedHit->monoCluster(), tsos0.localMomentum())==0 ||
           filter->isCompatible(DetId(matchedHit->stereoId()), matchedHit->stereoCluster(), tsos0.localMomentum())==0) passFilterHit0 = false;
     } else if (tid == typeid(SiStripRecHit2D)) {
       const SiStripRecHit2D* recHit = dynamic_cast<const SiStripRecHit2D *>(hit0->hit());
       if (filter->isCompatible(*recHit, tsos0.localMomentum())==0) passFilterHit0 = false;
     } else if (tid == typeid(ProjectedSiStripRecHit2D)) {
       const ProjectedSiStripRecHit2D* precHit = dynamic_cast<const ProjectedSiStripRecHit2D *>(hit0->hit());
       if (filter->isCompatible(precHit->originalHit(), tsos0.localMomentum())==0) passFilterHit0 = false;   //FIXME
     }
       }
       if (debugPair&&!passFilterHit0)  cout << "hit0 did not pass cluster shape filter" << endl;
       if (!passFilterHit0) continue;
       bool passFilterHit1 = true;
       if (//hit1->geographicalId().subdetId() > 2
       hit1->geographicalId().subdetId()==SiStripDetId::TIB 
       || hit1->geographicalId().subdetId()==SiStripDetId::TID
       //|| hit1->geographicalId().subdetId()==SiStripDetId::TOB
       //|| hit1->geographicalId().subdetId()==SiStripDetId::TEC
       ) {   
     const std::type_info &tid = typeid(*hit1->hit());
     if (tid == typeid(SiStripMatchedRecHit2D)) {
       const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D *>(hit1->hit());
       if (filter->isCompatible(DetId(matchedHit->monoId()), matchedHit->monoCluster(), tsos1.localMomentum())==0 ||
           filter->isCompatible(DetId(matchedHit->stereoId()), matchedHit->stereoCluster(), tsos1.localMomentum())==0) passFilterHit1 = false;
     } else if (tid == typeid(SiStripRecHit2D)) {
       const SiStripRecHit2D* recHit = dynamic_cast<const SiStripRecHit2D *>(hit1->hit());
       if (filter->isCompatible(*recHit, tsos1.localMomentum())==0) passFilterHit1 = false;
     } else if (tid == typeid(ProjectedSiStripRecHit2D)) {
       const ProjectedSiStripRecHit2D* precHit = dynamic_cast<const ProjectedSiStripRecHit2D *>(hit1->hit());
       if (filter->isCompatible(precHit->originalHit(), tsos1.localMomentum())==0) passFilterHit1 = false;  //FIXME
     }
       }
       if (debugPair&&!passFilterHit1)  cout << "hit1 did not pass cluster shape filter" << endl;
       if (!passFilterHit1) continue;
 
     } else {
       // not refit clone anyhow
       hit0.reset((BaseTrackerRecHit *)hit0->clone());
       hit1.reset((BaseTrackerRecHit *)hit1->clone());
     }
 
     //gc: create the RZ line for the pair
     SimpleLineRZ line(PixelRecoPointRZ(gp0.perp(),gp0.z()), PixelRecoPointRZ(gp1.perp(),gp1.z()));
     ThirdHitPredictionFromCircle predictionRPhi(gp0, gp1, extraHitRPhitolerance);
 
     //gc: this is the curvature of the two hits assuming the region
     Range pairCurvature = predictionRPhi.curvature(region.originRBound());
     //gc: intersect not only returns a bool but may change pairCurvature to intersection with curv
     if (!intersect(pairCurvature, Range(-curv, curv))) {
       if (debugPair) std::cout << "curvature cut: curv=" << curv 
                    << " gc=(" << pairCurvature.first << ", " << pairCurvature.second << ")" << std::endl;
       continue;
     }
 
     //gc: loop over all third layers compatible with the pair
     for(int il = 0; (il < size) & (!usePair); il++) {
 
       if (debugPair) 
     cout << "cosider layer: " << theLayers[il].name() << " for this pair. Location: " << theLayers[il].detLayer()->location() << endl;
 
       if (hitTree[il].empty()) {
     if (debugPair) {
       cout << "empty hitTree" << endl;
     }
     continue; // Don't bother if no hits
       }
 
       cacheHitPointer bestL2;
       float chi2FromThisLayer = std::numeric_limits<float>::max();
 
       const DetLayer *layer = theLayers[il].detLayer();
       bool barrelLayer = layer->location() == GeomDetEnumerators::barrel;
 
       LayerRZPredictions &predRZ = mapPred.find(theLayers[il].name())->second;
       predRZ.line.initPropagator(&line);
       
       //gc: this takes the z at R-thick/2 and R+thick/2 according to 
       //    the line from the two points and the adds the extra tolerance
       Range rzRange = predRZ.line();
 
       if (rzRange.first >= rzRange.second) {
     if (debugPair) {
       cout << "rzRange empty" << endl;
     }
         continue;
       }
       //gc: check that rzRange is compatible with detector bounds
       //    note that intersect may change rzRange to intersection with bounds
       if (!intersect(rzRange, predRZ.line.detSize())) {// theDetSize = Range(-maxZ, maxZ); 
     if (debugPair) {
       cout << "rzRange and detector do not intersect" << endl;
     }
     continue;
       }
       Range radius = barrelLayer ? predRZ.line.detRange() : rzRange;
 
       //gc: define the phi range of the hits
       Range phiRange;
       if (useFixedPreFiltering) { 
     //gc: in this case it takes as range the phi of the outer 
     //    hit +/- the phiPreFiltering value from cfg
         float phi0 = ip->outer()->globalPosition().phi();
         phiRange = Range(phi0 - dphi, phi0 + dphi);
       } else {  
     //gc: predictionRPhi uses the cosine rule to find the phi of the 3rd point at radius, assuming the pairCurvature range [-c,+c]
     if (pairCurvature.first<0. && pairCurvature.second<0.) {
       float phi12 = predictionRPhi.phi(pairCurvature.first,radius.second);
       float phi21 = predictionRPhi.phi(pairCurvature.second,radius.first);
       while(unlikely(phi12 <  phi21)) phi12 += float(2. * M_PI); 
       phiRange = Range(phi21,phi12);
     } else if (pairCurvature.first>=0. && pairCurvature.second>=0.) {
       float phi11 = predictionRPhi.phi(pairCurvature.first,radius.first);
       float phi22 = predictionRPhi.phi(pairCurvature.second,radius.second);
       while(unlikely(phi11 <  phi22)) phi11 += float(2. * M_PI); 
       phiRange = Range(phi22,phi11);
     } else {
       float phi12 = predictionRPhi.phi(pairCurvature.first,radius.second);
       float phi22 = predictionRPhi.phi(pairCurvature.second,radius.second);
       while(unlikely(phi12 <  phi22)) phi12 += float(2. * M_PI); 
       phiRange = Range(phi22,phi12);
     }
       }
       
       //gc: this is the place where hits in the compatible region are put in the foundNodes
       typedef RecHitsSortedInPhi::Hit Hit;
       foundNodes.clear(); // Now recover hits in bounding box...
       // This needs range -twoPi to +twoPi to work
       float prmin=phiRange.min(), prmax=phiRange.max(); //get contiguous range
       if ((prmax-prmin) > Geom::twoPi()) { 
     prmax=Geom::pi(); prmin = -Geom::pi();
       } else {
     while (prmax>maxphi) { prmin -= Geom::twoPi(); prmax -= Geom::twoPi();}
     while (prmin<minphi) { prmin += Geom::twoPi(); prmax += Geom::twoPi();}
       }
       
       if (debugPair) cout << "defining kd tree box" << endl;
 
       if (barrelLayer) {
     KDTreeBox phiZ(prmin-extraPhiKDBox, prmax+extraPhiKDBox,
                rzRange.min()-fnSigmaRZ*rzError[il]-extraZKDBox,
                rzRange.max()+fnSigmaRZ*rzError[il]+extraZKDBox);
     hitTree[il].search(phiZ, foundNodes);
 
     if (debugPair) cout << "kd tree box bounds, phi: " << prmin-extraPhiKDBox <<","<< prmax+extraPhiKDBox
                 << " z: "<< rzRange.min()-fnSigmaRZ*rzError[il]-extraZKDBox <<","<<rzRange.max()+fnSigmaRZ*rzError[il]+extraZKDBox
                 << " rzRange: " << rzRange.min() <<","<<rzRange.max()
                 << endl;
 
       } else {
     KDTreeBox phiR(prmin-extraPhiKDBox, prmax+extraPhiKDBox,
                rzRange.min()-fnSigmaRZ*rzError[il]-extraRKDBox,
                rzRange.max()+fnSigmaRZ*rzError[il]+extraRKDBox);
     hitTree[il].search(phiR, foundNodes);
 
     if (debugPair) cout << "kd tree box bounds, phi: " << prmin-extraPhiKDBox <<","<< prmax+extraPhiKDBox
                 << " r: "<< rzRange.min()-fnSigmaRZ*rzError[il]-extraRKDBox <<","<<rzRange.max()+fnSigmaRZ*rzError[il]+extraRKDBox
                 << " rzRange: " << rzRange.min() <<","<<rzRange.max()
                 << endl;
       }
 
       if (debugPair) cout << "kd tree box size: " << foundNodes.size() << endl;
       
 
       //gc: now we loop over the hits in the box for this layer
       for (std::vector<RecHitsSortedInPhi::HitIter>::iterator ih = foundNodes.begin();
        ih !=foundNodes.end() && !usePair; ++ih) {
 
     if (debugPair) std::cout << endl << "triplet candidate" << std::endl;
 
     const RecHitsSortedInPhi::HitIter KDdata = *ih;
 
     SeedingHitSet::ConstRecHitPointer oriHit2 = KDdata->hit();
     cacheHitPointer hit2;
 
     if (refitHits) {//fixme
 
       //fitting all 3 hits takes too much time... do it quickly only for 3rd hit
       GlobalVector initMomentum(oriHit2->globalPosition() - gp1);
       initMomentum *= (1./initMomentum.perp()); //set pT=1
       GlobalTrajectoryParameters kine = GlobalTrajectoryParameters(oriHit2->globalPosition(), initMomentum, 1, &*bfield);
       TrajectoryStateOnSurface state(kine,*oriHit2->surface());
       hit2.reset((SeedingHitSet::RecHitPointer)(cloner(*oriHit2,state)));
 
       //fixme add pixels
       bool passFilterHit2 = true;
       if (hit2->geographicalId().subdetId()==SiStripDetId::TIB 
           || hit2->geographicalId().subdetId()==SiStripDetId::TID
           // || hit2->geographicalId().subdetId()==SiStripDetId::TOB
           // || hit2->geographicalId().subdetId()==SiStripDetId::TEC
           ) {
         const std::type_info &tid = typeid(*hit2->hit());
         if (tid == typeid(SiStripMatchedRecHit2D)) {
           const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D *>(hit2->hit());
           if (filter->isCompatible(DetId(matchedHit->monoId()), matchedHit->monoCluster(), initMomentum)==0 ||
           filter->isCompatible(DetId(matchedHit->stereoId()), matchedHit->stereoCluster(), initMomentum)==0) passFilterHit2 = false;
         } else if (tid == typeid(SiStripRecHit2D)) {
           const SiStripRecHit2D* recHit = dynamic_cast<const SiStripRecHit2D *>(hit2->hit());
           if (filter->isCompatible(*recHit, initMomentum)==0) passFilterHit2 = false;
         } else if (tid == typeid(ProjectedSiStripRecHit2D)) {
           const ProjectedSiStripRecHit2D* precHit = dynamic_cast<const ProjectedSiStripRecHit2D *>(hit2->hit());
           if (filter->isCompatible(precHit->originalHit(), initMomentum)==0) passFilterHit2 = false;
         }
       }
       if (debugPair&&!passFilterHit2)  cout << "hit2 did not pass cluster shape filter" << endl;
       if (!passFilterHit2) continue;
       
       // fitting all 3 hits takes too much time :-(
       // TrajectoryStateOnSurface tsos0, tsos1, tsos2;
       // refit3Hits(hit0,hit1,hit2,tsos0,tsos1,tsos2,nomField,debugPair);
       // if (hit2->geographicalId().subdetId()==SiStripDetId::TIB 
       //     // || hit2->geographicalId().subdetId()==SiStripDetId::TOB
       //     // || hit2->geographicalId().subdetId()==SiStripDetId::TID
       //     // || hit2->geographicalId().subdetId()==SiStripDetId::TEC
       //     ) {
       //   const std::type_info &tid = typeid(*hit2->hit());
       //   if (tid == typeid(SiStripMatchedRecHit2D)) {
       //     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D *>(hit2->hit());
       //     if (filter->isCompatible(DetId(matchedHit->monoId()), matchedHit->monoCluster(), tsos2.localMomentum())==0 ||
       //          filter->isCompatible(DetId(matchedHit->stereoId()), matchedHit->stereoCluster(), tsos2.localMomentum())==0) continue;
       //   } else if (tid == typeid(SiStripRecHit2D)) {
       //     const SiStripRecHit2D* recHit = dynamic_cast<const SiStripRecHit2D *>(hit2->hit());
       //     if (filter->isCompatible(*recHit, tsos2.localMomentum())==0) continue;
       //   } else if (tid == typeid(ProjectedSiStripRecHit2D)) {
       //     const ProjectedSiStripRecHit2D* precHit = dynamic_cast<const ProjectedSiStripRecHit2D *>(hit2->hit());
       //     if (filter->isCompatible(precHit->originalHit(), tsos2.localMomentum())==0) continue;;
       //   }
       // }
 
     } else {
       // not refit clone anyhow
       hit2.reset((BaseTrackerRecHit*)oriHit2->clone());
     }
 
     //gc: add the chi2 cut
     vector<GlobalPoint> gp(3);
     vector<GlobalError> ge(3);
     vector<bool> bl(3);
     gp[0] = hit0->globalPosition();
     ge[0] = hit0->globalPositionError();
     int subid0 = hit0->geographicalId().subdetId();
     bl[0] = (subid0 == StripSubdetector::TIB || subid0 == StripSubdetector::TOB || subid0 == (int) PixelSubdetector::PixelBarrel);
     gp[1] = hit1->globalPosition();
     ge[1] = hit1->globalPositionError();
     int subid1 = hit1->geographicalId().subdetId();
     bl[1] = (subid1 == StripSubdetector::TIB || subid1 == StripSubdetector::TOB || subid1 == (int) PixelSubdetector::PixelBarrel);
     gp[2] = hit2->globalPosition();
     ge[2] = hit2->globalPositionError();
     int subid2 = hit2->geographicalId().subdetId();
     bl[2] = (subid2 == StripSubdetector::TIB || subid2 == StripSubdetector::TOB || subid2 == (int) PixelSubdetector::PixelBarrel);
     RZLine rzLine(gp,ge,bl);
     float  cottheta, intercept, covss, covii, covsi;
     rzLine.fit(cottheta, intercept, covss, covii, covsi);
     float chi2 = rzLine.chi2(cottheta, intercept);
 
     bool debugTriplet = debugPair && hit2->rawId()==debug_Id2;
     if (debugTriplet) {
       std::cout << endl << "triplet candidate in debug id" << std::endl;
       cout << "hit in id="<<hit2->rawId()<<" (from KDTree) with pos: " << KDdata->hit()->globalPosition()
            << " refitted: " << hit2->globalPosition() 
            << " chi2: " << chi2
            << endl;
       //cout << state << endl;
     }
     // should fix nan
     if ( (chi2 > maxChi2) | edm::isNotFinite(chi2) ) continue; 
 
     if (chi2VsPtCut) {
 
       FastCircle theCircle(hit2->globalPosition(),hit1->globalPosition(),hit0->globalPosition());
       float tesla0 = 0.1*nomField;
       float rho = theCircle.rho();
       float cm2GeV = 0.01 * 0.3*tesla0;
       float pt = cm2GeV * rho;
       if (debugTriplet) {
         std::cout << "triplet pT=" << pt << std::endl;
       }
       if (pt<region.ptMin()) continue;
       
       if (chi2_cuts.size()==4) {
         if ( ( pt>pt_interv[0] && pt<=pt_interv[1] && chi2 > chi2_cuts[0] ) ||
          ( pt>pt_interv[1] && pt<=pt_interv[2] && chi2 > chi2_cuts[1] ) ||
          ( pt>pt_interv[2] && pt<=pt_interv[3] && chi2 > chi2_cuts[2] ) ||
          ( pt>pt_interv[3] && chi2 > chi2_cuts[3] ) ) continue;     
       }
       
       // apparently this takes too much time...
       //      if (chi2_cuts.size()>1) {     
       //        int ncuts = chi2_cuts.size();
       //        if ( pt<=pt_interv[0] && chi2 > chi2_cuts[0] ) continue; 
       //        bool pass = true;
       //        for (int icut=1; icut<ncuts-1; icut++){
       //          if ( pt>pt_interv[icut-1] && pt<=pt_interv[icut] && chi2 > chi2_cuts[icut] ) pass=false;
       //        }
       //        if (!pass) continue;
       //        if ( pt>pt_interv[ncuts-2] && chi2 > chi2_cuts[ncuts-1] ) continue;         
       //        if (debug && hit0->rawId()==debug_Id0 && hit1->rawId()==debug_Id1 && hit2->rawId()==debug_Id2) {
       //          std::cout << "triplet passed chi2 vs pt cut" << std::endl;
       //        }
       //      } 
       
     }
     
     if (theMaxElement!=0 && result.size() >= theMaxElement) {
       result.clear();
       edm::LogError("TooManyTriplets")<<" number of triples exceed maximum. no triplets produced.";
       return;
     }
     if (debugPair) std::cout << "triplet made" << std::endl;
     //result.push_back(SeedingHitSet(hit0, hit1, hit2));
     /* no refit so keep only hit2
     assert(tripletFromThisLayer.empty());
     assert(hit0.isOwn()); assert(hit1.isOwn());assert(hit2.isOwn());
     tripletFromThisLayer.emplace_back(std::move(hit0));
     tripletFromThisLayer.emplace_back(std::move(hit1));
     tripletFromThisLayer.emplace_back(std::move(hit2));
     assert(hit0.isEmpty()); assert(hit1.isEmpty());assert(hit2.isEmpty());
     */
     bestL2 = std::move(hit2);
     chi2FromThisLayer = chi2;
     foundTripletsFromPair++;
     if (foundTripletsFromPair>=2) {
       usePair=true;
       if (debugPair) 
            std::cout << "using pair" << std::endl;
       break;
     }
       }//loop over hits in KDTree
 
       if (usePair) break;
       else {
     //if there is one triplet in more than one layer, try picking the one with best chi2
     if (chi2FromThisLayer<minChi2) {
       bestH2 = std::move(bestL2);
       minChi2 = chi2FromThisLayer;
     }
     /*
     else {
       if (!bestH2 && foundTripletsFromPair>0)
         std::cout << "what?? " <<  minChi2 << ' '  << chi2FromThisLayer << std::endl;
     }
     */
       }
 
     }//loop over layers
 
     if (foundTripletsFromPair==0) continue;
 
     //push back only (max) once per pair
     if (debugPair) std::cout << "Done seed #" << result.size() << std::endl;
     if (usePair) result.push_back(SeedingHitSet(ip->inner(), ip->outer())); 
     else { 
       assert(1==foundTripletsFromPair);
       assert(bestH2);
       result.emplace_back(&*hit0,&*hit1,&*bestH2); 
       assert(hit0.isOwn()); assert(hit1.isOwn());
       cache.emplace_back(const_cast<BaseTrackerRecHit*>(hit0.release()));
       cache.emplace_back(const_cast<BaseTrackerRecHit*>(hit1.release()));
       cache.emplace_back(std::move(bestH2));
       assert(hit0.empty()); assert(hit1.empty());assert(!bestH2);
     }
     // std::cout << (usePair ? "pair " : "triplet ") << minChi2 <<' ' << cache.size() << std::endl;  
 
 
   }//loop over pairs
   if (debug) {
     std::cout << "triplet size=" << result.size() << std::endl;
   }
 }

void MultiHitGeneratorFromChi2::init	(	const HitPairGenerator &	pairs,
		LayerCacheType *	layerCache
	)

overridevirtual

Implements MultiHitGeneratorFromPairAndLayers.

Definition at line 97 of file MultiHitGeneratorFromChi2.cc.

References HitPairGenerator::clone(), theLayerCache, and thePairGenerator.

 {
   thePairGenerator = pairs.clone();
   theLayerCache = layerCache;
 }

void MultiHitGeneratorFromChi2::initES ( const edm::EventSetup & es )

overridevirtual

Implements MultiHitGeneratorFromPairAndLayers.

Definition at line 104 of file MultiHitGeneratorFromChi2.cc.

References bfield, builder, builderName_, cloner, filter, filterName_, edm::EventSetup::get(), mfName_, nomField, MagneticField::nominalValue(), edm::ESHandle< class >::product(), and useSimpleMF_.

 {
 
   edm::ESHandle<MagneticField> bfield_h;
   if (useSimpleMF_) es.get<IdealMagneticFieldRecord>().get(mfName_, bfield_h);  
   else es.get<IdealMagneticFieldRecord>().get(bfield_h);
   bfield = bfield_h.product();
   nomField = bfield->nominalValue();
 
   edm::ESHandle<ClusterShapeHitFilter> filterHandle_;
   es.get<CkfComponentsRecord>().get(filterName_, filterHandle_);
   filter = filterHandle_.product();
 
   edm::ESHandle<TransientTrackingRecHitBuilder> builderH;
   es.get<TransientRecHitRecord>().get(builderName_, builderH);
   builder = (TkTransientTrackingRecHitBuilder const *)(builderH.product());
   cloner = (*builder).cloner();
 }

std::pair< float, float > MultiHitGeneratorFromChi2::mergePhiRanges	(	const std::pair< float, float > &	r1,
		const std::pair< float, float > &	r2
	)		const

private

Definition at line 655 of file MultiHitGeneratorFromChi2.cc.

References M_PI, bookConverter::max, and min().

 { float r2Min = r2.first;
   float r2Max = r2.second;
   while (r1.first - r2Min > +M_PI) r2Min += 2. * M_PI, r2Max += 2. * M_PI;
   while (r1.first - r2Min < -M_PI) r2Min -= 2. * M_PI, r2Max -= 2. * M_PI;
   //std::cout << "mergePhiRanges " << fabs(r1.first-r2Min) << " " <<  fabs(r1.second-r2Max) << endl;
   return std::make_pair(min(r1.first, r2Min), max(r1.second, r2Max));
 }

const HitPairGenerator& MultiHitGeneratorFromChi2::pairGenerator ( ) const

inline

Definition at line 46 of file MultiHitGeneratorFromChi2.h.

46 { return *thePairGenerator; }

MultiHitGeneratorFromChi2::thePairGenerator

HitPairGenerator * thePairGenerator

Definition: MultiHitGeneratorFromChi2.h:70

void MultiHitGeneratorFromChi2::refit2Hits	(	HitOwnPtr &	hit0,
		HitOwnPtr &	hit1,
		TrajectoryStateOnSurface &	tsos0,
		TrajectoryStateOnSurface &	tsos1,
		const TrackingRegion &	region,
		float	nomField,
		bool	isDebug
	)

private

Definition at line 665 of file MultiHitGeneratorFromChi2.cc.

References bfield, cloner, gather_cfg::cout, PV3DBase< T, PVType, FrameType >::eta(), TrajectoryStateOnSurface::globalMomentum(), nomField, TrackingRegion::origin(), p1, p2, PV3DBase< T, PVType, FrameType >::perp(), PV3DBase< T, PVType, FrameType >::phi(), EnergyCorrector::pt, lumiQueryAPI::q, mayown_ptr< T, N >::reset(), rho, FastCircle::rho(), funct::tan(), PV3DBase< T, PVType, FrameType >::theta(), x, PV3DBase< T, PVType, FrameType >::x(), FastCircle::x0(), y, PV3DBase< T, PVType, FrameType >::y(), and FastCircle::y0().

Referenced by hitSets().

                                                                                    {
 
   //these need to be sorted in R
   GlobalPoint gp0 = region.origin();
   GlobalPoint gp1 = hit1->globalPosition();
   GlobalPoint gp2 = hit2->globalPosition();
 
   if (isDebug) {
     cout << "positions before refitting: " << hit1->globalPosition() << " " << hit2->globalPosition() <<endl;
   }
 
   FastCircle theCircle(gp2,gp1,gp0);
   GlobalPoint cc(theCircle.x0(),theCircle.y0(),0);
   float tesla0 = 0.1*nomField;
   float rho = theCircle.rho();
   float cm2GeV = 0.01 * 0.3*tesla0;
   float pt = cm2GeV * rho;
 
   GlobalVector vec20 = gp2-gp0;
   //if (isDebug) { cout << "vec20.eta=" << vec20.eta() << endl; }
 
   GlobalVector p0( gp0.y()-cc.y(), -gp0.x()+cc.x(), 0. );
   p0 = p0*pt/p0.perp();
   GlobalVector p1( gp1.y()-cc.y(), -gp1.x()+cc.x(), 0. );
   p1 = p1*pt/p1.perp();
   GlobalVector p2( gp2.y()-cc.y(), -gp2.x()+cc.x(), 0. );
   p2 = p2*pt/p2.perp();
 
   //check sign according to scalar product
   if ( (p0.x()*(gp1.x()-gp0.x())+p0.y()*(gp1.y()-gp0.y()) ) < 0 ) {
     p0*=-1.;
     p1*=-1.;
     p2*=-1.;
   }
 
   //now set z component
   p0 = GlobalVector(p0.x(),p0.y(),p0.perp()/tan(vec20.theta()));
   p1 = GlobalVector(p1.x(),p1.y(),p1.perp()/tan(vec20.theta()));
   p2 = GlobalVector(p2.x(),p2.y(),p2.perp()/tan(vec20.theta()));
 
   //get charge from vectorial product
   TrackCharge q = 1;
   if ((gp1-cc).x()*p1.y() - (gp1-cc).y()*p1.x() > 0) q =-q;
 
   GlobalTrajectoryParameters kine1 = GlobalTrajectoryParameters(gp1, p1, q, &*bfield);
   state1 = TrajectoryStateOnSurface(kine1,*hit1->surface());
   hit1.reset((SeedingHitSet::RecHitPointer)(cloner(*hit1,state1)));
 
   GlobalTrajectoryParameters kine2 = GlobalTrajectoryParameters(gp2, p2, q, &*bfield);
   state2 = TrajectoryStateOnSurface(kine2,*hit2->surface());
   hit2.reset((SeedingHitSet::RecHitPointer)(cloner(*hit2,state2)));
 
   if (isDebug) {
     cout << "charge=" << q << endl;
     cout << "state1 pt=" << state1.globalMomentum().perp() << " eta=" << state1.globalMomentum().eta()  << " phi=" << state1.globalMomentum().phi() << endl;
     cout << "state2 pt=" << state2.globalMomentum().perp() << " eta=" << state2.globalMomentum().eta()  << " phi=" << state2.globalMomentum().phi() << endl;
     cout << "positions after refitting: " << hit1->globalPosition() << " " << hit2->globalPosition() <<endl;
   }
 
 }

void MultiHitGeneratorFromChi2::setSeedingLayers	(	SeedingLayerSetsHits::SeedingLayerSet	pairLayers,
		std::vector< SeedingLayerSetsHits::SeedingLayer >	thirdLayers
	)

overridevirtual

Implements MultiHitGeneratorFromPairAndLayers.

Definition at line 123 of file MultiHitGeneratorFromChi2.cc.

References HitPairGenerator::setSeedingLayers(), theLayers, and thePairGenerator.

                                                                                                           {
   thePairGenerator->setSeedingLayers(pairLayers);
   theLayers = thirdLayers;
 }