CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CAHitQuadrupletGenerator.cc
Go to the documentation of this file.
2 
4 
6 #include "LayerQuadruplets.h"
11 
15 
16 
18 
19 
20 
21 #include "CellularAutomaton.h"
22 
24 
26 
27 #include <functional>
28 
29 namespace
30 {
31 
32  template <typename T>
33  T sqr(T x)
34  {
35  return x*x;
36  }
37 }
38 
39 
40 
41 using namespace std;
42 using namespace ctfseeding;
43 
45 
47 extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")), //extra window in ThirdHitPredictionFromCircle range (divide by R to get phi)
48 maxChi2(cfg.getParameter<edm::ParameterSet>("maxChi2")),
49 fitFastCircle(cfg.getParameter<bool>("fitFastCircle")),
50 fitFastCircleChi2Cut(cfg.getParameter<bool>("fitFastCircleChi2Cut")),
51 useBendingCorrection(cfg.getParameter<bool>("useBendingCorrection")),
52 caThetaCut(cfg.getParameter<double>("CAThetaCut")),
53 caPhiCut(cfg.getParameter<double>("CAPhiCut")),
54 caHardPtCut(cfg.getParameter<double>("CAHardPtCut")),
55 caOnlyOneLastHitPerLayerFilter(cfg.getParameter<bool>("CAOnlyOneLastHitPerLayerFilter"))
56 {
57  if(needSeedingLayerSetsHits)
59 
60  if (cfg.exists("SeedComparitorPSet"))
61  {
62  edm::ParameterSet comparitorPSet =
63  cfg.getParameter<edm::ParameterSet>("SeedComparitorPSet");
64  std::string comparitorName = comparitorPSet.getParameter<std::string>("ComponentName");
65  if (comparitorName != "none")
66  {
67  theComparitor.reset(SeedComparitorFactory::get()->create(comparitorName, comparitorPSet, iC));
68  }
69  }
70 }
71 
73 {
74 }
75 
77  desc.add<double>("extraHitRPhitolerance", 0.1);
78  desc.add<bool>("fitFastCircle", false);
79  desc.add<bool>("fitFastCircleChi2Cut", false);
80  desc.add<bool>("useBendingCorrection", false);
81  desc.add<double>("CAThetaCut", 0.00125);
82  desc.add<double>("CAPhiCut", 10);
83  desc.add<double>("CAHardPtCut", 0);
84  desc.add<bool>("CAOnlyOneLastHitPerLayerFilter",false);
85  edm::ParameterSetDescription descMaxChi2;
86  descMaxChi2.add<double>("pt1", 0.2);
87  descMaxChi2.add<double>("pt2", 1.5);
88  descMaxChi2.add<double>("value1", 500);
89  descMaxChi2.add<double>("value2", 50);
90  descMaxChi2.add<bool>("enabled", true);
91  desc.add<edm::ParameterSetDescription>("maxChi2", descMaxChi2);
92 
93  edm::ParameterSetDescription descComparitor;
94  descComparitor.add<std::string>("ComponentName", "none");
95  descComparitor.setAllowAnything(); // until we have moved SeedComparitor too to EDProducers
96  desc.add<edm::ParameterSetDescription>("SeedComparitorPSet", descComparitor);
97 }
98 
100  if (theComparitor) theComparitor->init(ev, es);
101 }
102 
103 namespace {
104  template <typename T_HitDoublets, typename T_GeneratorOrPairsFunction>
105  void fillGraph(const SeedingLayerSetsHits& layers, CAGraph& g, T_HitDoublets& hitDoublets,
106  T_GeneratorOrPairsFunction generatorOrPairsFunction) {
107  for (unsigned int i = 0; i < layers.size(); i++)
108  {
109  for (unsigned int j = 0; j < 4; ++j)
110  {
111  auto vertexIndex = 0;
112  auto foundVertex = std::find(g.theLayers.begin(), g.theLayers.end(),
113  layers[i][j].name());
114  if (foundVertex == g.theLayers.end())
115  {
116  g.theLayers.emplace_back(layers[i][j].name(), layers[i][j].hits().size());
117  vertexIndex = g.theLayers.size() - 1;
118  }
119  else
120  {
121  vertexIndex = foundVertex - g.theLayers.begin();
122  }
123  if (j == 0)
124  {
125 
126  if (std::find(g.theRootLayers.begin(), g.theRootLayers.end(),
127  vertexIndex) == g.theRootLayers.end())
128  {
129  g.theRootLayers.emplace_back(vertexIndex);
130 
131  }
132 
133  }
134  else
135  {
136 
137  auto innerVertex = std::find(g.theLayers.begin(),
138  g.theLayers.end(), layers[i][j - 1].name());
139 
140  CALayerPair tmpInnerLayerPair(innerVertex - g.theLayers.begin(),
141  vertexIndex);
142 
143  if (std::find(g.theLayerPairs.begin(), g.theLayerPairs.end(),
144  tmpInnerLayerPair) == g.theLayerPairs.end())
145  {
146  const bool nonEmpty = generatorOrPairsFunction(layers[i][j-1], layers[i][j], hitDoublets);
147  if(nonEmpty) {
148  g.theLayerPairs.push_back(tmpInnerLayerPair);
149  g.theLayers[vertexIndex].theInnerLayers.push_back(
150  innerVertex - g.theLayers.begin());
151  innerVertex->theOuterLayers.push_back(vertexIndex);
152  g.theLayers[vertexIndex].theInnerLayerPairs.push_back(
153  g.theLayerPairs.size() - 1);
154  innerVertex->theOuterLayerPairs.push_back(
155  g.theLayerPairs.size() - 1);
156  }
157  }
158 
159  }
160 
161  }
162  }
163  }
164 }
165 
166 
169  const edm::EventSetup& es)
170 {
172  ev.getByToken(theSeedingLayerToken, hlayers);
173  const SeedingLayerSetsHits& layers = *hlayers;
174  if (layers.numberOfLayersInSet() != 4)
175  throw cms::Exception("Configuration")
176  << "CAHitQuadrupletsGenerator expects SeedingLayerSetsHits::numberOfLayersInSet() to be 4, got "
177  << layers.numberOfLayersInSet();
178 
179  CAGraph g;
180 
181 
182  std::vector<HitDoublets> hitDoublets;
183 
184 
185  HitPairGeneratorFromLayerPair thePairGenerator(0, 1, &theLayerCache);
186  fillGraph(layers, g, hitDoublets,
189  std::vector<HitDoublets>& hitDoublets) {
190  hitDoublets.emplace_back(thePairGenerator.doublets(region, ev, es, inner, outer));
191  return true;
192  });
193 
194  if (theComparitor)
195  theComparitor->init(ev, es);
196 
197  std::vector<const HitDoublets *> hitDoubletsPtr;
198  hitDoubletsPtr.reserve(hitDoublets.size());
199  for(const auto& e: hitDoublets)
200  hitDoubletsPtr.emplace_back(&e);
201 
202  hitQuadruplets(region, result, hitDoubletsPtr, g, es);
204 }
205 
208  const edm::EventSetup& es,
209  const SeedingLayerSetsHits& layers) {
210  CAGraph g;
211 
212  std::vector<const HitDoublets *> hitDoublets;
213 
214  auto layerPairEqual = [](const IntermediateHitDoublets::LayerPairHitDoublets& pair,
217  return pair.innerLayerIndex() == inner && pair.outerLayerIndex() == outer;
218  };
219  fillGraph(layers, g, hitDoublets,
220  [&](const SeedingLayerSetsHits::SeedingLayer& inner,
222  std::vector<const HitDoublets *>& hitDoublets) {
223  using namespace std::placeholders;
224  auto found = std::find_if(regionLayerPairs.begin(), regionLayerPairs.end(),
225  std::bind(layerPairEqual, _1, inner.index(), outer.index()));
226  if(found != regionLayerPairs.end()) {
227  hitDoublets.emplace_back(&(found->doublets()));
228  return true;
229  }
230  return false;
231  });
232 
233  hitQuadruplets(regionLayerPairs.region(), result, hitDoublets, g, es);
234 }
235 
238  std::vector<const HitDoublets *>& hitDoublets, const CAGraph& g,
239  const edm::EventSetup& es) {
240  //Retrieve tracker topology from geometry
242  es.get<TrackerTopologyRcd>().get(tTopoHand);
243  const TrackerTopology *tTopo=tTopoHand.product();
244 
245  const int numberOfHitsInNtuplet = 4;
246  std::vector<CACell::CAntuplet> foundQuadruplets;
247 
248  CellularAutomaton ca(g);
249 
250  ca.createAndConnectCells(hitDoublets, region, caThetaCut,
252 
253  ca.evolve(numberOfHitsInNtuplet);
254 
255  ca.findNtuplets(foundQuadruplets, numberOfHitsInNtuplet);
256 
257 
258  const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es);
259 
260  // re-used thoughout, need to be vectors because of RZLine interface
261  std::array<float, 4> bc_r;
262  std::array<float, 4> bc_z;
263  std::array<float, 4> bc_errZ2;
264  std::array<GlobalPoint, 4> gps;
265  std::array<GlobalError, 4> ges;
266  std::array<bool, 4> barrels;
267  unsigned int fourthLayerId = 0;
268  unsigned int previousfourthLayerId = 0;
269  int subDetId = 0;
270  int previousSubDetId = 0;
271  unsigned int sideId = 0;
272  unsigned int previousSideId = 0;
273  std::array<unsigned int, 2> previousCellIds ={{0,0}};
274  bool isTheSameTriplet = false;
275  bool isTheSameFourthLayer = false;
276  bool hasAlreadyPushedACandidate = false;
277  float selectedChi2 = std::numeric_limits<float>::max();
278 
279 
280  unsigned int numberOfFoundQuadruplets = foundQuadruplets.size();
281 
282  // Loop over quadruplets
283  for (unsigned int quadId = 0; quadId < numberOfFoundQuadruplets; ++quadId)
284  {
285 
286  auto isBarrel = [](const unsigned id) -> bool
287  {
288  return id == PixelSubdetector::PixelBarrel;
289  };
290  for(unsigned int i = 0; i< 3; ++i)
291  {
292  auto const& ahit = foundQuadruplets[quadId][i]->getInnerHit();
293  gps[i] = ahit->globalPosition();
294  ges[i] = ahit->globalPositionError();
295  barrels[i] = isBarrel(ahit->geographicalId().subdetId());
296  }
297 
298  auto const& ahit = foundQuadruplets[quadId][2]->getOuterHit();
299  gps[3] = ahit->globalPosition();
300  ges[3] = ahit->globalPositionError();
301  barrels[3] = isBarrel(ahit->geographicalId().subdetId());
302 
304  {
305  fourthLayerId = tTopo->layer(ahit->geographicalId());
306  sideId = tTopo->side(ahit->geographicalId());
307  subDetId = ahit->geographicalId().subdetId();
308  isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0]->getCellId() == previousCellIds[0]) && (foundQuadruplets[quadId][1]->getCellId() == previousCellIds[1]);
309  isTheSameFourthLayer = (quadId != 0) && (fourthLayerId == previousfourthLayerId) && (subDetId == previousSubDetId) && (sideId == previousSideId);
310 
311  previousCellIds = {{foundQuadruplets[quadId][0]->getCellId(), foundQuadruplets[quadId][1]->getCellId()}};
312  previousfourthLayerId = fourthLayerId;
313 
314 
315  if(!(isTheSameTriplet && isTheSameFourthLayer ))
316  {
317  selectedChi2 = std::numeric_limits<float>::max();
318  hasAlreadyPushedACandidate = false;
319  }
320 
321  }
322 
323 
324  // TODO:
325  // - 'line' is not used for anything
326  // - if we decide to always do the circle fit for 4 hits, we don't
327  // need ThirdHitPredictionFromCircle for the curvature; then we
328  // could remove extraHitRPhitolerance configuration parameter
329  PixelRecoLineRZ line(gps[0], gps[2]);
330  ThirdHitPredictionFromCircle predictionRPhi(gps[0], gps[2], extraHitRPhitolerance);
331  const float curvature = predictionRPhi.curvature(ThirdHitPredictionFromCircle::Vector2D(gps[1].x(), gps[1].y()));
332  const float abscurv = std::abs(curvature);
333  const float thisMaxChi2 = maxChi2Eval.value(abscurv);
334 
335  if (theComparitor)
336  {
337  SeedingHitSet tmpTriplet(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit());
338 
339 
340  if (!theComparitor->compatible(tmpTriplet) )
341  {
342  continue;
343  }
344  }
345 
346  float chi2 = std::numeric_limits<float>::quiet_NaN();
347  // TODO: Do we have any use case to not use bending correction?
349  {
350  // Following PixelFitterByConformalMappingAndLine
351  const float simpleCot = ( gps.back().z() - gps.front().z() ) / (gps.back().perp() - gps.front().perp() );
352  const float pt = 1.f / PixelRecoUtilities::inversePt(abscurv, es);
353  for (int i=0; i < 4; ++i)
354  {
355  const GlobalPoint & point = gps[i];
356  const GlobalError & error = ges[i];
357  bc_r[i] = sqrt( sqr(point.x() - region.origin().x()) + sqr(point.y() - region.origin().y()) );
358  bc_r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt, es)(bc_r[i]);
359  bc_z[i] = point.z() - region.origin().z();
360  bc_errZ2[i] = (barrels[i]) ? error.czz() : error.rerr(point)*sqr(simpleCot);
361  }
362  RZLine rzLine(bc_r, bc_z, bc_errZ2, RZLine::ErrZ2_tag());
363  chi2 = rzLine.chi2();
364  }
365  else
366  {
367  RZLine rzLine(gps, ges, barrels);
368  chi2 = rzLine.chi2();
369  }
370  if (edm::isNotFinite(chi2) || chi2 > thisMaxChi2)
371  {
372  continue;
373 
374  }
375  // TODO: Do we have any use case to not use circle fit? Maybe
376  // HLT where low-pT inefficiency is not a problem?
377  if (fitFastCircle)
378  {
379  FastCircleFit c(gps, ges);
380  chi2 += c.chi2();
381  if (edm::isNotFinite(chi2))
382  continue;
383  if (fitFastCircleChi2Cut && chi2 > thisMaxChi2)
384  continue;
385  }
386 
388  {
389  if (chi2 < selectedChi2)
390  {
391  selectedChi2 = chi2;
392 
393  if(hasAlreadyPushedACandidate)
394  {
395  result.pop_back();
396 
397  }
398  result.emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(),
399  foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit());
400  hasAlreadyPushedACandidate = true;
401 
402  }
403  }
404  else
405  {
406 
407  result.emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit());
408  }
409 
410  }
411 
412 }
413 
414 
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
unsigned short numberOfLayersInSet() const
Get number of layers in each SeedingLayerSets.
void findNtuplets(std::vector< CACell::CAntuplet > &, const unsigned int)
tuple cfg
Definition: looper.py:293
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
Definition: LayerTriplets.cc:4
void hitNtuplets(const IntermediateHitDoublets::RegionLayerSets &regionLayerPairs, OrderedHitSeeds &result, const edm::EventSetup &es, const SeedingLayerSetsHits &layers)
const_iterator begin() const
void createAndConnectCells(const std::vector< const HitDoublets * > &, const TrackingRegion &, const float, const float, const float)
GlobalPoint const & origin() const
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
bool isBarrel(GeomDetEnumerators::SubDetector m)
CAHitQuadrupletGenerator(const edm::ParameterSet &cfg, edm::ConsumesCollector &&iC, bool needSeedingLayerSetsHits=true)
const_iterator end() const
unsigned int side(const DetId &id) const
T y() const
Definition: PV3DBase.h:63
bool exists(std::string const &parameterName) const
checks if a parameter exists
bool ev
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:7
std::unique_ptr< SeedComparitor > theComparitor
The Signals That Services Can Subscribe To This is based on ActivityRegistry and is current per Services can connect to the signals distributed by the ActivityRegistry in order to monitor the activity of the application Each possible callback has some defined which we here list in angle e g
Definition: Activities.doc:4
#define constexpr
tuple result
Definition: mps_fire.py:84
Range curvature(double transverseIP) const
virtual void hitQuadruplets(const TrackingRegion &reg, OrderedHitSeeds &quadruplets, const edm::Event &ev, const edm::EventSetup &es)
from base class
T inversePt(T curvature, const edm::EventSetup &iSetup)
T x() const
Cartesian x coordinate.
bool isNotFinite(T x)
Definition: isFinite.h:10
T curvature(T InversePt, const edm::EventSetup &iSetup)
SeedingLayerSetsHits::LayerIndex innerLayerIndex() const
T sqrt(T t)
Definition: SSEVec.h:18
std::vector< CALayer > theLayers
Definition: CAGraph.h:67
T z() const
Definition: PV3DBase.h:64
void initEvent(const edm::Event &ev, const edm::EventSetup &es)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
ParameterDescriptionBase * add(U const &iLabel, T const &value)
std::vector< CALayerPair > theLayerPairs
Definition: CAGraph.h:68
HitDoublets doublets(const TrackingRegion &reg, const edm::Event &ev, const edm::EventSetup &es, Layers layers)
Definition: RZLine.h:12
const TrackingRegion & region() const
const QuantityDependsPt maxChi2
float chi2() const
Definition: FastCircleFit.h:66
T rerr(const GlobalPoint &aPoint) const
const T & get() const
Definition: EventSetup.h:56
T const * product() const
Definition: ESHandle.h:86
unsigned int layer(const DetId &id) const
SeedingLayerSetsHits::LayerIndex outerLayerIndex() const
Square< F >::type sqr(const F &f)
Definition: Square.h:13
std::vector< int > theRootLayers
Definition: CAGraph.h:69
static void fillDescriptions(edm::ParameterSetDescription &desc)
QuantityDependsPtEval evaluator(const edm::EventSetup &es) const
unsigned int subDetId[20]
long double T
T x() const
Definition: PV3DBase.h:62
float chi2() const
Definition: RZLine.h:97
unsigned short size() const
Get the number of SeedingLayerSets.
tuple size
Write out results.
edm::EDGetTokenT< SeedingLayerSetsHits > theSeedingLayerToken
T get(const Candidate &c)
Definition: component.h:55
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
void evolve(const unsigned int)