#include <PixelTripletLargeTipGenerator.h>

Inheritance diagram for PixelTripletLargeTipGenerator:

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

void	hitTriplets (const TrackingRegion &region, OrderedHitTriplets &trs, const edm::Event &ev, const edm::EventSetup &es, const HitDoublets &doublets, const std::vector< SeedingLayerSetsHits::SeedingLayer > &thirdLayers, std::vector< int > *tripletLastLayerIndex, LayerCacheType &layerCache)

void	hitTriplets (const TrackingRegion &region, OrderedHitTriplets &result, const edm::EventSetup &es, const HitDoublets &doublets, const RecHitsSortedInPhi *thirdHitMap, const std::vector< const DetLayer > &thirdLayerDetLayer, const int nThirdLayers) override

void	hitTriplets (const TrackingRegion &region, OrderedHitTriplets &result, const edm::EventSetup &es, const HitDoublets &doublets, const RecHitsSortedInPhi *thirdHitMap, const std::vector< const DetLayer > &thirdLayerDetLayer, const int nThirdLayers, std::vector< int > *tripletLastLayerIndex)

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

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

	~PixelTripletLargeTipGenerator () override

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 ()

Static Public Member Functions
static void	fillDescriptions (edm::ParameterSetDescription &desc)

static const char *	fillDescriptionsLabel ()

Static Public Member Functions inherited from HitTripletGeneratorFromPairAndLayers
static void	fillDescriptions (edm::ParameterSetDescription &desc)

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 18 of file PixelTripletLargeTipGenerator.h.

Member Typedef Documentation

typedef CombinedHitTripletGenerator::LayerCacheType PixelTripletLargeTipGenerator::LayerCacheType

private

Definition at line 19 of file PixelTripletLargeTipGenerator.h.

Constructor & Destructor Documentation

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

inline

Definition at line 22 of file PixelTripletLargeTipGenerator.h.

References looper::cfg, fillDescriptions(), and ~PixelTripletLargeTipGenerator().

23 : PixelTripletLargeTipGenerator(cfg, iC) {}

PixelTripletLargeTipGenerator::PixelTripletLargeTipGenerator

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

Definition: PixelTripletLargeTipGenerator.h:22

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

Definition at line 47 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 ( )

override

Definition at line 56 of file PixelTripletLargeTipGenerator.cc.

Referenced by PixelTripletLargeTipGenerator().

56 {}

Member Function Documentation

void PixelTripletLargeTipGenerator::fillDescriptions ( edm::ParameterSetDescription & desc )

static

Definition at line 58 of file PixelTripletLargeTipGenerator.cc.

References edm::ParameterSetDescription::add(), HitTripletGeneratorFromPairAndLayers::fillDescriptions(), FastTimerService_cff::range, and edm::second().

Referenced by PixelTripletLargeTipGenerator().

                                                                                      {
   HitTripletGeneratorFromPairAndLayers::fillDescriptions(desc);
   // Defaults for the extraHit*tolerance are set to 0 here since that
   // was already the case in practice in all offline occurrances.
   desc.add<double>("extraHitRPhitolerance", 0);  // default in old python was 0.032
   desc.add<double>("extraHitRZtolerance", 0);    // default in old python was 0.037
   desc.add<bool>("useMultScattering", true);
   desc.add<bool>("useBending", true);
   desc.add<bool>("useFixedPreFiltering", false);
   desc.add<double>("phiPreFiltering", 0.3);
 }

static const char* PixelTripletLargeTipGenerator::fillDescriptionsLabel ( )

inlinestatic

Definition at line 29 of file PixelTripletLargeTipGenerator.h.

References HLT_2018_cff::doublets, ev, hitTriplets(), HLT_2018_cff::region, and mps_fire::result.

29 { return "pixelTripletLargeTip"; }

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

overridevirtual

Implements HitTripletGeneratorFromPairAndLayers.

Definition at line 82 of file PixelTripletLargeTipGenerator.cc.

References HLT_2018_cff::doublets, HitTripletGeneratorFromPairAndLayers::theLayerCache, and HitTripletGeneratorFromPairAndLayers::thePairGenerator.

Referenced by fillDescriptionsLabel(), and hitTriplets().

                                                                                                                 {
   auto const& doublets = thePairGenerator->doublets(region, ev, es, pairLayers);
 
   if (doublets.empty())
     return;
 
   assert(theLayerCache);
   hitTriplets(region, result, ev, es, doublets, thirdLayers, nullptr, *theLayerCache);
 }

void PixelTripletLargeTipGenerator::hitTriplets	(	const TrackingRegion &	region,
		OrderedHitTriplets &	trs,
		const edm::Event &	ev,
		const edm::EventSetup &	es,
		const HitDoublets &	doublets,
		const std::vector< SeedingLayerSetsHits::SeedingLayer > &	thirdLayers,
		std::vector< int > *	tripletLastLayerIndex,
		LayerCacheType &	layerCache
	)

Definition at line 97 of file PixelTripletLargeTipGenerator.cc.

References hitTriplets(), and findQualityFiles::size.

                                                                             {
   int size = thirdLayers.size();
   const RecHitsSortedInPhi* thirdHitMap[size];
   vector<const DetLayer*> thirdLayerDetLayer(size, nullptr);
   for (int il = 0; il < size; ++il) {
     thirdHitMap[il] = &layerCache(thirdLayers[il], region, es);
     thirdLayerDetLayer[il] = thirdLayers[il].detLayer();
   }
   hitTriplets(region, result, es, doublets, thirdHitMap, thirdLayerDetLayer, size, tripletLastLayerIndex);
 }

void PixelTripletLargeTipGenerator::hitTriplets	(	const TrackingRegion &	region,
		OrderedHitTriplets &	result,
		const edm::EventSetup &	es,
		const HitDoublets &	doublets,
		const RecHitsSortedInPhi **	thirdHitMap,
		const std::vector< const DetLayer * > &	thirdLayerDetLayer,
		const int	nThirdLayers
	)

overridevirtual

Implements HitTripletGeneratorFromPairAndLayers.

Definition at line 115 of file PixelTripletLargeTipGenerator.cc.

References hitTriplets().

                                                                         {
   hitTriplets(region, result, es, doublets, thirdHitMap, thirdLayerDetLayer, nThirdLayers, nullptr);
 }

void PixelTripletLargeTipGenerator::hitTriplets	(	const TrackingRegion &	region,
		OrderedHitTriplets &	result,
		const edm::EventSetup &	es,
		const HitDoublets &	doublets,
		const RecHitsSortedInPhi **	thirdHitMap,
		const std::vector< const DetLayer * > &	thirdLayerDetLayer,
		const int	nThirdLayers,
		std::vector< int > *	tripletLastLayerIndex
	)

Definition at line 125 of file PixelTripletLargeTipGenerator.cc.

References angle(), checkPhiInRange(), constexpr, pfMETCorrectionType0_cfi::correction, PixelRecoUtilities::curvature(), declareDynArray, HitDoublets::detLayer(), dphi, relativeConstraints::empty, runTheMatrix::err, extraHitRPhitolerance, extraHitRZtolerance, f, dqmdumpme::first, dqmMemoryStats::float, fnSigmaRZ, Geom::ftwoPi(), edm::EventSetup::get(), HitDoublets::hit(), hfClusterShapes_cfi::hits, mps_fire::i, HitDoublets::inner, DetLayer::isBarrel(), mps_splice::line, M_PI, SiStripPI::max, ThirdHitPredictionFromCircle::HelixRZ::maxCurvature(), min(), normalizedPhi(), nSigmaPhi, nSigmaRZ, PixelRecoLineRZ::origin(), TrackingRegion::originRBound(), HitDoublets::outer, p3, PV3DBase< T, PVType, FrameType >::perp(), HitDoublets::phi(), GeometricSearchDet::position(), edm::ESHandle< T >::product(), proxim(), TrackingRegion::ptMin(), Basic2DVector< T >::r(), CosmicsPD_Skims::radius, edm::second(), DetLayer::seqNum(), OrderedHitTriplets::size(), HitDoublets::size(), std::swap(), HitTripletGeneratorFromPairAndLayers::theMaxElement, PbPb_ZMuSkimMuonDPG_cff::tracker, useBend, useFixedPreFiltering, useMScat, findQualityFiles::v, PV3DBase< T, PVType, FrameType >::x(), HitDoublets::x(), hybridSuperClusters_cfi::xi, PV3DBase< T, PVType, FrameType >::y(), HitDoublets::y(), PV3DBase< T, PVType, FrameType >::z(), and HitDoublets::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 outSeq = doublets.detLayer(HitDoublets::outer)->seqNum();
 
   using NodeInfo = KDTreeNodeInfo<unsigned int, 2>;
   std::vector<NodeInfo> layerTree;       // re-used throughout
   std::vector<unsigned int> foundNodes;  // re-used throughout
   foundNodes.reserve(100);
 
   declareDynArray(KDTreeLinkerAlgo<unsigned int>, nThirdLayers, hitTree);
   declareDynArray(LayerRZPredictions, nThirdLayers, mapPred);
 
   float rzError[nThirdLayers];  //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
 
   for (int il = 0; il < nThirdLayers; il++) {
     auto const& hits = *thirdHitMap[il];
 
     const DetLayer* layer = thirdLayerDetLayer[il];
     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),
                            *thirdLayerDetLayer[il],
                            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 < nThirdLayers; il++) {
       if (hitTree[il].empty())
         continue;  // Don't bother if no hits
       const DetLayer* layer = thirdLayerDetLayer[il];
       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;
 
       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;
         // limit error to 90 degree
         constexpr float maxPhiErr = 0.5 * M_PI;
         float phiErr = nSigmaPhi * hits.drphi[KDdata] * ir;
         phiErr = std::min(maxPhiErr, phiErr);
         if (!checkPhiInRange(p3_phi, rangeRPhi.first * ir - phiErr, rangeRPhi.second * ir + phiErr, maxPhiErr))
           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();
           if (tripletLastLayerIndex)
             tripletLastLayerIndex->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());
         // to bookkeep the triplets and 3rd layers in triplet EDProducer
         if (tripletLastLayerIndex)
           tripletLastLayerIndex->push_back(il);
       }
     }
   }
   // std::cout << "found triplets " << result.size() << std::endl;
 }