CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TrajectorySeedProducer.cc
Go to the documentation of this file.
2 
7 
10 
13 
17 
22 
26 
29 
32 
33 //Propagator withMaterial
35 
36 
37 #include <unordered_set>
38 
45 
46 template class SeedingTree<TrackingLayer>;
47 template class SeedingNode<TrackingLayer>;
48 
51  magneticFieldMap(nullptr),
52  trackerGeometry(nullptr),
53  trackerTopology(nullptr),
54  theRegionProducer(nullptr)
55 {
56 
57  produces<TrajectorySeedCollection>();
58 
59  // tokens
60  recHitCombinationsToken = consumes<FastTrackerRecHitCombinationCollection>(conf.getParameter<edm::InputTag>("recHitCombinations"));
61  measurementTrackerEventToken = consumes<MeasurementTrackerEvent>(conf.getParameter<edm::InputTag>("MeasurementTrackerEvent"));
62  hitMasks_exists = conf.exists("hitMasks");
63  if (hitMasks_exists){
64  edm::InputTag hitMasksTag = conf.getParameter<edm::InputTag>("hitMasks");
65  hitMasksToken = consumes<std::vector<bool> >(hitMasksTag);
66  }
67 
68  // The smallest number of hits for a track candidate
69  minLayersCrossed = conf.getParameter<unsigned int>("minLayersCrossed");
70 
71  // read Layers
72  std::vector<std::string> layerStringList = conf.getParameter<std::vector<std::string>>("layerList");
73  for(auto it=layerStringList.cbegin(); it < layerStringList.cend(); ++it)
74  {
75  std::vector<TrackingLayer> trackingLayerList;
76  std::string line = *it;
78  while (pos != std::string::npos)
79  {
80  pos=line.find("+");
81  std::string layer = line.substr(0, pos);
83 
84  trackingLayerList.push_back(layerSpec);
85  line=line.substr(pos+1,std::string::npos);
86  }
87  _seedingTree.insert(trackingLayerList);
88  seedingLayers.push_back(std::move(trackingLayerList));
89  }
90 
91  // region producer
92  if(conf.exists("RegionFactoryPSet")){
93  edm::ParameterSet regfactoryPSet = conf.getParameter<edm::ParameterSet>("RegionFactoryPSet");
94  std::string regfactoryName = regfactoryPSet.getParameter<std::string>("ComponentName");
95  theRegionProducer.reset(TrackingRegionProducerFactory::get()->create(regfactoryName,regfactoryPSet, consumesCollector()));
96 
97  // seed creator
98  const edm::ParameterSet & seedCreatorPSet = conf.getParameter<edm::ParameterSet>("SeedCreatorPSet");
99  std::string seedCreatorName = seedCreatorPSet.getParameter<std::string>("ComponentName");
100  seedCreator.reset(SeedCreatorFactory::get()->create( seedCreatorName, seedCreatorPSet));
101  }
102 }
103 
104 bool
106 {
107 
108  const DetLayer * innerLayer =
110  const DetLayer * outerLayer =
112 
114 
115  for(Regions::const_iterator ir=regions.begin(); ir < regions.end(); ++ir){
116 
117  auto const & gs = outerHit.hit()->globalState();
118  auto loc = gs.position-(*ir)->origin().basicVector();
119  const HitRZCompatibility * checkRZ = (*ir)->checkRZ(innerLayer, outerHit.hit(), *es_, outerLayer,
120  loc.perp(),gs.position.z(),gs.errorR,gs.errorZ);
121 
122  float u = innerLayer->isBarrel() ? loc.perp() : gs.position.z();
123  float v = innerLayer->isBarrel() ? gs.position.z() : loc.perp();
124  float dv = innerLayer->isBarrel() ? gs.errorZ : gs.errorR;
125  constexpr float nSigmaRZ = 3.46410161514f;
126  Range allowed = checkRZ->range(u);
127  float vErr = nSigmaRZ * dv;
128  Range hitRZ(v-vErr, v+vErr);
129  Range crossRange = allowed.intersection(hitRZ);
130 
131  if( ! crossRange.empty()){
132  seedCreator->init(**ir,*es_,0);
133  return true;}
134 
135  }
136  return false;
137 }
138 
140  const std::vector<TrajectorySeedHitCandidate>& trackerRecHits,
141  std::vector<int>& hitIndicesInTree,
142  const SeedingNode<TrackingLayer>* node, unsigned int trackerHit
143  ) const
144 {
145  if (!node->getParent() || hitIndicesInTree[node->getParent()->getIndex()]>=0)
146  {
147  if (hitIndicesInTree[node->getIndex()]<0)
148  {
149  const TrajectorySeedHitCandidate& currentTrackerHit = trackerRecHits[trackerHit];
150  if (!isHitOnLayer(currentTrackerHit,node->getData()))
151  {
152  return nullptr;
153  }
154  if (!passHitTuplesCuts(*node,trackerRecHits,hitIndicesInTree,currentTrackerHit))
155  {
156  return nullptr;
157  }
158  hitIndicesInTree[node->getIndex()]=trackerHit;
159  if (node->getChildrenSize()==0)
160  {
161  return node;
162  }
163  return nullptr;
164  }
165  else
166  {
167  for (unsigned int ichild = 0; ichild<node->getChildrenSize(); ++ichild)
168  {
169  const SeedingNode<TrackingLayer>* seed = insertHit(trackerRecHits,hitIndicesInTree,node->getChild(ichild),trackerHit);
170  if (seed)
171  {
172  return seed;
173  }
174  }
175  }
176  }
177  return nullptr;
178 }
179 
180 
181 std::vector<unsigned int> TrajectorySeedProducer::iterateHits(
182  unsigned int start,
183  const std::vector<TrajectorySeedHitCandidate>& trackerRecHits,
184  std::vector<int> hitIndicesInTree,
185  bool processSkippedHits
186  ) const
187 {
188  for (unsigned int irecHit = start; irecHit<trackerRecHits.size(); ++irecHit)
189  {
190  unsigned int currentHitIndex=irecHit;
191 
192  for (unsigned int inext=currentHitIndex+1; inext< trackerRecHits.size(); ++inext)
193  {
194  //if multiple hits are on the same layer -> follow all possibilities by recusion
195  if (trackerRecHits[currentHitIndex].getTrackingLayer()==trackerRecHits[inext].getTrackingLayer())
196  {
197  if (processSkippedHits)
198  {
199  //recusively call the method again with hit 'inext' but skip all following on the same layer though 'processSkippedHits=false'
200  std::vector<unsigned int> seedHits = iterateHits(
201  inext,
202  trackerRecHits,
203  hitIndicesInTree,
204  false
205  );
206  if (seedHits.size()>0)
207  {
208  return seedHits;
209  }
210  }
211  irecHit+=1;
212  }
213  else
214  {
215  break;
216  }
217  }
218 
219  processSkippedHits=true;
220 
221  const SeedingNode<TrackingLayer>* seedNode = nullptr;
222  for (unsigned int iroot=0; seedNode==nullptr && iroot<_seedingTree.numberOfRoots(); ++iroot)
223  {
224  seedNode=insertHit(trackerRecHits,hitIndicesInTree,_seedingTree.getRoot(iroot), currentHitIndex);
225  }
226  if (seedNode)
227  {
228  std::vector<unsigned int> seedIndices(seedNode->getDepth()+1);
229  while (seedNode)
230  {
231  seedIndices[seedNode->getDepth()]=hitIndicesInTree[seedNode->getIndex()];
232  seedNode=seedNode->getParent();
233  }
234  return seedIndices;
235  }
236 
237  }
238 
239  return std::vector<unsigned int>();
240 
241 }
242 
243 void
245 {
246 
247  // services
248  edm::ESHandle<MagneticField> magneticFieldHandle;
249  edm::ESHandle<TrackerGeometry> trackerGeometryHandle;
250  edm::ESHandle<MagneticFieldMap> magneticFieldMapHandle;
251  edm::ESHandle<TrackerTopology> trackerTopologyHandle;
252 
253  es.get<IdealMagneticFieldRecord>().get(magneticFieldHandle);
254  es.get<TrackerDigiGeometryRecord>().get(trackerGeometryHandle);
255  es.get<MagneticFieldMapRecord>().get(magneticFieldMapHandle);
256  es.get<TrackerTopologyRcd>().get(trackerTopologyHandle);
257 
258  magneticField = &(*magneticFieldHandle);
259  trackerGeometry = &(*trackerGeometryHandle);
260  magneticFieldMap = &(*magneticFieldMapHandle);
261  trackerTopology = &(*trackerTopologyHandle);
262 
263  es_ = &es;
264 
265  // hit masks
266  std::unique_ptr<HitMaskHelper> hitMaskHelper;
267  if (hitMasks_exists){
269  e.getByToken(hitMasksToken,hitMasks);
270  hitMaskHelper.reset(new HitMaskHelper(hitMasks.product()));
271  }
272 
273  // hit combinations
275  e.getByToken(recHitCombinationsToken, recHitCombinations);
276 
277  // measurement tracker event (only used to access services)
278  edm::Handle<MeasurementTrackerEvent> measurementTrackerEventHandle;
279  e.getByToken(measurementTrackerEventToken,measurementTrackerEventHandle);
280  measurementTrackerEvent = measurementTrackerEventHandle.product();
281 
282  // output
283  std::unique_ptr<TrajectorySeedCollection> output(new TrajectorySeedCollection());
284 
285  // produce the regions;
286  if(!theRegionProducer){
287  edm::LogWarning("TrajectorySeedProducer") << " RegionFactory is not initialised properly, please check your configuration. Producing empty seed collection" << std::endl;
288  e.put(std::move(output));
289  return;
290  }
291 
292  regions = theRegionProducer->regions(e,es);
293 
294 
295  for ( unsigned icomb=0; icomb<recHitCombinations->size(); ++icomb)
296  {
297 
298 
299  FastTrackerRecHitCombination recHitCombination = (*recHitCombinations)[icomb];
300 
301  TrajectorySeedHitCandidate previousTrackerHit;
302  TrajectorySeedHitCandidate currentTrackerHit;
303  unsigned int layersCrossed=0;
304 
305 
306  std::vector<TrajectorySeedHitCandidate> trackerRecHits;
307  for (const auto & _hit : recHitCombination )
308  {
309  // skip masked hits
310  if(hitMaskHelper && hitMaskHelper->mask(_hit.get()))
311  continue;
312 
313  previousTrackerHit=currentTrackerHit;
314 
315  currentTrackerHit = TrajectorySeedHitCandidate(_hit.get(),trackerGeometry,trackerTopology);
316 
317  if (!currentTrackerHit.isOnTheSameLayer(previousTrackerHit))
318  {
319  ++layersCrossed;
320  }
321  if (_seedingTree.getSingleSet().find(currentTrackerHit.getTrackingLayer())!=_seedingTree.getSingleSet().end())
322  {
323  //add only the hits which are actually on the requested layers
324  trackerRecHits.push_back(std::move(currentTrackerHit));
325  }
326  }
327 
328  if ( layersCrossed < minLayersCrossed)
329  {
330  continue;
331  }
332 
333  // set the combination index
334 
335  //A SeedingNode is associated by its index to this list. The list stores the indices of the hits in 'trackerRecHits'
336  /* example
337  SeedingNode | hit index | hit
338  -------------------------------------------------------------------------------
339  index= 0: [BPix1] | hitIndicesInTree[0] (=1) | trackerRecHits[1]
340  index= 1: -- [BPix2] | hitIndicesInTree[1] (=3) | trackerRecHits[3]
341  index= 2: -- -- [BPix3] | hitIndicesInTree[2] (=4) | trackerRecHits[4]
342  index= 3: -- -- [FPix1_pos] | hitIndicesInTree[3] (=6) | trackerRecHits[6]
343  index= 4: -- -- [FPix1_neg] | hitIndicesInTree[4] (=7) | trackerRecHits[7]
344 
345  The implementation has been chosen such that the tree only needs to be build once upon construction.
346  */
347 
348  // find the first combination of hits,
349  // compatible with the seedinglayer,
350  // and with one of the tracking regions
351  std::vector<int> hitIndicesInTree(_seedingTree.numberOfNodes(),-1);
352  std::vector<unsigned int> seedHitNumbers = iterateHits(0,trackerRecHits,hitIndicesInTree,true);
353 
354  // create a seed from those hits
355  if (seedHitNumbers.size()>1){
356  // temporary hit copies to allow setting the recHitCombinationIndex
358  for(unsigned iIndex = 0;iIndex < seedHitNumbers.size();++iIndex){
359  seedHits.push_back(trackerRecHits[seedHitNumbers[iIndex]].hit()->clone());
360  }
362 
363  // the actual seed creation
364  seedCreator->makeSeed(*output,SeedingHitSet(&seedHits[0],&seedHits[1],
365  seedHits.size() >=3 ? &seedHits[2] : nullptr,
366  seedHits.size() >=4 ? &seedHits[3] : nullptr));
367  }
368  }
369  e.put(std::move(output));
370 
371 }
bool passHitTuplesCuts(const SeedingNode< TrackingLayer > &seedingNode, const std::vector< TrajectorySeedHitCandidate > &trackerRecHits, const std::vector< int > &hitIndicesInTree, const TrajectorySeedHitCandidate &currentTrackerHit) const
method checks if a SimTrack fulfills the quality requirements.
T getParameter(std::string const &) const
tuple start
Check for commandline option errors.
Definition: dqm_diff.py:58
const MagneticField * magneticField
TrajectorySeedProducer(const edm::ParameterSet &conf)
TrackingRecHitGlobalState globalState() const
std::unique_ptr< TrackingRegionProducer > theRegionProducer
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:464
edm::EDGetTokenT< std::vector< bool > > hitMasksToken
virtual Range range(const float &rORz) const =0
const SingleSet & getSingleSet() const
Definition: SeedingTree.h:204
size_type size() const
Definition: OwnVector.h:254
bool exists(std::string const &parameterName) const
checks if a parameter exists
virtual void produce(edm::Event &e, const edm::EventSetup &es)
const DATA & getData() const
Definition: SeedingTree.h:136
#define nullptr
unsigned int getChildrenSize() const
Definition: SeedingTree.h:121
uint16_t size_type
#define constexpr
tuple node
Definition: Node.py:50
unsigned int getDepth() const
Definition: SeedingTree.h:84
const SeedingNode * getParent() const
Definition: SeedingTree.h:112
bool insert(const std::vector< DATA > &dataList)
Definition: SeedingTree.h:179
unsigned int numberOfNodes() const
Definition: SeedingTree.h:227
void push_back(D *&d)
Definition: OwnVector.h:280
unsigned int numberOfRoots() const
Definition: SeedingTree.h:222
const TrackerTopology * trackerTopology
std::vector< TrajectorySeed > TrajectorySeedCollection
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:120
const GeomDet * det() const
const edm::EventSetup * es_
ConsumesCollector consumesCollector()
Use a ConsumesCollector to gather consumes information from helper functions.
T z() const
Definition: PV3DBase.h:64
def move
Definition: eostools.py:510
virtual std::vector< unsigned int > iterateHits(unsigned int start, const std::vector< TrajectorySeedHitCandidate > &trackerRecHits, std::vector< int > hitIndicesInTree, bool processSkippedHits) const
method tries to insert all hits into the tree structure.
const DetLayer * detLayer(const DetId &id) const
obsolete method. Use idToLayer() instead.
DetId geographicalId() const
The label of this GeomDet.
Definition: GeomDet.h:77
static TrackingLayer createFromString(std::string layerSpecification)
bool isHitOnLayer(const TrajectorySeedHitCandidate &trackerRecHit, const TrackingLayer &layer) const
const MeasurementTracker & measurementTracker() const
tuple conf
Definition: dbtoconf.py:185
bool pass2HitsCuts(const TrajectorySeedHitCandidate &hit1, const TrajectorySeedHitCandidate &hit2) const
bool isBarrel() const
Definition: DetLayer.h:32
std::vector< FastTrackerRecHitRef > FastTrackerRecHitCombination
PixelRecoRange< float > Range
T const * product() const
Definition: Handle.h:81
virtual const Surface::PositionType & position() const
Returns position of the surface.
const T & get() const
Definition: EventSetup.h:56
std::vector< std::vector< TrackingLayer > > seedingLayers
TEveGeoShape * clone(const TEveElement *element, TEveElement *parent)
Definition: eve_macros.cc:135
edm::EDGetTokenT< FastTrackerRecHitCombinationCollection > recHitCombinationsToken
const MeasurementTrackerEvent * measurementTrackerEvent
const SeedingNode< DATA > * getRoot(unsigned int i) const
Definition: SeedingTree.h:232
void setRecHitCombinationIndex(edm::OwnVector< T > &recHits, int32_t icomb)
edm::EDGetTokenT< MeasurementTrackerEvent > measurementTrackerEventToken
const TrackerGeometry * trackerGeometry
const SeedingNode< TrackingLayer > * insertHit(const std::vector< TrajectorySeedHitCandidate > &trackerRecHits, std::vector< int > &hitIndicesInTree, const SeedingNode< TrackingLayer > *node, unsigned int trackerHit) const
const FastTrackerRecHit * hit() const
The Hit itself.
unsigned int getIndex() const
Definition: SeedingTree.h:107
const TrackingLayer & getTrackingLayer() const
SeedingTree< TrackingLayer > _seedingTree
bool isOnTheSameLayer(const TrajectorySeedHitCandidate &other) const
Check if two hits are on the same layer of the same subdetector.
const GeometricSearchTracker * geometricSearchTracker() const
std::unique_ptr< SeedCreator > seedCreator
SurfaceDeformation * create(int type, const std::vector< double > &params)
const MagneticFieldMap * magneticFieldMap
T get(const Candidate &c)
Definition: component.h:55
const SeedingNode< DATA > * getChild(unsigned int ichild) const
Definition: SeedingTree.h:126